In what follows, the argument of the functions $G_i$ is implied to be $z$ in the formulas for $\tilde{a}_{gg,Q}^{(3),+} (z)$ and $\tilde{a}_{gg,Q}^{(3),\text{reg}}(z)$, and it is implied to be $y$ in the functions $\Phi_i$ which follow. Such arguments are omitted in the interest of brevity. Defining for the inverse Mellin transform of $\tilde{a}_{gg,Q}^{(3)}$

\begin{eqnarray}
\tilde{a}_{gg,Q}^{(3)} (N) &=& \tilde{a}_{gg,Q}^{(3),\delta}  + \int_{0}^{1} \text{d} z \ \left( z^{N-1} - 1 \right) \ \tilde{a}_{gg,Q}^{(3),+} (z)
+ \int_{0}^{1} \text{d} z \ z^{N-1} \ \tilde{a}_{gg,Q}^{(3),\text{reg}} (z) ~,
\end{eqnarray}

we find

with the polynomials

The functions $\Phi_1,\ldots,\Phi_8$, which appear as arguments of a further integral, are:

The polynomials $R_i$ are:

In the Appendices of [198], a set of relations and evaluations of iterated integrals can be found, which were used for this calculation.  We present further relations in Appendices <ref> and <ref>.

§.§.§ Numerical results

A numerical evaluation of the size of this correction compared to the total $\mathcal{O}(T_F^2)$ is shown in Figure <ref>, depicting the ratio between the two-mass contribution and the total contribution to $O(T_F^2)$. Their appreciable size of $O(0.4)$ makes it necessary to include them in numerical studies.

§.§.§ Summary

We computed the two-mass contributions to the polarized OMEs $A_{gg,Q}^{(3)}$ and $A_{Qq}^{(3),PS}$ in the Larin scheme at $\mathcal O(a_s^3)$ in semi-analytic form in $z$-space and, for $A_{gg,Q}^{(3)}$, analytically in $N$-space. We show that their numerical contribution is of the same order of magnitude as the single-mass contribution and conclude that two-mass effects should be taken into account in the variable flavour number scheme.

The ratio of the two mass contribution to the total
contribution at $O(T_F^2)$ for the
polarized massive OME $A_{gg,Q}^{(3)}$ as a function of the momentum
fraction $z$ and the virtuality $\mu^2$.
Dashed line: $\mu^2 = 50~\GeV^2$;
Dash–dotted line: $\mu^2 = 100~\GeV^2$;
Full line: $\mu^2 = 1000~\GeV^2$.
For the values of $m_c$ and $m_b$ we refer to the on–shell heavy
quark masses
$m_c = 1.59~\GeV$ and $m_b = 4.78~\GeV$ [286, 287], from [198].

§ THE LOGARITHMIC SINGLE-MASS CONTRIBUTIONS TO THE POLARIZED ASYMPTOTIC $\MATHCAL{O}(A_S^3)$ WILSON COEFFICIENTS IN DEEPLY INELASTIC SCATTERING

Heavy flavour contributions to the structure functions in deep-inelastic scattering need to be taken into account in precision studies of DIS. They arise when quark masses are not neglected in the DIS process. These heavy flavour contributions can be attributed to the Wilson coefficients

\begin{eqnarray}
\mathcal{C}_i^{S,PS,NS}\Big(N,N_F+1,\frac{Q^2}{\mu^2},\frac{m^2}{\mu^2}\Big)
&=&C_i^{S,PS,NS}\Big(N,N_F,\frac{Q^2}{\mu^2} \Big) + H_i^{S,PS}(N,N_F+1,\frac{Q^2}{\mu^2},\frac{m^2}{\mu^2}\Big)
\nonumber\\&&
+ L_i^{S,PS,NS}(N,N_F+1,\frac{Q^2}{\mu^2},\frac{m^2}{\mu^2}\Big) ~,
\end{eqnarray}

where the symbol $N_F+1$ refers to $N_F$ massless flavours and one additional massive flavour [159, 158]. In this notation, the Wilson coefficients $H_i$ identify the contributions where the off-shell photon couples to a heavy quark, and $L_i$ to those where it couples to a light quark. The light flavour Wilson coefficients are denoted by $C_i$.

In the asymptotic limit $Q^2\gg m^2$, where $m$ is the mass of the heavy quark, the asymptotic form of the Wilson coefficients can be obtained through the factorization formula

\begin{eqnarray}
\mathcal{C}_j^{S,PS,NS,asymp} \Big(N,N_F+1,\frac{Q^2}{\mu^2},\frac{m^2}{\mu^2}\Big)
\nonumber\\&&
\hspace{-5cm}
\sum_i A_{ij}^{S,PS,NS}\Big(N,N_F+1,\frac{m^2}{\mu^2}\Big) C_i^{S,PS,NS}\Big(N,N_F+1,\frac{Q^2}{\mu^2}\Big) +\mathcal{O}\Big(\frac{m^2}{Q^2}\Big)
\end{eqnarray}

with $A_{i,j}^{S,PS,NS}$ the matrix elements of twist-2 operators between partonic states, which contain the contributions of heavy quarks in the loops:

\begin{equation}
A_{i,j}^{S,PS,NS}\Big(N,N_F+1,\frac{m^2}{\mu^2}\Big) = \langle j| O_i^{S,NS} |j \rangle, \qquad j=q,g.
\end{equation}

Both the light Wilson coefficients $C_i$ and the OMEs $A_{ij}$ can be calculated perturbatively in a series in $a_s$,

\begin{eqnarray}
C_i\Big(N,\frac{Q^2}{\mu^2}\Big) &=& \sum_{k=0}^\infty a_s^k \ C_i^{(k)}\Big(N,\frac{Q^2}{\mu^2}\Big) ~,\\
A_{ij}^{S,PS,NS}\Big(N,\frac{m^2}{\mu^2}\Big) &=& \sum_{k=0}^\infty a_s^k \ A_{ij}^{(k),S,PS,NS}\Big(N,\frac{m^2}{\mu^2}\Big) ~.
\end{eqnarray}

After the renormalization procedure is carried out, in our case in the $\MS$ scheme for the strong coupling constant and the on-mass-shell scheme for the mass of the heavy quark, the light Wilson coefficients and the OMEs take the form

\begin{eqnarray}
A_{ij}^{(k)} &=& \sum_{\ell=0}^k A_{ij}^{(k,\ell)} \ln^\ell\frac{m^2}{\mu^2} ~,\\
C_i^{(k)} &=& \sum_{\ell=0}^k C_i^{(k,\ell)} \ln^\ell\frac{Q^2}{\mu^2} ~.
\label{eq:logexpansion}
\end{eqnarray}

Treating $\gamma_5$, we work in the Larin scheme.

With the definitions

\begin{eqnarray}
\tilde f(N_f) &\equiv & \frac{f(N_F)}{N_F} ~,\\
\hat f(N_F) &\equiv & f(N_F+1)-f(N_F) ~,
\end{eqnarray}

the heavy quark Wilson coefficients take the asymptotic form [158]

\begin{eqnarray}
\label{eqWIL1}
L_{q,g_1}^{\sf NS,Q}(N_F+1) &=&
a_s^2 \Bigl[A_{qq,Q}^{(2), {\sf NS}}(N_F+1) +
\hat{C}^{(2), {\sf NS}}_{q,g_1}(N_F)\Bigr]
\N\\
a_s^3 \Bigl[A_{qq,Q}^{(3), {\sf NS}}(N_F+1)
+  A_{qq,Q}^{(2), {\sf NS}}(N_F+1) C_{q,g_1}^{(1), {\sf NS}}(N_F+1)
\N \\
&& \hspace*{5mm}
+ \hat{C}^{(3), {\sf NS}}_{q,g_1}(N_F)\Bigr]~,  \\
\label{eqWIL2}
L_{q,g_1}^{\sf PS}(N_F+1) &=&
a_s^3 \Bigl[~A_{qq,Q}^{(3), {\sf PS}}(N_F+1)
+  A_{gq,Q}^{(2)}(N_F+1)~~N_F\Ctil_{g,g_1}^{(1)}(N_F+1) \N \\
&& \hspace*{5mm}
+ N_F \hat{\Ctil}^{(3), {\sf PS}}_{q,g_1}(N_F)\Bigr]~,
\\
\label{eqWIL3}
L_{g,g_1}^{\sf S}(N_F+1) &=&
a_s^2 A_{gg,Q}^{(1)}(N_F+1)N_F \Ctil_{g,g_1}^{(1)}(N_F+1)
\N\\ &+&
a_s^3 \Bigl[~A_{qg,Q}^{(3)}(N_F+1)
+  A_{gg,Q}^{(1)}(N_F+1)~~N_F\Ctil_{g,g_1}^{(2)}(N_F+1)
\N\\ && \hspace*{5mm}
+  A_{gg,Q}^{(2)}(N_F+1)~~N_F\Ctil_{g,g_1}^{(1)}(N_F+1)
\N\\ && \hspace*{5mm}
+  ~A^{(1)}_{Qg}(N_F+1)~~N_F\Ctil_{q,g_1}^{(2), {\sf PS}}(N_F+1)
+ N_F \hat{\Ctil}^{(3)}_{g,g_1}(N_F)\Bigr]~,
\\ \N \\
\label{eq:WILPS}
H_{q,g_1}^{\sf PS}(N_F+1)
&=& a_s^2 \Bigl[~A_{Qq}^{(2), {\sf PS}}(N_F+1)
+~\Ctil_{q,g_1}^{(2), {\sf PS}}(N_F+1)\Bigr]
\\
&+& a_s^3 \Bigl[~A_{Qq}^{(3), {\sf PS}}(N_F+1)
+~\Ctil_{q,g_1}^{(3), {\sf PS}}(N_F+1) \N\\
&& %\hspace*{-20mm}
+ A_{gq,Q}^{(2)}(N_F+1)~\Ctil_{g,g_1}^{(1)}(N_F+1)
+ A_{Qq}^{(2), {\sf PS}}(N_F+1)~C_{q,g_1}^{(1), {\sf NS}}(N_F+1)
\Bigr]~,       \label{eqWIL4}
\N\\
\label{eq:WILS}
H_{g,g_1}^{\sf S}(N_F+1) &=& a_s \Bigl[~A_{Qg}^{(1)}(N_F+1)
+~\Ctil^{(1)}_{g,g_1}(N_F+1) \Bigr] \N\\
&+& a_s^2 \Bigl[~A_{Qg}^{(2)}(N_F+1)
+~A_{Qg}^{(1)}(N_F+1)~C^{(1), {\sf NS}}_{q,g_1}(N_F+1)\N\\ &&
\hspace*{5mm}
+~\Ctil^{(2)}_{g,g_1}(N_F+1) \Bigr]
\N\\ &+&
a_s^3 \Bigl[~A_{Qg}^{(3)}(N_F+1)
+~A_{Qg}^{(2)}(N_F+1)~C^{(1), {\sf NS}}_{q,g_1}(N_F+1)
\N\\ &&
\hspace*{5mm}
\N\\ && \hspace*{5mm}
C^{(2), {\sf NS}}_{q,g_1}(N_F+1)
+~\Ctil^{(2), {\sf PS}}_{q,g_1}(N_F+1)\Bigr\}
\N\\ && \hspace*{5mm}
+~\Ctil^{(3)}_{g,g_1}(N_F+1) \Bigr] ~.
\label{eqWIL5}
\end{eqnarray}

In the unpolarized case, the asymptotic form of the heavy flavour Wilson coefficients has been computed with this method to $\mathcal O(a^3)$ in [183, 165] and in the case of the structure function $g_1$ the non-singlet Wilson coefficient $L_{q,g_1}^{NS}$ was calculated in [184]; see also [290].

We briefly review how the logarithmic terms in the massless Wilson coefficients are obtained following [291]. From the renormalization group equation

\begin{equation}
\sum_i\Big[\frac{\partial}{\partial \ln\mu^2} +\beta(a_s) \frac{\partial}{\partial a_s} -\gamma_{ij} \Big] C_i\Big(a_s,\frac{Q^2}{\mu^2}\Big) = 0
\label{eq:rge}
\end{equation}

\begin{equation}
\beta(a_s) = - \sum_{k=0}^\infty a_s^{k+2} \beta_k ,
\end{equation}

inserting the ansatz (<ref>) into (<ref>) one obtains the explicit form of $C_i^{(k,\ell)}$ for $\ell>0$ in terms of lower-order quantities. For instance, one can write for the non-singlet Wilson coefficient in $N$-space

\begin{eqnarray}
C_q^{NS(0,0)} &=& 1 \\
C_q^{NS(1,1)} &=& -C_q^{NS(0, 0)} \gamma_{qq}^{NS(0)} \\
C_q^{NS(2,2)} &=& \frac{1}{2} \big(-\beta_0 C_q^{NS(1, 1)} - c_q^{NS(1, 1)} \gamma_{qq}^{NS(0)} \big) \\
C_q^{NS(2,1)} &=& -\beta_0 C_q^{NS(1, 0)} - C_q^{NS(1, 0)} \gamma_{qq}^{NS(0)} - C_q^{NS(0, 0)} \gamma_{qq}^{NS(1)}.
\end{eqnarray}

Formulas for $C_g^{S(k,\ell)}$ and $C_q^{S(k,\ell)}$ can be derived in a very similar way.

In the case of the structure function $g_1$, the massless Wilson coefficients are not currently known to 3-loop accuracy. In the formulas reported in Appendix <ref> and Appendix <ref>, the terms $C_i^{(3)}$ are all left symbolic. However, they cannot affect the logarithmic terms $\ln(Q^2/m^2)$.

The form of the renormalized massive OMEs has been derived in [158, 175]. In [199] we have collected the explicit form of the renormalized OMEs in the Larin scheme, in terms of harmonic sums [246, 247, 255, 256]. These formulas, which we do not repeat here due to their length, have been collected using the Mathematica packages HarmonicSums and Sigma.

The explicit results on the logarithmic corrections are reported in Appendix <ref> and Appendix <ref>.

§ N$^3$LO SCHEME-INVARIANT EVOLUTION OF THE NON-SINGLET STRUCTURE FUNCTIONS $F_2^{NS}$ AND $G_1^{NS}$

The measurement of DIS structure functions provides an ideal framework for the measurement of the strong coupling constant $a_s(M_Z)$, which is known today with a precision of $\mathcal O(1\%)$. Fits of the experimental data are typically performed by parametrizing the non-perturbative parton distribution functions through an appropriate functional form with free parameters, and by building the theoretical prediction for the structure functions through a convolution with the DIS Wilson coefficients. A global fit then delivers the parameters and the measurement of $a_s(M_Z)$ by an error minimization procedure, along with the respective uncertainties, see e.g. [211] for a study of $F_2$ and [212] for a study of $g_1$.

In the non-singlet case, it is also possible to directly compare the structure functions at two virtualities, considering an experimentally determined input at a starting scale $Q_0^2$ and fitting data at $Q^2>Q_0^2$. In the massless case, the relationship between the structure function at two different virtualities only depends on one parameter, namely $a_s$, and is scheme-invariant, i.e. independent on the factorization scheme which defines the factorization of the structure function into the Wilson coefficients and the parton distribution functions [41, 223], see also [292] for the explicit treatment of the singlet case. In principle, this makes it possible to perform a one-parameter fit for $a_s$, potentially reducing the uncertainty on its measurement.

This type of framework could be extended by considering massive quarks, by importing precision measurements of the masses $m_b$, $m_c$ from other sources. In [200], we completed the formalism of [223] by also considering the asymptotic heavy Wilson coefficients, and perform a numerical study of their effects. An implementation of the formalism has been developed in a numerical code which performs the evolution in $N$-space and obtains the structure functions in momentum fraction space by performing the inverse Mellin transform through a numerical integration in the complex $N$-plane. The evolution of the PDFs in $N$-space has been considered before in [293]; the analyses of [211, 212] have also been performed in $N$-space.

In an analysis on experimental data, it would be important to apply the cuts $W^2>15\text{~GeV}^2$, $Q^2>10\text{~GeV}^2$, with $W$ the invariant mass of the hadrons in the final state, in order to exclude higher-twist effects from the experimental sample and ensure that the description given by the asymptotic Wilson coefficients is sufficiently accurate.

§.§ Flavour decomposition

The non-singlet flavour distributions can be written as [293]

\begin{equation}
v_{k^2-1}^\pm = \sum_{l=1}^k (q_l\pm \bar q_l) - k(q_k\pm \bar q_k) ,
\end{equation}

with $q_i$ the quark distributions. For three light flavours, they are

\begin{eqnarray}
v_0^\pm &=& 0 \\
v_3^\pm &=& (u\pm\bar u) - (d\pm\bar d) \\
v_8^\pm &=& (u\pm\bar u) + (d\pm\bar d) - 2(s\pm\bar s)
\end{eqnarray}

and in general

\begin{eqnarray}
q_i+\bar q_i &=& \frac{1}{N_F} \Sigma - \frac{1}{i}v_{i^2-1}^+ + \sum_{l=i+1}^{N_F} \frac{1}{l(l-1)} v_{l^2-1}^+, \\
\Sigma &=& \sum_{l=1}^{N_F}(q_l+\bar q_l) ,
\end{eqnarray}

and the nucleon structure functions are given by

\begin{eqnarray}
F_2^p &=& x \Big[ \frac{2}{9} \Sigma + \frac{1}{6} v_3^+ + \frac{1}{18} v_8^+ \Big] \\
F_2^d &=& \frac{1}{2} [ F_2^p + P_2^n ] = x \Big[ \frac{2}{9} \Sigma + \frac{1}{18} v_8^+ \Big]
\end{eqnarray}

with a similar relation for $g_1^{p,d}$. A projection on the singlet distribution would require charged current structure functions [134]

\begin{equation}
\frac{1}{2}[W_2^{p,+} + W_2^{p,-}] = x \Sigma
\end{equation}

with the index $\pm$ indicating the exchange of a $W^+$ or $W^-$ boson. We will consider the flavour non-singlet combinations

\begin{eqnarray}
F_2^{NS}   &=& F_2^p - F_2^d = \frac{1}{6} x C_q^{NS,+} \otimes v_3^+ , \\
x g_1^{NS} &=& \frac{1}{6} x \Delta C_q^{NS,+} \otimes \Delta v_3^+ .
\end{eqnarray}

§.§ The non-singlet evolution

We derive a relation, valid in $N$-space,

\begin{equation}
F^{NS}(Q^2) = E_{NS}(Q^2,Q_0^2) F_{NS}(Q_0^2),
\end{equation}

where $F^{NS}$ refers to $F_2^{NS}$ or to $g_1^{NS}$, solving the evolution equation

\begin{equation}
\label{eq:nsEvolution}
\frac{d}{dt}\ln\left[F^{\NS}(Q^{2})\right]=\frac{d}{dt}\ln\left[C^{\NS}(Q^{2})\right]+\frac{d}{dt}\ln\left[q^{\NS}(Q^{2})\right] .
\end{equation}

The Wilson coefficient is given by

\begin{equation}
\label{eq:WilCoeff}
C(Q^{2})=1+\sum_{k=1}^{\infty}a_{s}^{k}(Q^{2}) C_k, \quad C_k = c_k + h_k(L_c,L_b).
\end{equation}

Here $c_k$ denote the expansion coefficients of the massless Wilson coefficients and $h_k$ of the massive Wilson coefficient, with

\begin{eqnarray}
L_c=\ln\frac{Q^2}{m_c^2}, \quad L_b=\ln\frac{Q^2}{m_b^2} .
\end{eqnarray}

In the non-singlet case the heavy flavor corrections contribute from $\mathcal{O}(a_s^2)$ onward. One has

\begin{eqnarray}
\\
\end{eqnarray}

where $\hat{h}$ denote the single mass and $\hat{\hat{h}}$ the double mass contributions.

One may rewrite the differential operator

\begin{equation}
\frac{d}{d \ln(Q^2)} = \frac{d a_s(Q^2)}{d\ln(Q^2)} \frac{d}{d a_s(Q^2)}
\end{equation}

\begin{equation}
\frac{d a_s}{d\ln(Q^2)} = -\sum_{k=0}^\infty \beta_k a_s^{k+2} .
\label{eq:AS}
\end{equation}

The evolution equation for the non-singlet quark density is

\begin{equation}
\label{eq:AP}
\frac{d}{dt}\ln\left[q^{\NS}(Q^{2})\right]= \sum_{k=0}^\infty P_{k}a^{k+1}(Q^{2}) ,
\end{equation}

where $\beta_k$ are the expansion coefficients of the QCD-$\beta$ function and $P_{k,qq}^{NS}\equiv P_k^{NS}$ are the splitting functions. The anomalous dimensions are related to the splitting functions by[Our normalizations are such that a factor of two has to be applied to those given in [69, 70]]

\begin{equation}
\gamma_{ij,(k)}(N)=-\int_{0}^{1}dxx^{N-1}P_{ij,(k)}(x) ,
\end{equation}

where $\gamma_{ij,(k)}(N)$ are the expansion coefficients of the non-singlet anomalous dimensions.

One obtains to $N^3LO$

\begin{eqnarray}
E_{NS}(Q^2,Q_0^2) & =&\left(\frac{a}{ a_0 }\right)^{-\frac{ P_0 }{2\beta_0 }}
\Biggl\{
+\frac{a- a_0 }{2\beta_0 ^2} \biggl\{
\Bigl[1 + a^2  C_2(Q^2) - a_0 ^2  C_2(Q_0^2) \Bigr] \bigl(2 \beta_0 ^2  C_1 -\beta_0   P_1 +\beta_1  P_0 \bigr)
\nonumber\\&&
-\frac{\bigl(a^2- a_0 ^2\bigr)}{4 \beta_0 ^3} \bigl(2\beta_0 ^2  C_1 -\beta_0   P_1 +\beta_1   P_0 \bigr) \Bigl[2\beta_0 ^3 C_1 ^2+\beta_0 ^2 P_2 - \beta_0 \beta_1 P_1 + \bigl(\beta_1 ^2-\beta_0  \beta_2 \bigr) P_0 \Bigr]
\nonumber\\&&
+\frac{\bigl(a^2+a  a_0 + a_0 ^2\bigr)}{3 \beta_0 ^2} \Bigl[2\beta_0 ^4  C_1 ^3-\beta_0 ^3  P_3 +\beta_0 ^2 \beta_1 P_2 + \bigl(\beta_0 ^2 \beta_2 -\beta_0  \beta_1 ^2\bigr) P_1
\nonumber\\&&
+ \bigl(\beta_0 ^2 \beta_3 -2 \beta_0  \beta_1  \beta_2 +\beta_1 ^3\bigr) P_0 \Bigr]
+\frac{a- a_0 }{4 \beta_0 ^2} \bigl(2\beta_0 ^2  C_1 -\beta_0   P_1 +\beta_1   P_0 \bigr)^2
\nonumber\\&&
+\frac{(a- a_0 )^2}{24 \beta_0 ^4} \bigl(2\beta_0 ^2  C_1 -\beta_0   P_1 +\beta_1   P_0 \bigr)^3
-\frac{a+ a_0 }{2 \beta_0 } \Bigl[2\beta_0 ^3  C_1 ^2+\beta_0 ^2  P_2 -\beta_0  \beta_1 P_1
\nonumber\\&&
+ P_0  \bigl(\beta_1 ^2-\beta_0  \beta_2 \bigr)
\Bigr]
\biggr\}
+a^2  C_2(Q^2) - a_0 ^2  C_2(Q_0^2)
- C_1  \Bigl[a^3  C_2(Q^2) - a_0 ^3  C_2(Q_0^2) \Bigr]
\nonumber\\&&
+ a^3 C_3(Q^2) - a_0 ^3  C_3(Q_0^2)
\Biggr\} .
\label{eq:ENS}
\end{eqnarray}

Here we used the notation $a=a_s(Q^2)$, $a_0=a_s(Q_0^2)$ and $P_i=P_{qq,(i)}^{NS}$, considered in $N$-space.

The heavy quark contributions to the Wilson coefficients are given by [181, 73, 184, 158]

\begin{eqnarray}
\label{H2}
\hat{h}_2^{(Q)} &=& - \frac{\beta_{0,Q}}{4} P_{qq,(0)} \ln^2\left(\frac{Q^2}{m^2}\right)
+ \frac{1}{2} \hat{P}_{qq,(1)}^{\rm NS} \ln\left(\frac{Q^2}{m^2}\right)
+ a_{qq}^{(2), \rm NS} + \frac{\beta_{0,Q}}{4} \zeta_2 P_{qq,(0)} + \hat{C_q}^{(2),\rm NS}
\\
\label{H3}
\hat{h}_3^{(Q)} &=& -\frac{1}{6} P_{qq,(0)} \beta_{0,Q} \left(\beta_0 + 2 \beta_{0,Q}\right)
\ln^3\left(\frac{Q^2}{m^2}\right)
+ \frac{1}{4} \Biggl[-2 P_{qq,(1)}^{\rm NS} \beta_{0,Q} + 2 \hat{P}_{qq,(1)}^{\rm NS} \left(\beta_0 +
\beta_{0,Q}\right)
\nonumber\\ &&
- \beta_{1,Q} P_{qq,(0)} \Biggr] \ln^2\left(\frac{Q^2}{m^2}\right)
- \frac{1}{2} \Biggl[- \hat{P}_{qq,(2)}^{\rm NS} - \left(4 a_{qq,Q}^{(2),\rm NS} + \zeta_2 \beta_{0,Q}
P_{qq,(0)}\right)  \left(\beta_0 +  \beta_{0,Q}\right)
\nonumber\\ &&
- P_{qq,(0)} \beta_{1,Q}^{(1)}\Biggr]
\ln\left(\frac{Q^2}{m^2}\right)
+ 4 \bar{a}_{qq,Q}^{(2),\rm NS}(\beta_0+\beta_{0,Q}) + P_{qq,(0)} \beta_{1,Q}^{(2)}
+ \frac{1}{6} P_{qq,(0)} \beta_0 \beta_{0,Q} \zeta_3
\nonumber\\ &&
+ \frac{1}{4} P_{qq,(1)}^{\rm NS} \beta_{0,Q} \zeta_2
- 2 \delta m_1^{(1)} \beta_{0,Q} P_{qq,(0)} - \delta m_1^{(0)} \hat{P}_{qq,(1)}^{\rm NS} + 2 \delta m_1^{(-1)}
a_{qq,Q}^{(2),\rm NS} + a_{qq,Q}^{(3),\rm NS}
\nonumber\\ && + \Biggl[- \frac{\beta_{0,Q}}{4} P_{qq,(0)} \ln^2\left(\frac{Q^2}{m^2}\right)
+ \frac{1}{2} \hat{P}_{qq,(1)}^{\rm NS} \ln\left(\frac{Q^2}{m^2}\right)
+ a_{qq}^{(2), \rm NS} + \frac{\beta_{0,Q}}{4} \zeta_2 P_{qq,(0)}\Biggr] C_q^{(1),\rm NS}
\nonumber\\ &&
+ \hat{C}_q^{(3),\rm NS}.
\end{eqnarray}

and the two-mass contribution by [121]

\begin{eqnarray}
\label{H23}
\hat{\hat{h}}_3^{\rm NS} &=& P_{qq}^{(0)} \beta_{0,Q}^2\left[\frac{2}{3} \left(L_c^3 + L_b^3\right) + \frac{1}{2}
\left(L_c^2 L_b +L_c L_b^2 \right)\right] - \beta_{0,Q} \hat{P}_{qq}^{(1),\rm NS} \left(L_c^2 + L_b^2\right)
-\left[4 a_{qq,Q}^{(2),\rm NS} \beta_{0,Q} \right.
\nonumber\\ && \left.
- \frac{1}{2} \beta_{0,Q}^2 P_{qq}^{(0)} \zeta_2\right] (L_c + L_b)
+ 8 \bar{a}_{qq,Q}^{(2),\rm NS} \beta_{0,Q} + \tilde{a}^{(3),\rm NS}_{qq,Q}(m_c,m_b,Q^2).
\end{eqnarray}

The two-mass term is the same in the unpolarized and polarized case. We employed the definition

\begin{eqnarray}
\hat{f}(x,N_F) = f(x,N_F+1) - f(x,N_F).
\end{eqnarray}

The perturbative solution for $a_s(Q^2)$ is given in the $\overline{\sf MS}$ scheme by [294]

\begin{eqnarray}
a_s(Q^2) &=& \frac{1}{\beta_0 L}
- \frac{\beta_1}{\beta_0^3 L^2} \ln(L)
+ \frac{1}{\beta_0^3 L^3} \Biggl[\frac{\beta_1^2}{\beta_0^2 }(\ln^2(L) - \ln(L) -1) + \frac{\beta_2}{\beta_0}\Biggr]
\nonumber\\ &&
+ \frac{1}{\beta_0^4 L^4} \Biggl[\frac{\beta_1^3}{\beta_0^3 }\Biggl(-\ln^3(L) + \frac{5}{2} \ln^2(L) + 2 \ln(L) -
\frac{1}{2}\Biggr) - 3 \frac{\beta_1 \beta_2}{\beta_0^2} \ln(L) + \frac{\beta_3}{2 \beta_0}\Biggr],
\label{eq:asSol}
\end{eqnarray}

with $L = \ln(Q^2/\Lambda_{\rm QCD}^2)$. Here the integration constant for solving (<ref>) is chosen by
$(\beta_1/\beta_0^2) \ln(\beta_0)$ [40]. The expansion coefficients of the
$\beta$-function to N$^3$LO were calculated in [36, 35, 93, 94, 95, 96, 97, 98].
The flavor matching conditions were given in [294].
The expansion coefficients of the renormalized mass were given in  [234, 235].
The constant and $O(\ep)$ parts of the massive unrenormalized operator matrix elements at $O(a_s^k)$
are denoted by $a_{ij}^{(k)}$ and $\bar{a}_{ij}^{(k)}$, respectively, cf. [159, 160, 161, 295].

In the numerical evaluation shown below, an approximate form is used for the three-loop massless Wilson coefficients. For all other Wilson coefficients, the analytic Mellin-space representations are used. After reduction to a basis of independent harmonic sums [255], the objects depend on 32 harmonic sums [246, 247] up to weight 6; with weight-6 sums appearing only in the 3-loop Wilson coefficient.

The harmonic sums in $N$ space are calculated in the complex plane after representing them using the Mellin transforms of harmonic polylogarithms, and computing their asymptotic expansions [247, 256, 257]. Together with exact step relations in $N$, an accurate numerical evaluation of the harmonic sums becomes possible in the complex plane. The analytic continuation of harmonic sums needed for the anomalous dimensions up to 3 loops has first been performed in [258].

In the case of $F_2^{NS}$, the relevant splitting functions are $P_{qq,(k)}^{NS,+}$, whereas for $g_1^{NS}$ they are $P_{qq,(k)}^{NS,-}$. The massless Wilson coefficients have been calculated in [40, 62, 63, 64, 71] and [88, 89] respectively. For the four-loop splitting functions, which are not currently known in analytic form, we employ in numerical illustrations below the Padé approximant

\begin{eqnarray}
\label{eq:PADE}
P_{qq}^{3, \pm, \rm NS}(N) \approx \frac{P_{qq}^{2, \pm, \rm NS}(N)^2}{P_{qq}^{1, \pm, \rm NS}(N)} .
\end{eqnarray}

Low moments of the four-loop splitting functions have been calculated in [296, 297, 298, 299, 300, 301]. A comparison of these exact moments to the Padé approximant is shown in Table <ref>.

A previous analysis [211] has showed that a $100\%$ error on $P_{qq,(3)}^{NS}$ would determine an error of 2 MeV on $\Lambda_{QCD}$, well below the experimental error currently of $\delta\Lambda_{QCD}=26\text{~MeV}$.

The leading small-$x$ terms of $P_3^{NS,+}$ and $P_2^{NS,-}$ have been studied in [302, 303, 304]; the leading large-$N_F$ behaviour of the splitting functions has also been given in [305].

1|c$\delta \gamma^{+,\rm NS}$
1|c|$\delta \gamma^{-,\rm NS}$

2     0.208822541      1      0.0

4     0.123728742      3      0.147102092

6     0.087155544      5      0.101634935

8     0.064949195      7      0.074593595

10     0.049680399      9      0.056598595

12     0.038394815     11      0.043633919

14     0.029638565     13      0.033767853

16     0.022602035     15      0.025956941

The relative error comparing the exact moments of the four–loop anomalous dimensions, $\gamma^{(3),\pm, \rm
NS}$, with the Padé approximation (<ref>).

§.§ Numerical results

In our numerical illustration we employ the values of the charm and bottom quarks of $m_c = 1.59\text{~GeV}$ [286] and $m_b = 4.78\text{~GeV}$ [306].

The input structure function in the unpolarized case is built from the non-singlet parton distribution [211]

\begin{eqnarray}
xq^{\rm NS}(x,Q_0^2) &=&
\frac{1}{3} \biggl[
0.262 \ x^{0.298} (1-x)^{4.032} (1+6.042 \sqrt{x} +35.49 x)
\nonumber\\&&
~~~~~-1.085 \ x^{0.5} (1-x)^{5.921} (1-3.618 \sqrt{x} +16.41 x) \biggr]
\end{eqnarray}

at $Q_0^2=4\text{~GeV}^2$; for the polarized case we use a fit of the structure function of [212] at $Q_0^2 = 10~\text{GeV}^2$.

We also employ for the purpose of illustration the quantity

\begin{eqnarray}
F_2^{\rm h}(N,Q^2) = \left[E_{\rm NS} - \left. E_{\rm NS}\right|_{\rm h=0}\right] F_2(N,Q_0^2).
\end{eqnarray}

In Fig. <ref> we show the scheme-invariant evolution of the non-singlet structure functions $F_2^{NS}$ and $xg_1^{NS}$, in the kinematic region $Q^2\in [10, 10^4]\text{~GeV}^2$. In Fig. <ref> we expand the representation for the region of larger values of $x$.
In Fig. <ref> we illustrate the relative effect of the scale evolution in $Q^2$ both for $F_2^{NS}$ and $xg_1^{NS}$ comparing to the starting scale $Q_0^2$.
In Fig. <ref> we show the ratio of the results obtained at leading order (LO), next-to-leading order (NLO), and next-to-next-to leading order (NNLO) to the N$^3$LO results at $Q^2 = 100~\GeV^2$.
In Fig. <ref> we illustrate the relative size of the heavy flavor parts for the same region in $Q^2$ in the unpolarized and polarized cases.  In the important region $x \leq 0.4$ the heavy flavor corrections reach the size of $\sim 1\%$.
In Fig. <ref> we illustrate the effect of the half difference if putting $P_{qq}^{3, \pm, \rm NS} = 2 {P_{qq}^{2, \pm, \rm NS}}^2/ P_{qq}^{2, \pm, \rm NS}$ and $P_{qq}^{3, \pm, \rm NS} = 0$ for both $F_2^{\rm NS}$ and $xg_1^{\rm NS}$. This rescaled correction is in the sub–percent range. Moreover, the impact on $\Lambda_{\rm QCD}$ comes from the slope in $Q^2$ which is seen to be rather small.

Left: The structure function $F_2^{\rm NS}$ at N$^3$LO. Right: The structure function
$xg_1^{\rm NS}$ at N$^3$LO. Full lines: $Q^2 = 10~\GeV^2$; dashed lines: $100~\GeV^2$;
dash-dotted lines: $1000~\GeV^2$;
dotted lines: $10000~\GeV^2$.

The same as Figure 1 but with expanded large $x$ region.

Left: The relative contribution of $F_2^{\rm NS}$ in the evolution from $Q^2 = 10~\GeV^2$ to $10000~\GeV^2$.
Right: The same for the structure function $xg_1^{\rm NS}$.

Left: The relative contributions from LO (dotted lines), NLO (dashed lines) and NNLO (full lines)
to the structure function
$F_2^{\rm NS}$ at N$^3$LO at $Q^2 = 100~\GeV^2$ as an example. Right: The same for the structure function $xg_1^{\rm NS}$.

Left: The relative contribution of the heavy flavor contributions due
to $c$ and $b$ quarks to the structure function $F_2^{\rm NS}$ at N$^3$LO;
dashed lines: $100~\GeV^2$;
dashed-dotted lines: $1000~\GeV^2$;
dotted lines: $10000~\GeV^2$. Right: The same
for the structure function $xg_1^{\rm NS}$ at N$^3$LO.

The effect of the variation of $P_{qq}^{3, \pm, \rm NS}$ around the value in Eq. (<ref>) by $\pm
100\%$. Dashed lines: $Q^2 = 100~\GeV^2$, Dash-dotted lines: $Q^2 = 1000~\GeV^2$; Dotted lines: $Q^2 = 10000~\GeV^2$.
Left: $F_2^{\rm NS}$; Right: $xg_1^{\rm NS}$.

§ HYPERGEOMETRIC FUNCTIONS AND DIFFERENTIAL SYSTEMS

Hypergeometric functions are known to occur in the computation of Feynman integrals [229, 253, 160, 308, 307]. In the mathematical literature, hypergeometric functions of various kinds have been studied by a number of authors [309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325].

A major methodological advance in the calculation of Feynman diagrams was the development of the method of integration by parts [326, 327, 328, 329, 330, 331, 332, 333, 334], used together with systems of differential equations obeyed by sets of Feynman integrals [335, 336, 337, 338, 254, 339, 264, 340]. Systems obtained as the result of IBP reduction can be decoupled [342, 341] using methods available in the package OreSys [343]. In the case of partial differential equations, there exists the method of Janet bases [344].

Calculations in which hypergeometric functions have played a role are among others [345, 346, 347, 348, 349] at one loop, multileg level, and [350, 351, 352, 353] at higher loop order.

Other authors [355, 354] have emphasized the relationship between Feynman integrals and Gelfand-Kapranov-Zelevinsky hypergeometric systems [356, 357, 358, 359], which have solutions in terms of hypergeometric functions.

In [201] we systematized the study of a class of differential systems obeyed by certain classes of hypergeometric functions and investigated in some concrete examples how the series expansion of hypergeometric functions with respect to a parameter is obtained in the context of difference ring methods [360, 361, 362, 363, 364, 365, 366, 368, 367] implemented in the package Sigma [202, 203, 204] and EvaluateMultiSums [282].

We follow here the exposition of [201], in which we present our package HypSeries which employs strategies useful for the classifications of systems of differential equations having hypergeometric solutions.

§.§ Differential equations

Hypergeometric functions satisfy differential equations and have integral representations. For instance, the function $_2F_1$ has the property

\begin{eqnarray}
\pFq{2}{1}{a_1,a_2}{b_1}{z} = \frac{\Gamma(b_1)}{\Gamma(a_1) \Gamma(b_1-a_1)} \int_0^1 dx x^{a_1-1} (1-x)^{b_1-a_1-1}
\end{eqnarray}

and satisfies the differential equation

\begin{eqnarray}
x(1-x) \frac{d^2}{dx^2} + (c -(a+b+1)x) \frac{d}{dx} - ab ~.
\label{eq:D2F1g}
\end{eqnarray}

For the function $_{3}F_2$ one obtains

\begin{eqnarray}
x^2(1-x) \frac{d^3}{dx^3} + x (A_{{2}}+B_{{2}} x) \frac{d^2}{dx^2}
+ (A_{{1}}+B_{{1}} x) \frac{d}{dx} + C,
\label{eq:D3F2}
\end{eqnarray}

with $A_{{2}} = b_1+b_2 +1, B_{{2}} = -(3 + a_1 + a_2 +a_3), A_{{1}} = b_1b_2,
B_{{1}} = -(a_2 a_1 + a_3 a_1 + a_2 a_3 + a_1 +a_2 +a_3 +1), C = - a_1 a_2 a_3$,
while $_{p+1}F_p$ can be represented by the iterative integral

\begin{eqnarray}
\int_0^1 dx x^{a-1} (1-x)^{b-1}
\pFq{p}{q}{a_1,...,a_p}{b_1,...b_q}{x z} = \frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)} \pFq{p+1}{q+1}
\end{eqnarray}

and satisfies the differential equation

\begin{eqnarray}
x^p (1-x) \frac{d^{p+1}}{dx^{p+1}} + \sum_{k = 1}^{p} x^{k-1}(A_k + B_k x) \frac{d^k}{dx^k} + C
\end{eqnarray}

since it is annihilated by the differential operator

\begin{eqnarray}
x \frac{d}{dx} \Biggl( x \frac{d}{dx} + b_1 -1 \Biggr) ... \Biggl( x \frac{d}{dx} + b_q -1 \Biggr)
- x \Biggl( x \frac{d}{dx} + a_1  \Biggr) ... \Biggl( x \frac{d}{dx} + a_p \Biggr).
\label{eq:PFQ1}
\end{eqnarray}

The products of the differential operators in (<ref>) $\vartheta = x (d/dx) \equiv x \partial_x$,
can be written in the following form

\begin{eqnarray}
\vartheta &=& x \partial_x
\\
\vartheta^2 &=& x \partial_x + x^2 \partial_x^2
\\
\vartheta^3 &=& x \partial_x + 3 x^{{2}} \partial_x^2 + x^3 \partial_x^3
\\
\vartheta^4 &=& x \partial_x + 7 x^{{2}} \partial_x^2 + 6 x^3 \partial_x^3+ x^4 \partial_x^4
\\
\vartheta^5 &=& x \partial_x + 15 x^{{2}} \partial_x^2 + 25 x^3 \partial_x^3+ 10 x^4 \partial_x^4 + x^5
\partial_x^5,~~\text{etc.}
\end{eqnarray}

One can parametrize the equations obeyed by the Horn hypergeometric functions [312, 313, 314, 315, 317] $F_1$ to $F_4$, $G_1$ to $G_3$, and $H_1$ to $H_7$, including the Appell functions [312, 313] as follows:

\begin{eqnarray}
+ (b x+c) \partial_x
+x (d+e x) \partial_x^2
+f y \partial_y
+(g y+h x y) \partial_{x,y}^2
+j y^2 \partial_y^2 &=& 0
\label{APP1}
\\
+ (b_1 y+ c_1) \partial_y
+y (d_1+e_1
y) \partial_y^2
+f_1 x \partial_x
+ (g_1  x+ h_1 x y) \partial_{x,y}^2
+j_1 x^2 \partial_x^2 &=& 0,
\label{APP2}
\end{eqnarray}

with the example of the Appell $F_1$ function

\begin{eqnarray}
F_1 &:& x(1-x)\partial_x^2 +y(1-x) \partial_{x,y}^2 + (A+B x)\partial_x + C y\partial_y + D
\\
F_1 &:& y(1-y)\partial_y^2 +x(1-y) \partial_{x,y}^2 + (A +B' y)\partial_y + C' x\partial_x + D'.
\end{eqnarray}

In physical applications, the functions $S_1$ and $S_2$ appeared in [350, 351]. They are annihilated by the differential operators

\begin{eqnarray}
S_1 &:&
a+(c+b x) \partial_x +x (d+e x) \partial_x^2 +x^2 (l+p x) \partial_x^3+f y \partial_y +x
(q+r x) y \partial_x^2 \partial_y +j y^2 \partial_y^2
\nonumber\\ &&
+s x y^2 \partial_x \partial_y^2 +(g y+h x y)
\partial_{x,y}^2
\\
a_1+f_1 x \partial_x + j_1 x^2 \partial_x^2 +(c_1+b_1 y) \partial_y+y (d_1+e_1 y) \partial_y^2 +(g_1 x+h_1 x
y) \partial_{x,y}^2
\\
S_2 &:&
a+(c+b x) \partial_x +x (d+e x) \partial_x^2 +x^2 (l+p x) \partial_x^3 +f y \partial_y +x
(q+r x) y \partial_x^2 \partial_y
+j y^2 \partial_y^2
\nonumber\\ &&
+s x y^2 \partial_x \partial_y^2 +(g y+h x y)
\partial_{x,y}^2
\\
a_1+f_1 x \partial_x +(c_1+b_1 y) \partial_y +j_1 x^2 \partial_x^2 \partial_y +y (d_1+e_1 y)
\partial_y^2+q_1 x y \partial_x \partial_y^2 +p_1 y^2 \partial_y^3
\nonumber\\ &&
+(g_1 x+h_1 x y) \partial_{x,y}^2.
\end{eqnarray}

For the Kampé de Fériet function

\begin{eqnarray}
F^{p;q;k}_{l;m;n} \left[ \begin{array}{cc}
(a_p) ; (b_q) ; (c_k) & \\
& x, y \\
(\alpha_l) ; (\beta_m) ; (\gamma_n) &
\end{array}
\right]
&=& \sum_{r,s=0}^\infty \, \frac{\prod_{j=1}^p (a_j)_{r+s}
\prod_{j=1}^q (b_j)_r \prod_{j=1}^k (c_{{j}})_s}{\prod_{j=1}^l (\alpha_j)_{r+s} \prod_{j=1}^m (\beta_j)_r \prod_{j=1}^n
(\gamma_j)_s } \frac{x^r}{r!} \frac{y^s}{s!} \\
&=& \sum_{r,s=0}^\infty f[r,s] x^r y^s
\end{eqnarray}

one obtains the following annihilating differential operators [314, 315]

\begin{eqnarray}
\prod_{j=1}^p (x \partial_x + y \partial_y +a_j) \prod_{j=1}^q (x \partial_x +b_j) -{\partial_x} \prod_{j=1}^l (x \partial_x + y \partial_y -1+\alpha_j) \prod_{j=1}^m ({x \partial_x} -1+\beta_j) = 0,
\\
\prod_{j=1}^p (x \partial_x +y \partial_y +a_j) \prod_{j=1}^k (y \partial_y +c_j) -{\partial_y} \prod_{j=1}^l (x
\partial_x + y \partial_y -1+\alpha_j) \prod_{j=1}^n (y \partial_y -1+\gamma_j) = 0.
\end{eqnarray}

The differential operators for the triple hypergeometric series [321] read

\begin{eqnarray}
D_{3,1} &=& A+(B_0+B_1 x) \partial_x +x (E_0+E_1 x) \partial_x^2 + C_1 y \partial_y +{F_1} y^2 \partial_y^2 +(H_0+H_1 x) y
\partial_{x,y}^2
\nonumber\\ &&
+D_1 z \partial_z
+G_1 z^2 \partial_z^2 +(L_0+L_1 x) z \partial_{x,z}^2 + S_1 y z \partial_{y,z}^2
\\
D_{3,2} &=& A'+{B_1'} x \partial_x +E_1' x^2 \partial_x^2 +(C_0'+C_1' y) \partial_y +y (F_0'+F_1' y)
\partial_y^2+x
(H_2'+H_1' y) \partial_{x,y}^2
\nonumber\\ &&
+ D_1' z \partial_z
+ G_1' z^2 \partial_z^2 + L_1' x z \partial_{x,z}^2
+(S_0'+S_1' y) z \partial_{y,z}^2
\\
D_{3,3} &=&
A''+ B_1'' x \partial_x + {E''_1} x^2 \partial_x^2 + C_1'' y \partial_y + {F_1''} y^2 \partial_y^2
+ H_1'' x y \partial_{x,y}^2 + (D''_0 + D_1'' z) \partial_z
\nonumber\\ &&
+ z (G_0'' + G_1'' z) \partial_z^2
+x {(L_2'' + L''_1 z)} \partial_{x,z}^2 + y (S_2''+ S_1'' z) \partial_{y,z}^2 .
\end{eqnarray}

The differential operators for the quadruple versions are given by

\begin{eqnarray}
D_{4,1} &=&
A + E_1 t \partial_t +L_1 t^2 \partial_t^2 +(B_0+B_1 x)\partial_x +x (F_0+F_1 x) \partial_x^2+t
(P_0+P_1 x) \partial_{t,x}^2 +C_1 y \partial_y
\nonumber\\ &&
+G_1 y^2 \partial_y^2 +R_1 t y \partial_{t,y}^2 + (M_0+M_1 x) y
\partial_{x,y}^2 + D_1 z \partial_z + H_1 z^2 \partial_z^2 + S_1 t z \partial_{t,z}^2
\nonumber\\ &&
+(N_0+N_1 x) z
\partial_{x,z}^2 +Q_1 y z \partial_{y,z}^2
\\
D_{4,2} &=&
A'+E_1' t \partial_t + L_1' t^2 \partial_t^2 + B_1' x \partial_x +F_1' x^2 \partial_x^2 + P_1' t x
\partial_{t,x}^2+(C_0'+ C_1' y) \partial_y
+y (G_0' + G_1' y) \partial_y^2
\nonumber\\ &&
+ t (R_0'+R_1' y) \partial_{t,y}^2 + x (M_2' + M_1' y)
\partial_{x,y}^2 + D_1' z \partial_z +H_1' z^2 \partial_z^2 +S_1' t z \partial_{t,z}^2 + N_1' x z
\partial_{x,z}^2
\nonumber\\ &&
+ (Q_0'+Q_1' y) z \partial_{y,z}^2
\\
D_{4,3} &=&
A''+E_1'' t \partial_t + L_1'' t^2 \partial_t^2 + B_1'' x \partial_x + F_1'' x^2 \partial_x^2
+ P_1'' t x \partial_{t,x}^2  + C_1'' y \partial_y + G_1'' y^2 \partial_y^2 + R_1'' t y \partial_{t,y}^2
\nonumber\\ &&
+ M_1'' x
y \partial_{x,y}^2
+(D_0'' + D_1'' z) \partial_z + z (H_0''+H_1'' z) \partial_z^2 + t (S_0''+S_1'' z) \partial_{t,z}^2+x
(N''_2 + N_1'' z) \partial_{x,z}^2
\nonumber\\ &&
+ y (Q_2''+Q_1'' z) \partial_{y,z}^2
\\
D_{4,4} &=&
A''' + (E_0''' + E_1''' t) \partial_t + t (L_0''' + L_1''' t) \partial_t^2 + B_1''' x \partial_x
+ F_1''' x^2 \partial_x^2 + (P_2''' + P_1''' t) x \partial_{t,x}^2
\nonumber\\ &&
+ C_1''' y \partial_y + G_1''' y^2
\partial_y^2 + (R_2''' + R_1''' t) y \partial_{t,y}^2 + M_1''' x y \partial_{x,y}^2 + D_1''' z \partial_z + H_1''' z^2
\partial_z^2
\nonumber\\ &&
+ (S_2'''+ S_1''' t) z \partial_{t,z}^2 + N_1''' x z \partial_{x,z}^2 + Q_1''' y z \partial_{y,z}^2.
\end{eqnarray}

They cover the functions $K_i,~~i = 1 ... 21$ of Refs. [318, 319].

§.§ Recursions

The formal power series ansatz

\begin{eqnarray}
\sum_{k_1, ..., k_n = 0}^\infty f[k_1, ..., k_n] x_1^{k_1} ... x_n^{k_n}
\label{eq:TAYLOR}
\end{eqnarray}

with $f[k_1, ..., k_n]$ hypergeometric induces a recurrence relation on the expansion coefficients, deriving from the differential equations of Section <ref>. (In addition, hypergeometric functions obey contiguous relations in their parameters [369, 307, 370]).

We list here the recursions corresponding to the differential systems from the previous Section:

\begin{eqnarray}
_2F_1 &:&
(C + n (1 - n + B_1)) f[n] + (1 + n) (n + A_1) f[1 + n] = 0
\label{eq:R2F1}
\\
_3F_2 &:&
\big( n B_1 + (n-1)n B_2 +C-(n-2) (n-1) n \big) f[n]
\nonumber\\ &&
+ \big( (n+1) A_1 + n (n+1) A_2 +  (n-1)n(n+1)\big) f[n+1] = 0
\label{eq:R3F2}
\\
_{p+1}F_p &:&
\Biggl[\frac{C}{n!} - \frac{1}{(n-p-1)!} + \sum_{k=1}^p \frac{B_k}{(n-k)!} \Biggr] f[n]
\nonumber\\ &&
+ (n+1) \Biggl[ \frac{1}{(n-p)!} +  \sum_{k=1}^p \frac{A_k}{(n-k+1)!} \Biggr] f[n+1] = 0.
\end{eqnarray}

In the two–variable cases the expansion coefficients of the  Horn–type functions obey

\begin{eqnarray}
&&\Big[a+b m+e (m-1) m+n \big(f+h m+j (n-1)\big)\Big] f[m,n]
\nonumber\\&&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
+(1+m) (c+d m+g n) f[1+m,n]=0,
\\
&&\Big[a_1+f_1 m+j_1 (m-1) m+n \big(b_1+h_1 m+e_1 (n-1)\big)\Big] f[m,n]
\nonumber\\&&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
+(1+n) (c_1+g_1 m+d_1 n) f[m,1+n]=0,
\end{eqnarray}

and for the $S_1$-functions one has

\begin{eqnarray}
&&\Big[a+b m+e (m-1) m+f n+h m n+j (n-1) n+(m-2) (m-1) m p+(m-1) m n r
\nonumber\\&&
+m (n-1) n s\Big] f[m,n]
+(1+m) \Big[c+g n+m \big(d+l (m-1)+n q\big)\Big] f[1+m,n]=0,
\\
&&\Big[a_1+f_1 m+j_1 (m-1) m+n \big(b_1+h_1 m+e_1 (n-1)\big)\Big] f[m,n]
\nonumber\\&&
+(1+n) (c_1+g_1 m+d_1 n) f[m,1+n]=0,
\end{eqnarray}

as, likewise, for the

\begin{eqnarray}
&&\Big[a+b m+e (m-1) m+f n+h m n+j (n-1) n+(m-2) (m-1) m p+(m-1) m n r+
\nonumber\\&&
m (n-1) n s\Big] f[m,n]+(1+m) \Big[c+g n+m \big(d+l (-1+m)+n q\big)\Big] f[1+m,n]=0,
\\
&&\Big[a_1+f_1 m+n \big(b_1+h_1 m+e_1 {(n-1)}\big)\Big] f[m,n]+(1+n) \Big[c_1+n (d_1+(-1+n) p_1)
\nonumber\\&&
+m \big(g_1+j_1 (m-1)+n q_1\big)\Big] f[m,1+n]=0.
\end{eqnarray}

For the expansion coefficients $f[r,s]$ of the Kampé de Fériet functions the recurrences read

\begin{eqnarray}
\prod_{j=1}^p (r+s+a_j) \prod_{j=1}^q (r+b_j) f[r,s] - (r+1) \prod_{j=1}^l (r+s+\alpha_j) \prod_{j=1}^m (r+\beta_j)
f[r+1,s] = 0,
\\
\prod_{j=1}^p (r+s+a_j) \prod_{j=1}^k (s+c_j) f[r,s] - (s+1) \prod_{j=1}^l (r+s+\alpha_j) \prod_{j=1}^n (s+\gamma_j)
f[r,s+1] =0.
\end{eqnarray}

The coefficients in the 3-variable cases obey

\begin{eqnarray}
&&\Big[A+ B_1  m+ E_1  (m-1) m+ C_1  n+ H_1  m n+ F_1  (n-1) n+ D_1  p+ L_1  m p+ G_1  (p-1) p
\nonumber\\&&
+n p  S_1 \Big] f[m,n,p]+(1+m) ( B_0 + E_0  m+ H_0  n+ L_0  p) f[1+m,n,p]=0,
\\
&&\Big[ A' + B'_1  m+ E'_1  (m-1) m+ C'_1  n+ H'_1  m n+ F'_1  (n-1) n+ D'_1  p+ L'_1  m p+ G'_1  (p-1) p
\nonumber\\&&
+n p  S'_1 \Big] f[m,n,p]+(1+n) ( C'_0 + H'_2  m+ F'_0  n+p  S'_0 ) f[m,1+n,p]=0,
\\
&&\Big[ A'' + B''_1  m+ E''_1  (m-1) m+ C''_1  n+ H''_1  m n+ F''_1  (n-1) n+ D''_1  p+ L''_1  m p+ G''_1  (p-1) p
\nonumber\\&&
+n p  S''_1 \Big] f[m,n,p]+(1+p) ( D''_0 + L''_2  m+ G''_0  p+n  S''_2 ) f[m,n,1+p]=0.
\end{eqnarray}

For the 4-variable systems one has

\begin{eqnarray}
&&\Big[A+ B_1  m+ F_1  (m-1) m+ C_1  n+m  M_1  n+ G_1  (n-1) n+ D_1  p+m  N_1  p+ H_1  (p-1) p
\nonumber\\&&
+ E_1  q+m  P_1  q+ L_1  (q-1) q+n p  Q_1 +n q  R_1 +p q  S_1 \Big] f[m,n,p,q]+(1+m) ( B_0 + F_0  m
\nonumber\\&&
+ M_0  n+ N_0  p+ P_0  q) f[1+m,n,p,q]=0,
\\[2mm]
&&\Big[ A' + B'_1  m+ F'_1  (m-1) m+ C'_1  n+m  M'_1  n+ G'_1  (n-1) n+ D'_1  p+m  N'_1  p+ H'_1  (p-1) p
\nonumber\\&&
+ E'_1  q+m  P'_1  q+ L'_1  (q-1) q+n p  Q'_1 +n q  R'_1 +p q  S'_1 \Big] f[m,n,p,q]+(1+n) ( C'_0 +m  M'_2
\nonumber\\&&
+ G'_0  n+p  Q'_0 +q  R'_0 ) f[m,1+n,p,q]=0,
\\[2mm]
&&\Big[ A'' + B''_1  m+ F''_1  (m-1) m+ C''_1  n+m  M''_1  n+ G''_1  (n-1) n+ D''_1  p+m  N''_1  p+ H''_1  (p-1) p
\nonumber\\&&
+ E''_1  q+m  P''_1  q+ L''_1  (q-1) q+n p  Q''_1 +n q  R''_1 +p q  S''_1 \Big] f[m,n,p,q]+(1+p) ( D''_0 +m  N''_2
\nonumber\\&&
+ H''_0  p+n  Q''_2 +q  S''_0 ) f[m,n,1+p,q]=0,
\\[2mm]
&&\Big[ A''' + B'''_1  m+ F'''_1  (m-1) m+ C'''_1  n+m  M'''_1  n+ G'''_1  (n-1) n+ D'''_1  p+m  N'''_1  p
\nonumber\\&&
+ H'''_1  (p-1) p+ E'''_1  q+m  P'''_1  q+ L'''_1  (q-1) q+n p  Q'''_1 +n q  R'''_1 +p q  S'''_1 \Big] f[m,n,p,q]
\nonumber\\&&
+(1+q) ( E'''_0 +m  P'''_2 + L'''_0  q+n  R'''_2 +p  S'''_2 ) f[m,n,p,1+q]=0.
\end{eqnarray}

§.§ The Solution of the Recursions

Given a field $\mathbb K$ with characteristic zero, a multiple hypergeometric series is a function

\begin{equation}
f(x_1,\ldots,x_r) = \sum_{n_i\geq0} A(n_1,\ldots,n_r) x_1^{n_1}\cdots x_r^{n_r}
\label{eq:hyp-series}
\end{equation}

where $A:\mathbb N^r\to \mathbb K$, is hypergeometric, i.e.

\begin{equation}
s_i\,A(n_1,\ldots,n_i,\ldots,n_r) = t_i\,A(n_1,\ldots,n_i+1,\ldots,n_r) , \qquad i=1,\ldots,r
\label{eq:hyp-ratio}
\end{equation}

for polynomials $s_i,t_i\in \mathbb K[n_1,\ldots,n_r]$ being coprime.

It was remarked in the previous sections how hypergeometric series satisfy differential equations which can be mapped into difference equations for the coefficients of the Taylor expansion. In general, not many algorithms exist to solve such a multivariate system of difference equations having the target solution space being that of hypergeometric functions; we refer the reader to Section <ref> for an exposition of one direction more targeted toward the class of rational functions.

However, if a system of differential equations induces a system of recurrences of the form (<ref>), then a method to investigate its solutions is possible using Sigma. Let us concentrate on a system of linear differential equations of the form

\begin{equation}
\Big[
s_i \Bigl(x_1\frac{\partial}{\partial x_i}, \ldots, x_i \frac{\partial}{\partial x_i}, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr)
- \frac{1}{x_i} t_i \Bigl(x_1 \frac{\partial}{\partial x_1}, \ldots, x_i\frac{\partial}{\partial x_i}-1, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr)
\Big] f(x_1,\ldots,x_r)
= 0.
\label{eq:DEsystem}
\end{equation}

One has the property

\begin{equation}
x_i \frac{\partial}{\partial x_i} x_1^{n_1}\cdots x_r^{n_r} = n_i x_1^{n_1}\cdots x_r^{n_r}
\end{equation}

which implies that for a polynomial $p(n_1,\ldots,n_r)$ we have

\begin{equation}
p \Big(x_1\frac{\partial}{\partial x_1}, \ldots, x_i \frac{\partial}{\partial x_i}, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr) x_1^{n_1}\cdots x_r^{n_r} = p(n_1,\ldots,n_r) x_1^{n_1}\cdots x_r^{n_r} .
\end{equation}

\begin{align}
\label{eq:t-term}
&\Big[ s_i \Bigl(x_1\frac{\partial}{\partial x_1}, \ldots, x_i \frac{\partial}{\partial x_i}, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr) \Big] f(x_1,\ldots,x_r)
\nonumber\\
&\qquad \qquad \qquad \qquad \qquad = \sum_{n_i\ge 0} s_i(n_1,\ldots,n_i,\ldots,n_r)  A(n_1,\ldots,n_i,\ldots,n_r) x_1^{n_1}\cdots x_r^{n_r}	\\
&\Big[ t_i \Bigl(x_1 \frac{\partial}{\partial x_1}, \ldots, x_i\frac{\partial}{\partial x_i}-1, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr) \Big] f(x_1,\ldots,x_r)
\nonumber\\
&\qquad \qquad \qquad \qquad \qquad = \sum_{n_i\ge 0} t_i(n_1,\ldots,n_i-1,\ldots,n_r)  A(n_1,\ldots,n_i,\ldots,n_r) x_1^{n_1}\cdots x_r^{n_r}
\end{align}

and therefore, dividing the second equation by $x_i$ from the left,

\begin{align}
\label{eq:s-term}
&\Big[ \frac{1}{x_i} t_i \Bigl(x_1 \frac{\partial}{\partial x_1}, \ldots, x_i\frac{\partial}{\partial x_i}-1, \ldots, x_r \frac{\partial}{\partial x_r} \Bigr) \Big] f(x_1,\ldots,x_r)
\nonumber\\
& \qquad \qquad \qquad \qquad \qquad = \sum_{n_i\ge 0}
t_i(n_1,\ldots,n_i-1,\ldots,n_r)  A(n_1,\ldots,n_i,\ldots,n_r)
x_1^{n_1}\cdots x_i^{n_i-1} \cdots x_r^{n_r}.
\end{align}

The coefficient of the term $x_1^{n_1}\cdots x_i^{n_i}\cdots x_r^{n_r}$ in
(<ref>) and in
(<ref>) is respectively

\begin{equation}
s_i(n_1,\ldots,n_i,\ldots,n_r) A(n_1,\ldots,n_i,\ldots,n_r)
\end{equation}

\begin{equation}
t_i(n_1,\ldots,n_i,\ldots,n_r) A(n_1,\ldots,n_i+1,\ldots,n_r).
\end{equation}

This shows, due to (<ref>), that (<ref>) holds.

For example, for the case of the Gauss hypergeometric function

\begin{equation}
_2F_1(a,b;c;x) = \sum_{n\ge 0} \frac{(a)_n (b)_n}{(c)_n n!}x^n
\end{equation}

one has

\begin{eqnarray}
A(n) &=& \frac{(a)_n (b)_n}{(c)_n n!} \\
s(n) &=& (a+n)(b+n) \\
t(n) &=& (n+1)(c+n)
\end{eqnarray}

and the differential equation obeyed by $_2F_1(a,b;c;x)$ is, from (<ref>),

\begin{equation}
\Big[ \Big(a+ x \frac{\partial}{\partial x} \Big) \Big(b+ x \frac{\partial}{\partial x} \Big) - \frac{1}{x} \Big(x \frac{\partial}{\partial x}\Big)	\Big( x \frac{\partial}{\partial x} -1 +c \Big) \Big] {}_2F_1(a,b;c;x) = 0.
\end{equation}

To summarize, a system of equations of the type (<ref>) induces first-order recurrences, which can be studied using the methods of Sigma as described next.

§.§ An algorithm for hypergeometric products

We consider a system (<ref>) and look for a hypergeometric solution $A$. The univariate case is well-known and treatable under the methods of Sigma. To be explicit, for $r=1$, the polynomials $s_1,t_1$ will be nonzero for all $k>\lambda_1,k,\lambda_1\in \mathbb N$. Then, using the hypergeometric property, one writes

\begin{equation}\label{Equ:Prod1}
\begin{split}
\end{split}
\end{equation}

and a formula for $A$ is obtained as a hypergeometric product $\prod_{k=\lambda_1+1}^{n_1}\frac{s_1(k-1)}{t_1(k-1)}$. For the multivariate case, we seek a generalization and a solution in the form of a nested product.

One calls a sequence non-trivial if its zero points vanish on a polynomial in $\mathbb K[n_1,\ldots,n_r]$, [371]. For such a sequence, the system (<ref>) implies a compatibility condition: calling

\begin{equation}
R_i = \frac{s_i}{t_i} \in K(n_1,\ldots,n_r)
\end{equation}

then for there to be solutions it must also be that <cit.>

\begin{equation}
\frac{R_i(n_1,\dots,n_j+1,\dots,n_r)}{R_i(n_1,\dots,n_j,\dots,n_r)}=\frac{R_j(n_1,\dots,n_i+1,\dots,n_r)}{R_j(n_1,\dots,n_i,\dots,n_r)}.
\label{Equ:CompatibilityProp}
\end{equation}

Relevant to this context, the Ore-Sato theorem [372, 373, 374] states what a general form of the hypergeometric solution, if it exists, must have; it is a product of factorials and hypergeometric terms. Here we specialize to the additional assumption that $s_i,t_i\neq0$ for all $(n_1,\ldots,n_r)\in \mathbb N^r$ with $n_i>\lambda_i$. (In general no algorithm exists [385] to find the $\lambda_i$, although it is often possible in practically occurring cases).

Under these assumptions, one can write

\begin{equation}
A(n_1,\ldots,n_i+1,\ldots,n_r) = R_i(n_1,\dots,n_r)\,A(n_1,\ldots,n_i,\ldots,n_r) , \qquad i=1,\ldots,r
\end{equation}

\begin{equation}
\end{equation}

for any $(n_1,\dots,n_r)\in\mathbb N^r$ with $n_i\geq\lambda_i$ and some $c\in\mathbb K\setminus\{0\}$.

One can derive under these assumptions a formula for $A$ as the nested product

\begin{equation}
c\,\left(\prod_{k=\lambda_1}^{n_1}h_1(k,n_2,\dots,n_r)\right)\left(\prod_{k=\lambda_2}^{n_2}h_2(k,n_3\dots,n_r)\right)\dots \left(\prod_{k=\lambda_r}^{n_r}h_r(k)\right)
\label{Equ:FinalProdForm}
\end{equation}

with $c=A(\lambda_1,\dots,\lambda_r)\in\mathbb K\setminus\{0\}$ in a recursive manner. The case $r=1$ is treated in (<ref>). Otherwise one can write

\begin{equation}
\end{equation}

\begin{equation}
\end{equation}

Then one considers

\begin{equation}
\end{equation}

which satisfies

\begin{equation}
A'(n_2,\ldots,n_i+1,\ldots,n_r) = R_i(\lambda_1,n_2,\dots,n_r)\,A'(n_2,\ldots,n_i,\ldots,n_r) , \qquad i=2,\ldots,r
\end{equation}

with $R_i(\lambda_1,n_2,\dots,n_r)\in\mathbb K[n_2,\dots,n_r]$ and obtains

\begin{equation}
\dots \left(\prod_{k=\lambda_r}^{n_r}h_r(k)\right),
\end{equation}

$c=A'(\lambda_2,\dots,\lambda_r)=A(\lambda_1,\lambda_2,\dots,\lambda_r)\in\mathbb K\setminus\{0\}$
and $h_i(x,n_{i+1},\dots,n_r)\in\mathbb K(x,n_{i},\dots,n_r)$ with $2\leq i\leq r$.

Proceeding in the same way by recursion one obtains (<ref>). We show now the results of the application of the algorithm in some examples.

§.§.§ Examples

We consider the differential equation
(<ref>) which leads to the recurrence (<ref>) for the expansion coefficient $f[n]$. The recursion is
of order one and is solved for $f[n] = 0$. Sigma obtains the product solution

\begin{eqnarray}
f[n] =
\frac{\prod_{i_1=1}^n \big(
-3 i_1
-B_1 i_1
\big)}{n! (A_1)_n}
\equiv \frac{\prod_{i_1=1}^n \big[-C +B_1(1-i_1) + (1-i_1)(2-i_1)
\big]}{n! (A_1)_n},
\label{eq:T1}
\nonumber\\
\end{eqnarray}

One can factorize the product in (<ref>) in terms of Pochhammer symbols by

\begin{eqnarray}
f[n] = \frac{(\alpha_1)_n (\alpha_2)_n}{(A_1)_n n!},
\end{eqnarray}

\begin{eqnarray}
\alpha_{1(2)} = -\frac{1}{2}(1+B_1) \mp  \frac{1}{2} \sqrt{(1+B_1)^2 + 4 C}.
\end{eqnarray}

By replacing $A_1,B_1$ and $C$ by

\begin{eqnarray}
C \rightarrow -a b,~~A_1 \rightarrow c,~~B_1 \rightarrow -1 - a - b
\end{eqnarray}

one obtains

\begin{eqnarray}
f[n] = \frac{(a)_n (b)_n}{(c)_n n!}.
\end{eqnarray}

In the case of the hypergeometric function $_3F_2$, the differential equation (<ref>) implies the recurrence (<ref>) for $f[n]$ with $f[n]=1$, which has the solution

\begin{eqnarray}
f[n] = \frac{\prod_{i_1=1}^n
+B_1 \big(
-B_2 \big(
\big)
{ \big(1-i_1\big) }
\big)
\big(2-i_1
\big)
\big(1-i_1\big)]}
-A_2 \big(
\big)
\big(1-i_1\big)]}.
\label{eq:3F2a}
\end{eqnarray}

Eq. (<ref>) can be rewritten in terms of radicals by

\begin{eqnarray}
f[n] = \frac{
(\alpha_3)_n}{ {n!} \big(
+\frac{A_2 }{2}
\big)_n \big(
+\frac{A_2 }{2}
\big)_n},
\end{eqnarray}

\begin{eqnarray}
\alpha_1 &=&
}{6 \sqrt[3]{2}}
-\frac{i \sqrt[3]{z_1
}}{2 \sqrt[3]{2} \sqrt{3}}
+\frac{\sqrt[3]{-2} z_3}{3 \sqrt[3]{z_1
\\
\alpha_2 &=&
}}{3 \sqrt[3]{2}}
-\frac{\sqrt[3]{2} z_3}{3 \sqrt[3]{z_1
\\
\alpha_3 &=&
}}{6 \sqrt[3]{2}}
+\frac{i \sqrt[3]{z_1
}}{2 \sqrt[3]{2} \sqrt{3}}
-\frac{(-1)^{2/3} \sqrt[3]{2} z_3}{3 \sqrt[3]{z_1
\\
z_1 &=&  27 C + (3 + B_2) (9 B_1 + B_2 (3 + 2 B_2))
\\
z_2 &=&  \sqrt{-4 (3 + 3 B_1 +
B_2 (3 + B_2))^3 + (27 C + (3 + B_2) (9 B_1 +
B_2 (3 + 2 B_2)))^2}
\\
z_3 &=&  3 + 3 B_1 + B_2 (3 + B_2)
\\
z_4 &=&  6 + B_2
\\
z_5 &=&  \sqrt{-4 A_1 + (A_2-1)^2}.
\end{eqnarray}

After performing the replacements

\begin{eqnarray} A_{{2}} &\rightarrow& b_1+b_2 +1,~~ B_{{2}} \rightarrow -(3 + a_1 + a_2 +a_3),~~A_{{1}} \rightarrow b_1 b_2,
\nonumber\\
B_{{1}} &\rightarrow& -(a_2 a_1 + a_3 a_1 + a_2 a_3 + a_1 +
a_2 +a_3 +1),~~C \rightarrow - a_1 a_2 a_3
\label{eq:VIETA}
\end{eqnarray}

one obtains

\begin{eqnarray}
f[n] = \frac{(a_1)_n (a_2)_n (a_3)_n}{(b_1)_n (b_2)_n n!}.
\end{eqnarray}

The replacements are related to Vieta's theorem [375] for the roots $r_i$ of the algebraic equation

\begin{eqnarray}
x^n + \sum_{k = 1}^{n} a_{n-k} x^{n-k} = 0,
\end{eqnarray}

which obey

\begin{eqnarray}
-a_{n-1} &=& r_1 + ... + r_n
\nonumber\\
a_{n-2} &=& r_1 (r_2 + ... + r_n) + r_2(r_3+ ... r_n) + ... r_{n-1} r_n
\nonumber\\
\nonumber\\
(-1)^n a_0 &=& r_1 ... r_n.
\end{eqnarray}

For the function $_{p+1}F_p$ the product solution is

\begin{eqnarray}
f[n] =
\prod_{i=1}^n \frac{\frac{1}{(i-p-2)!} - \sum_{k=1}^p \frac{B_k}{(i-k-1)!} - \frac{C}{(i-1)!}}
{\frac{i}{(i-p-1)!} + \sum_{k=1}^p \frac{A_k i}{(i-k)!} } .
\end{eqnarray}

In the multivariate case, we list the product solutions for the Horn–type functions, from Eqs. (<ref>), (<ref>),

\begin{eqnarray}
f_{\rm H}[m,n] &=&
\Biggl[
\prod_{i_1=1}^m \frac{-a
-2 e
-f n
+h n
+j n
-j n^2
-b i_1
+3 e i_1
-h n i_1
-e i_1^2
+g n
+d i_1
\big) i_1}\Biggr]
\nonumber\\&& \times
\Biggl[\prod_{i_1=1}^n \frac{-a_1
-2 e_1
-b_1 i_1
+3 e_1 i_1
-e_1 i_1^2
+d_1 i_1
\big) i_1}\Biggr]
\label{eq:HORN2}
\end{eqnarray}

and for the functions $S_1$ and $S_2$ one has

\begin{eqnarray}
f_{\rm S_1}[m,n] &=&
\Biggl[
\prod_{i_1=1}^n \frac{- a_1
+ b_1
-2  e_1
- b_1  i_1
+3  e_1  i_1
- e_1  i_1^2
- d_1
+ d_1  i_1
\big) i_1}\Biggr]
\nonumber\\&& \times
\Biggl[
\prod_{i_1=1}^m
\bigg(
\frac{1}{\big(
+2 l
+g n
-n q
+d i_1
-3 l i_1
+n q i_1
+l i_1^2
\big) i_1}	(-a
\nonumber\\&&
-2 e
-f n
+h n
+j n
-j n^2
+6 p
-2 n r
-n s
+n^2 s
-b i_1
\nonumber\\&&
+3 e i_1
-h n i_1
-11 p i_1
+3 n r i_1
+n s i_1
-n^2 s i_1
-e i_1^2
+6 p i_1^2
-n r i_1^2
\nonumber\\&&
-p i_1^3
\bigg)
\bigg]
\\
f_{\rm S_2}[m,n] &=&
\bigg(
\prod_{i_1=1}^n \frac{- a_1
+ b_1
-2  e_1
- b_1  i_1
+3  e_1  i_1
- e_1  i_1^2
+2  p_1
+d_1 i_1
-3  p_1  i_1
+ p_1  i_1^2
\big) i_1}\bigg)
\nonumber\\&& \times
\bigg[
\prod_{i_1=1}^m
\bigg(
\frac{1}{\big(
+2 l
+g n
-n q
+d i_1
-3 l i_1
+n q i_1
+l i_1^2
\big) i_1}
-2 e
-f n
\nonumber\\&&
+h n
+j n
-j n^2
+6 p
-2 n r
-n s
+n^2 s
-b i_1
+3 e i_1
-h n i_1
-11 p i_1
+3 n r i_1
\nonumber\\&&
+n s i_1
-n^2 s i_1
-e i_1^2
+6 p i_1^2
-n r i_1^2
-p i_1^3
\bigg)\Bigr].
\end{eqnarray}

Eq. (<ref>) can be rewritten as

\begin{eqnarray}
f[n,m] &=& \frac{\displaystyle (-1)^{m+n}}
{\displaystyle m! n! \Biggl(
\frac{c_1}{d_1}\Biggr)_n \Biggl(
\frac{c}{d}
+\frac{g n}{d}
\Biggr)_m}
\Biggl(
\frac{e}{d}\Biggr)^m \Biggl(
\frac{e_1}{d_1}\Biggr)^n \Biggl(
+\frac{b_1}{2 e_1}
-\frac{r_1}{2 e_1}
\Biggr)_n
\nonumber\\ && \times
\Biggl(
+\frac{b_1}{2 e_1}
+\frac{r_1}{2 e_1}
\Biggr)_n \Biggl(
+\frac{b}{2 e}
+\frac{h n}{2 e}
-\frac{r_2}{2 e}
\Biggr)_m \Biggl(
+\frac{b}{2 e}
+\frac{h n}{2 e}
+\frac{r_2}{2 e}
\Biggr)_m,
\nonumber\\
\label{eq:FMN}
\end{eqnarray}

\begin{eqnarray}\label{Equ:AlgebraicExt}
r_1 &=& \sqrt{(b_1 - e_1)^2 - 4 a_1 e_1}
\\
r_2 &=&
\sqrt{(b - 3 e + h n)^2 - 4 e (a - b + 2 e + f n - h n - j n + j n^2)}.
\end{eqnarray}

They can be factorized into the usual Pochhammer representation if appropriate replacements are obeyed.

In the tri–variate cases one obtains

\begin{eqnarray}
f[m,n,p] &=&
\bigg(
\prod_{i_1=1}^m \frac{-A
+ B_1
-2  E_1
- B_1 i_1
+3  E_1 i_1
- E_1 i_1^2
+E_0  i_1
\big) i_1}
\bigg)
\nonumber\\&& \times
\bigg[\prod_{i_1=1}^n
\bigg(
\frac{1}{\big(
- F'_0
+ H'_2  m
+ F'_0  i_1
\big) i_1}
(- A'
+ C'_1
-2  F'_1
- B'_1  m
\nonumber\\&&
+ E'_1  m
+ H'_1  m
- E'_1  m^2
- C'_1  i_1
+3  F'_1  i_1
- H'_1  m i_1
- F'_1  i_1^2
\bigg)
\bigg]
\nonumber\\&& \times
\bigg[
\prod_{i_1=1}^p
\bigg(
\frac{1}{\big(
- G''_0
+ L''_2  m
+n  S''_2
+ G''_0  i_1
\big) i_1}
+ D''_1
-2  G''_1
- B''_1  m
+ E''_1  m
\nonumber\\&&
+ L''_1  m
- E''_1  m^2
- C''_1  n
+ F''_1  n
- H''_1  m n
- F''_1  n^2
+n  S''_1
- D''_1  i_1
+3  G''_1  i_1
- L''_1  m i_1
\nonumber\\&&
-n  S''_1  i_1
- G''_1  i_1^2
\bigg)
\bigg].
\end{eqnarray}

Finally, in the four–variable case the product solution reads

\begin{eqnarray}
f[m,n,p,q] &=&	\bigg(
\prod_{i_1=1}^m
\frac{-A
-2 F_1
-B_1 i_1
+3 F_1 i_1
-F_1 i_1^2
+F_0 i_1
\big) i_1}
\bigg)
\nonumber\\&&\times
\bigg[
\prod_{i_1=1}^n
\frac{1}{\big(
+m M'_2
+G'_0 i_1
\big) i_1}
-2 G'_1
-B'_1 m
+F'_1 m
\nonumber\\&&
-F'_1 m^2
+m M'_1
-C'_1 i_1
+3 G'_1 i_1
-m M'_1 i_1
-G'_1 i_1^2
\bigg]
\nonumber\\&&\times
\bigg[
\prod_{i_1=1}^p
\frac{1}{\big(
+m N''_2
+n Q''_2
+H''_0 i_1
\big) i_1}
-2 H''_1
-B''_1 m
+F''_1 m
\nonumber\\&&
-F''_1 m^2
-C''_1 n
+G''_1 n
-m M''_1 n
-G''_1 n^2
+m N''_1
+n Q''_1
-D''_1 i_1
+3 H''_1 i_1
\nonumber\\&&
-m N''_1 i_1
-n Q''_1 i_1
-H''_1 i_1^2
\bigg]
\nonumber\\&&\times
\bigg[
\prod_{i_1=1}^q \frac{1}{\big(
+m P'''_2
+n R'''_2
+p S'''_2
+L'''_0 i_1
\big) i_1}
-2 L'''_1
-B'''_1 m
\nonumber\\&&
+F'''_1 m
-F'''_1 m^2
-C'''_1 n
+G'''_1 n
-m M'''_1 n
-G'''_1 n^2
-D'''_1 p
+H'''_1 p
-m N'''_1 p
\nonumber\\&&
-H'''_1 p^2
+m P'''_1
-n p Q'''_1
+n R'''_1
+p S'''_1
-E'''_1 i_1
+3 L'''_1 i_1
-m P'''_1 i_1
-n R'''_1 i_1
\nonumber\\&&
-p S'''_1 i_1
-L'''_1 i_1^2
\bigg].
\end{eqnarray}

§.§ Computing the expansion in $\ep$

Algorithms for the series expansion of hypergeometric series would be of interest for their applicability to physics. One may try to obtain the $\ep$-expansion of the hypergeometric series as nested sums using the difference ring methods of EvaluateMultiSums, but this is not always possible. An alternative is to try to obtain the series expansion of the summand and to sum the terms of the expansion. For an explanation of an algorithms which can work on nested products of the type discussed in Section <ref>, we refer to <cit.>. For these methods to be applicable, the convergence region of the hypergeometric function needs to be understood, as well as the conditions under which the summation quantifier commutes with the differential operator in $\ep$.

Here we reproduce from [201] two examples in which these methods are employed to obtain the $\ep$ expansion of hypergeometric functions.

§.§.§ Example 1

Consider for example the system of equations

\begin{eqnarray}
\Bigg[ (x-1) y \partial_{x,y}^2+ \Big[x \Big(2 \varepsilon +\frac{7}{2}\Big)-\varepsilon +1\Big]\partial_x + (x-1) x \partial_x^2
+y (2 \varepsilon +1) \partial_y +\frac{3}{2} (2 \varepsilon +1) \Bigg] f(x,y) &= 0, \\
\Bigg[ x (y-1) \partial_{x,y}^2 +x (4-\varepsilon ) \partial_x+ \Big[y \Big(\frac{13}{2}-\varepsilon \Big)-\varepsilon +1\Big]\partial_y
+(y-1) y \partial_y^2+\frac{3 (4-\varepsilon )}{2}
\Bigg]f(x,y) &= 0,
\end{eqnarray}

for which we search for a solution of the form (<ref>) with $r=2$ where $x_1=x$ and $x_2=y$.
Computing a first–order recurrence system of $A(n_1,n_2)=A(m,n)$ and solving it by the method presented in Section <ref> provides the solution

\begin{equation}
f(x,y) = \sum_{m,n= 0}^\infty A(m,n) = \sum_{m,n= 0}^\infty \frac{x^m y^n \big(
\frac{3}{2}\big)_{m
} (4-\varepsilon )_n (1+2 \varepsilon )_m}{m! n! (-1+\varepsilon )_{m
\label{eq:example-summand}
\end{equation}

A series expansion of the summand $A(m,n)$ in (<ref>) up to $\mathcal O(\varepsilon^0)$ gives

\begin{eqnarray}
A(m,n)&=&	-\frac{1}{6} \frac{x^m y^n (3+n)! \big(
\frac{3}{2}\big)_{m
}}{n! (-2
)! \varepsilon}
\frac{1}{36} \bigg[
-\frac{1}{(1+n) (2+n) (3+n) (m
) (-1
\nonumber\\&& \times
\big(
-30 n
+17 n^2
+97 n^3
+79 n^4
+17 n^5
+m^2 \big(
36+115 n+84 n^2+17 n^3\big)
\nonumber\\&&
+m \big(
36+89 n
+218 n^2
+163 n^3+34 n^4\big)
\big)
-12 S_1({m})
+6 S_1({n})
+6 S_1({m+n})
\bigg]
\nonumber\\&& \times
\frac{x^m y^n (3+n)! \big(
\frac{3}{2}\big)_{m
}}{n! (-2
+\mathcal O(\varepsilon).
\end{eqnarray}

A series expansion of (<ref>) in the region $0<x<\sqrt{y}$, $0<y<\frac{1}{2}$,

\begin{equation}
f(x,y) = \frac{1}{\varepsilon} f_{-1}(x,y) + f_0(x,y) +\mathcal O(\varepsilon)
\end{equation}

is possible using EvaluateMultiSum and results in an expression involving the sums

\begin{eqnarray}
R_0 &=& \sum_{i=1}^{\infty } \frac{x^i \big(\frac{3}{2}\big)_i}{i!} = -1+\frac{1}{(1-x)^{3/2}}
\\
R_1 &=& \sum_{i=1}^{\infty } \frac{y^i \big(\frac{3}{2}\big)_i}{i!} = -1+\frac{1}{(1-y)^{3/2}}
\end{eqnarray}

at $\mathcal O(\varepsilon^{-1})$. The function $f_{-1}(x,y)$ reads

\begin{eqnarray}
f_{-1}(x,y) &=& -\frac{15 x^6 }{4 (x-y)^4(1-x)^{7/2}}
-\frac{15 y^3}{64 (x-y)^4 (1-y)^{13/2}} \big[
y^3 \big(
160+80 y-10 y^2+y^3\big)
\nonumber\\&&
-x y^2 \big(
576+176 y
-64 y^2+5 y^3\big)
+x^3 \big(
-320+120 y-36 y^2+5 y^3\big)
\nonumber\\&&
+3 x^2 y \big(
240+8 y-22 y^2+5 y^3\big)
\big] .
\end{eqnarray}

In addition, one encounters at $\mathcal O(\varepsilon^0)$ the sums

\begin{eqnarray}
% R2
R_2 &=& \sum_{i=1}^{\infty } \frac{x^i \big(
\frac{3}{2}\big)_i}{(1+2 i)^2 i!} = -1 + \frac{1}{\sqrt{x}} \arcsin\big( \sqrt{x}\big)
\\ % R3
R_3 &=& \sum_{i=1}^{\infty } \frac{y^i \big(
\frac{3}{2}\big)_i}{i i!} = -2
+2 \ln(2)
+2 \frac{1}{\sqrt{1-y}}
-2 \HA_{-1}\big(
\sqrt{1-y}\big)
\\ % R4
R_4 &=& \sum_{i_1=1}^{\infty } \frac{x^{i_1} \big(\frac{3}{2}\big)_{i_1}
\sum_{i_2=1}^{i_1} \frac{1}{1+2 i_2}
= \frac{1}{2} \frac{\HA_1(x)}{(1-x)^{3/2}}
\\ % R5
R_5 &=& \sum_{i_1=1}^{\infty } \frac{x^{i_1}
\big(\frac{3}{2}\big)_{i_1}}{i_1!}
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
\nonumber\\&=&
\frac{y}{(1-x) (y-x)}
+\frac{y}{2(1-x)^{3/2} \sqrt{1-y}} \Big[
i \pi
-\HA_0\big( \sqrt{1-y} -\sqrt{1-x} \big)
\nonumber\\&&
-2 \HA_{-1}\big( \sqrt{1-y}\big)
+\HA_0\big( \sqrt{1-y}+\sqrt{1-x} \big)
\Big]
\\ % R6
R_6 &=& \sum_{i_1=1}^{\infty } \frac{x^{i_1} y^{i_1}
\big(\big(
\frac{3}{2}\big)_{i_1}\big)^2}{\big(
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
\nonumber\\&=&
\frac{1}{2}\int_0^1 dt \bigg[
\frac{-1+t}{\pi  (-1
)} \Big[
\frac{1}{(1
-(1-t) x
\big[
4 E(x -t x)
-2 [1
-(1-t) x
] K(x-t x)
\big]
\nonumber\\&&
-\frac{4 E(x y)
+2 (-1
+x y
) K(x y)
+x y
\Big] \frac{1}{\sqrt{t}}
\bigg]
\\ % R7
R_7 &=& \sum_{i_1=1}^{\infty } \frac{x^{i_1} y^{i_1}
\big(\big(
\frac{3}{2}\big)_{i_1}\big)^2}{\big(
i_1!\big)^2 \big(
1+2 i_1\big)^2}
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
{ \frac{1}{\pi} } \int_0^1 dt \frac{t-1}{\sqrt{t}(t+y-1) } \big[ K(x(1-t))
-K(x y) \big]
\nonumber\\&&
\\ % R8
R_8 &=& \sum_{i_1=1}^{\infty } \frac{y^{i_1}
\big(\frac{3}{2}\big)_{i_1}}{i_1! \big(
1+2 i_1\big)}
\sum_{i_2=1}^{i_1} \frac{x^{i_2} y^{-i_2}}{i_2}
-2 \frac{\HA_0\big(
\sqrt{1-x}+\sqrt{1-y}\big)}{\sqrt{1-y}}
+2 \frac{\HA_{-1}\big(
\sqrt{1-y}\big)}{\sqrt{1-y}}
\\ % R9
R_9 &=& \sum_{i_1=1}^{\infty } \frac{x^{i_1}
\big(\frac{3}{2}\big)_{i_1}}{i_1!} \big(
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
\big)
\big(
\sum_{i_2=1}^{i_1} \frac{y^{i_2} \big(
\frac{3}{2}\big)_{i_2}}{i_2!}
\big)
\nonumber\\&=&
\frac{1}{(1-y)^{5/2}} \Big[ \left(\frac{1}{(1-x)^{3/2}}-1\right) \left((1-y)^{3/2}-1\right)
\biggl(\sqrt{1-y} y
\Big(\HA_{-1}\big(\sqrt{1-y}\big)
\nonumber\\&&
-\frac{\HA_0(y)}{2}+\frac{i \pi }{2}\Big)-y+1\biggr) \Big]
+\sum_{i_ 1=1}^{\infty } \bigg\{
\frac{1}{\pi
(1-y)^2 \sqrt{1- y} \Gamma \big(
1+i_ 1\big) \Gamma \big(
2+i_ 1\big)}
\nonumber\\&& \times
\bigg[ 4 x^{i_ 1} y^{1+i_ 1} \Big[
+\sqrt{1-y} y \Big(
\frac{i \pi }{2}
-\frac{1}{2} \HA_ 0(y)
+\HA_ {-1}\big(
\sqrt{1-y}\big)
\Big)
\Big] \Gamma \big(
\frac{3}{2}+i_ 1\big)
\nonumber\\&& \times
\Gamma \big(
\frac{5}{2}+i_ 1\big)
\,_ 2F_ 1\Big(
-\frac{1}{2},1+i_ 1;2+i_ 1;y\Big) \bigg]
\sqrt{\pi -\pi  y} \Gamma \big(
1+i_ 1\big)}
\nonumber\\&& \times
\Big[2 x^{i_ 1} \Gamma \big(
\frac{3}{2}+i_ 1\big)
\,_ 2F_ 1\Big(
-\frac{1}{2},1+i_ 1;2+i_ 1;y\Big)
\,_ 2F_ 1\Big(
1,2+i_ 1;\frac{5}{2}+i_ 1;\frac{1}{y}\Big) \Big]
\nonumber\\&&
-\frac{1}{(-1+y) \sqrt{\pi -\pi  y} \big(
3+2 i_ 1\big)}
\Big[ 2 \sqrt{\pi } x^{i_ 1} y^{-1-i_ 1} \big(
\nonumber\\&& \times
\,_ 2F_ 1\Big(
_ 1;\frac{5}{2}+i_ 1;\frac{1}{y}
\Big)
\big(1+i_ 1\big) \Big]
\bigg\}
\\ % R10
R_{10} &=& \sum_{i_1=1}^{\infty } \frac{y^{i_1} \big(
\frac{3}{2}\big)_{i_1} S_1\big({i_1}\big) i_1}{i_1!}
-3 \ln(2) y \frac{1}{(1-y)^{5/2}}
+\frac{3}{2} y \frac{\HA_1(y)}{(1-y)^{5/2}}
\nonumber\\&&
+3 y \frac{\HA_{-1}\big(
\sqrt{1-y}\big)}{(1-y)^{5/2}}
+y \big(
3-2 \sqrt{1
\Big] \frac{1}{(1-y)^{5/2}},
\end{eqnarray}

as well as the combination

\begin{eqnarray}
R_{11} &=&\big( 1-\big(1-x\big)^{3/2} \big)
\sum_{i_1=1}^{\infty } \frac{y^{i_1}
\big(\frac{3}{2}\big)_{i_1}}{i_1! \big(
1+2 i_1\big)}
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
\nonumber\\&&
\sum_{i_1=1}^{\infty } \frac{y^{i_1}
\big(\frac{3}{2}\big)_{i_1}}{i_1! \big(
1+2 i_1\big)}\Big(
\sum_{i_2=1}^{i_1} \frac{y^{-i_2} i_2!}{\big(
\frac{3}{2}\big)_{i_2}}
\Big)
\Big(
\sum_{i_2=1}^{i_1} \frac{x^{i_2} \big(
\frac{3}{2}\big)_{i_2}}{i_2!}
\Big)
\nonumber\\&=&
\frac{1}{4} \bigg\{
\Big[
-2 \sqrt{1-x}
+2 x \sqrt{1-x}
-3 x F_1\Big(
\frac{5}{2};\frac{1}{2},1;2;x y,x
\Big)
\big(1-x\big)^{3/2}
\Big]
\big[i \pi  y
\nonumber\\&&
+2 \sqrt{1-y}
-y \HA_0(y)
+2 y \HA_{-1} \big(\sqrt{1-y}\big)
\big]\bigg\} \frac{1}{\sqrt{1-y}}
\nonumber\\&&
\sum_{i_1=1}^{\infty } \Big[ \frac{x^{1+ i_1} \big(
1-x\big)^{3/2} \Gamma \big(
\frac{1}{2}+i_1\big)  }{\sqrt{\pi } y \Gamma \big(
\,_2F_1\Big(
\nonumber\\&& \times
\,_2F_1\Big(
1,\frac{5}{2}+i_1;2+i_1;x\Big) \Big].
\end{eqnarray}

The harmonic polylogarithms [266] are defined by Eqs. (<ref>)-(<ref>).

One can further employ the relations

\begin{eqnarray}
_2F_1\Big(\frac{3}{2},\frac{3}{2};1,z\Big) &=& \frac{2(z-1)K(z)+4E(z)}{\pi(z-1)^2}
\\
K(z) &=& \int_0^1 \frac{1}{\sqrt{(1-t^2)(1-z t^2)}} dt = \frac{\pi}{2} \,_2F_1\Bigr(\frac{1}{2},\frac{1}{2};1;z  \Bigl)
\\
E(z) &=& \int_0^1 \frac{\sqrt{1-z t^2}}{\sqrt{1-t^2}} dt = \frac{\pi}{2} \,_2F_1\Bigl(-\frac{1}{2},\frac{1}{2},1,z  \Bigr).
\end{eqnarray}

The function $f_0(x,y)$ reads

\begin{eqnarray}
f_0(x,y) &=&
\frac{5 R_{10} y^2}{16 (1-y)^4 (x
)^4} \Bigl[
y^3 \big(
160+80 y-10 y^2+y^3\big)
-x y^2 \big(
576+176 y-64 y^2+5 y^3\big)
\nonumber\\&&
+x^3 \big(
-320+120 y-36 y^2+5 y^3\big)
+3 x^2 y \big(
240+8 y-22 y^2+5 y^3\big)
\Bigr]
\nonumber\\&&
-\frac{15 R_8 y^3}{32 (1-y)^6 (x
)^4} \Bigl[
y^3 \big(
160+80 y-10 y^2+y^3\big)
-x y^2 \big(
576+176 y-64 y^2+5 y^3\big)
\nonumber\\&&
+x^3 \big(
-320+120 y-36 y^2+5 y^3\big)
+3 x^2 y \big(
240+8 y-22 y^2+5 y^3\big)
\Bigr]
\nonumber\\&&
+\frac{1}{128 (1-x)^2 (1-y)^6 y (x
)^4} \big(
-960 x^7 (-1+y)^7
+y^5 \big(
128-1344 y+1536 y^2
\nonumber\\&&
-4240 y^3+4110 y^4-223 y^5+33
+8 x^6 y \big(
4-24 y+60 y^2-4880 y^3+1860 y^4
\nonumber\\&&
-564 y^5+79 y^6\big)
+x^5 y \big(
-832+2752 y-1920 y^2+76160 y^3+69280 y^4-8772
\nonumber\\&&
+2427 y^6-495 y^7\big)
-x y^4 \big(
512-4544 y+3264 y^2-34480 y^3+4240 y^4+3363 y^5
\nonumber\\&&
y^6+66 y^7\big)
+x^3 y^2 \big(
-512+384 y+11008 y^2+83680 y^3+169980 y^4+6287
\nonumber\\&&
+5931 y^6+747 y^7-305 y^8\big)
+x^2 y^3 \big(
768-4736 y-2272 y^2-84288 y^3-56570 y^4
\nonumber\\&&
y^5-3549 y^6+387 y^7+33 y^8\big)
+x^4 y \big(
128+1984 y-9792 y^2-29440 y^3
\nonumber\\&&
-180320 y^4-63768
y^5+10536 y^6-8583 y^7+2055 y^8\big)
\big)
\nonumber\\&&
\biggl[
\frac{1}{128 (1-x)^2 (1-y)^5 y (x
)^5} \Bigl[
960 x^7 (-1+y)^6
-y^6 \big(
128+4800 y-4640 y^2
\nonumber\\&&
+4480 y^3-270 y^4+33
+2 x y^5 \big(
320+10944 y-4016 y^2+4664 y^3+1749 y^4
\nonumber\\&&
y^5+33 y^6\big)
+x^6 y \big(
-448+1920 y-12800 y^2+5440 y^3+1920 y^4-628
y^5+65 y^6\big)
\nonumber\\&&
-x^2 y^4 \big(
1280+38080 y+20128 y^2-6304 y^3+18546
y^4-4204 y^5+406 y^6+33 y^7\big)
\nonumber\\&&
+2 x^3 y^3 \big(
640+15040 y+32080 y^2-9896 y^3+11559 y^4-3791
y^5-491 y^6+169 y^7\big)
\nonumber\\&&
-5 x^4 y^2 \big(
128+1856 y+11712 y^2+1856 y^3-2076 y^4+1727
y^5-2058 y^6+448 y^7\big)
\nonumber\\&&
+2 x^5 y \big(
64+320 y+8000 y^2+15360 y^3-12080 y^4+5164
y^5-4170 y^6+935 y^7\big)
\Bigr]
\nonumber\\&&
-\frac{15 R_5 x^6 (1-y)^2}{2 (1-x)^2 y (x
\biggr]
+R_0 \biggl[
\frac{1}{16 (1-x)^2 (1-y)^6 y (x
)^5} \Bigl[
-120 x^8 (-1+y)^7
\nonumber\\&&
+y^6 \big(
-32+384 y+2988 y^2+140 y^3-15 y^4\big)
+5 x y^5 \big(
32-352 y-2676 y^2-1772 y^3
\nonumber\\&&
-95 y^4+12 y^5\big)
+5 x^2 y^4 \big(
-64+608 y+4688 y^2+7516 y^3+1723 y^4+100
y^5-18 y^6\big)
\nonumber\\&&
+x^6 y \big(
20-624 y+3324 y^2+8040 y^3+18380 y^4-6720
y^5+2164 y^6-329 y^7\big)
\nonumber\\&&
+5 x^3 y^3 \big(
64-448 y-4056 y^2-12396 y^3-6809 y^4-612
y^5-10 y^6+12 y^7\big)
\nonumber\\&&
-5 x^4 y^2 \big(
32-64 y-1888 y^2-9120 y^3-11664 y^4-1436
y^5-158 y^6+40 y^7+3 y^8\big)
\nonumber\\&&
+x^5 y \big(
32+416 y-2848 y^2-12000 y^3-46800 y^4-11880
y^5+430 y^6-200 y^7+85 y^8\big)
\nonumber\\&&
+x^7 \big(
-60+304 y-444 y^2-360 y^3-2780 y^4-1080
y^5+1536 y^6-701 y^7+120 y^8\big)
\Bigr]
\nonumber\\&&
+\frac{15 R_3
x^6}{4 (1-x)^2 (x
-\frac{15 R_1 x^6 (1-y)}{2 (1-x)^2 y (x
\biggr]
+\frac{15 R_9 x^6 (1-y)^2}{2 (1-x)^2 y (x
\nonumber\\&&
+\frac{15 R_6 x^6 (1-y) (-1
+(2+x) y
)}{4 (1-x)^2 y (x
+\frac{15 R_3 x^6}{4 (1-x)^2 (x
+\frac{15 R_2 x^6 (1-y)}{4 (1-x)^2 y (x
\nonumber\\&&
+\frac{15 R_4 x^6 (1-y)}{2 (1-x)^2 y (x
-\frac{15 R_7 x^6 (1-y)}{4 (1-x)^2 y (x
-\frac{15 R_{11} x^6 (1-y) }{2 y (x
-y)^4 (1-x)^{7/2}}.
\end{eqnarray}

The sums $R_i$ could be treated using the methods of [252] which are encoded in HarmonicSums.

§.§.§ Example 2

Consider for example the system of equations

\begin{eqnarray}
1+\varepsilon +(2-x+\varepsilon ) \partial_x+2 x (1+x)
\partial_x^2=0\\
2-\varepsilon +(1-2 y+2 \varepsilon ) \partial_y+y (3+y) \partial_y^2=0 ~.
\end{eqnarray}

We can write its solution as

\begin{equation}
\mathcal F(x,y) = \sum_{x,y\ge 0}A(m,n) x^m y^n
\end{equation}

\begin{equation}
A(m,n) = \bigg(
\prod_{i_1=1}^m \frac{-6
+7 i_1
-2 i_1^2
\varepsilon
+2 i_1
\big) i_1}\bigg) \prod_{i_1=1}^n \frac{-6
+5 i_1
+2 \varepsilon
+3 i_1
\big) i_1}.
\end{equation}

The quantity $A(m,n)$ can also be expressed as

\begin{eqnarray}
A(m,n) &=& \frac{(-1)^m  \big(
-\frac{3}{4}-\frac{1}{4} \sqrt{1
-8 \varepsilon
}\big)_m \big(
\frac{1}{4} \big(
-8 \varepsilon
}\big)\big)_m }{\big(
1+\frac{\varepsilon }{2}\big)_m  \Gamma (1+m) }
\nonumber\\&&\times
\frac{ (-1)^n 3^{-n} \big(
-\frac{3}{2}-\frac{1}{2} \sqrt{1
+4 \varepsilon
}\big)_n \big(
\frac{1}{2} \big(
+4 \varepsilon
\frac{1}{3}+\frac{2 \varepsilon }{3}\big)_n \Gamma (1+n)}
\end{eqnarray}

and $\mathcal F(x,y)$ can be rewritten as

\begin{eqnarray}
\mathcal F(x,y) &=& \Big( \sum_{m\ge 0} x^m f_1(m,\varepsilon) \Big) \Big( \sum_{n\ge 0}y^n f_2(n,\varepsilon) \Big).
\\
&=& F_1(x,\varepsilon) F_2(y,\varepsilon).
\end{eqnarray}

Expanding $F_1$ and $F_2$ in a series in $\varepsilon$ using EvaluateMultiSums, one can write an expression containing infinite (nested) sums. These are rewritten as iterated integrals following [252]. Two of the sums are written in semi-analytic form as definite integrals by writing part of the summand as the Mellin transform of a function. For example, we encounter the sum

\begin{equation}
s_1=\sum_{i=1}^{\infty } \frac{(-1)^i x^i \big(
-\frac{3}{2}+i\big)! \big(
\sum_{j=1}^i \frac{1}{1+2 j}\big) S_1({i})}{i i!}.
\end{equation}

By isolating the term $i=1$ and applying the Legendre duplication formula

\begin{equation}
\Gamma\Big(z+\frac{1}{2}\Big) = \sqrt{\pi} \frac{\Gamma(2z)}{2^{2z-1} \Gamma(z)}
\end{equation}

and the identity

\begin{equation}
\Gamma(2z) = \frac{1}{2} \binom{2z}{z} \Gamma(z) \Gamma(z+1)
\end{equation}

we write

\begin{eqnarray}
s_1 &=& -\frac{1}{3} x \sqrt{\pi }
\sum_{i=1}^{\infty } \frac{(-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i}
\binom{2 i}{i}}{(1+i)^3 (3+2 i)}
\sum_{i=1}^{\infty } \frac{(-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i}
\binom{2 i}{i}
\sum_{j=1}^i \frac{1}{1+2 j}}{(1+i)^3}
\nonumber\\&&
\sum_{i=1}^{\infty } \frac{(-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i}
\binom{2 i}{i} S_1({i})}{(1+i)^2 (3+2 i)}
\nonumber\\&&
\sum_{i=1}^{\infty } \frac{(-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i}
\binom{2 i}{i} \big(
\sum_{j=1}^i \frac{1}{1+2 j}\big) S_1({i})}{(1+i)^2}.
\end{eqnarray}

The first three sums are treated following [252]. The fourth sum can be written as

\begin{eqnarray}
t_1 &=& \sum_{i=1}^{\infty } \frac{(-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i}
\binom{2 i}{i} \big(
\sum_{j=1}^i \frac{1}{1+2 j}\big) S_1({i})}{(1+i)^2}
\nonumber\\&=&
\sum_{i=1}^\infty \frac{(-1)^{1+i}2^{-2i}\sqrt{\pi} {x^{1+i}} \binom{2i}{i}}{(1+i)^2} \int_0^1 dz \bigg\{ (z^i-1) \frac{1}{2 (-1+z)} \bigg[
+2 z
1+\sqrt{z}\big) G\Big(
\frac{\sqrt{\tau }}{1-\tau };z\Big)
\nonumber\\&&
+2 \sqrt{z} \big(1-\ln (2)\big)
\bigg] \bigg\}
\nonumber\\&=&
\int_0^1 dz \bigg\{ { \frac{1}{2(z-1)} } \bigg[
+2 z
1+\sqrt{z}\big) G\Big(
\frac{\sqrt{\tau }}{1-\tau };z\Big)
+2 \sqrt{z} \big(1-\ln (2)\big)
\bigg]
\nonumber\\&& \times
\sum_{i=1}^\infty \frac{\big(
-1+z^i\big) (-1)^{1+i} 2^{-2 i} \sqrt{\pi } x^{1+i} \binom{2 i}{i}}{(1+i)^2}
\bigg\}
\nonumber\\&=&
\int_0^1 dz \bigg\{ { \frac{1}{2(z-1)} }
\bigg[
+2 z
1+\sqrt{z}\big) G\Big(
\frac{\sqrt{\tau }}{1-\tau };z\Big)
+2 \sqrt{z} \big(1-\ln (2)\big)
\bigg]
\bigg[
-\frac{4}{z} \Big[
\nonumber\\&&
+x z
-z \sqrt{1+x}
+z \HA_0\Big(
\frac{1}{2} \big(
\frac{1}{2} \Big(
1+\sqrt{1+x z}\big)\Big)
\Big] \sqrt{\pi }
\bigg]
\bigg\}.
\end{eqnarray}

The $\varepsilon$ expansion of $F_1(x,\varepsilon)$ and $F_2(x,\varepsilon)$ then can be written by

\begin{eqnarray}
F_1(x,\varepsilon) &=& 1
\frac{x}{2}
+\varepsilon  \biggl\{
+\frac{1}{4} x \big(
-9+4 \sqrt{1
+\frac{1}{2} (-2+x) \HA_0(x)
\nonumber\\&&
+\frac{1}{2} (2-x)  G_{3}(x)
\biggr\}
+\varepsilon ^2 \Biggl\{
\frac{1}{8} \Bigl[
20 \big(
+x \big(
+4 x
+20 \sqrt{1+x}
\big)
\Bigr]
\nonumber\\&&
+(-2+x) \Bigl(
-\frac{1}{4} \HA_0(x)^2
+\frac{1}{4} \HA_{0,-1}(x)
+\frac{ G_{11}(x)}{4}
-\frac{ G_{12}(x)}{4}
+\frac{ G_{5}(x)^2}{8}
\Bigr)
\nonumber\\&&
+\frac{1}{2} (2-x) \bigl( G_{8}(x)
+ G_{9}(x)
\bigr)
+\frac{1}{4} (-4+13 x) \HA_0(x)
+\frac{1}{2} (1+x)^{3/2} \HA_{-1}(x)
\nonumber\\&&
-\frac{13 x}{4}
+\frac{1}{2} (-2+x) \HA_0(x)
\biggr]  G_{3}(x)
\frac{1}{2} (1+x)^{3/2}
+\frac{1}{4} (2-x) \HA_0(x)
\biggr]  G_{5}(x)
\Biggr\}
\nonumber\\&&
+\varepsilon ^3 \Biggl\{
\frac{ G_{5}(x)^2}{16} \Bigl[
x \big(
17-4 \sqrt{1
-4 \big(
\Bigr]
+ G_{11}(x) \Bigl[
\frac{1}{8} \big(
x \big(
13-6 \sqrt{1
\nonumber\\&&
-6 \sqrt{1+x}
\big)
+\frac{1}{8} (-2+x) \HA_0(x)
\Bigr]
+ G_{12}(x) \biggl[
\frac{1}{8} \Bigl[
2 \big(
+x \big(
-17+2 \sqrt{1
\Bigr]
\nonumber\\&&
+\frac{3}{8} (-2+x) \HA_0(x)
\biggr]
+ G_{3}(x) \biggl[
\frac{9}{2}
-\frac{105 x}{8}
+\frac{1}{4} (-14+17 x) \HA_0(x)
+\frac{1}{4} (2-x) \HA_0(x)^2
\biggr]
\nonumber\\&&
+ G_{5}(x) \Biggl[
(-2+x) \Bigl(
-\frac{1}{16} \HA_0(x)^2
+\frac{1}{8} \HA_{0,-1}(x)
+\frac{3  G_{11}(x)}{8}
-\frac{ G_{12}(x)}{8}
\Bigr)
\nonumber\\&&
+\frac{1}{4} \biggl[
-6 x
-2 x^2
+(1+x) (5+16 x) \sqrt{1+x}
\biggr]
+\frac{1}{8} \biggl[
x \big(
-13+6 \sqrt{1
\nonumber\\&&
+6 \sqrt{1+x}
\biggr] \HA_0(x)
-\frac{1}{4} (1+x)^{3/2} \HA_{-1}(x)
\Biggr]
+(-2+x) \big(
\frac{1}{12} \HA_0(x)^3
-\frac{3}{8} \HA_{0,0,-1}(x)
\nonumber\\&&
-\frac{1}{8} \HA_{0,-1,-1}(x)
-\frac{ G_{18}(x)}{8}
-\frac{5  G_{19}(x)}{8}
-\frac{3  G_{20}(x)}{8}
+\frac{ G_{21}(x)}{8}
-\frac{3  G_{22}(x)}{4}
+\frac{ G_{23}(x)}{4}
\nonumber\\&&
-\frac{1}{24}  G_{5}(x)^3
\big)
+\frac{1}{2} (2-x) ( G_{14}(x)
+ G_{15}(x)
+ G_{16}(x)
+ G_{17}(x)
\nonumber\\&&
+\frac{1}{240} \Biggl\{
-60 (-2+x)
\int_0^1 dz \biggl\{
-\frac{1}{(-1+z) z}2 \sqrt{\pi } \Bigl[
+2 z
1+\sqrt{z}\big) G_1(z)
\nonumber\\&&
-2 \sqrt{z} \bigl(-1+\ln (2)\bigr)
\Bigr]
\biggl[-1
-\sqrt{1+x} z
+x z
+z \HA_0\Big(
\frac{1}{2} \big(
\nonumber\\&&
\frac{1}{2} \big(
\biggr] \biggr\}
+180 x
\int_0^1 dz \biggl\{
\frac{1}{4 \sqrt{1
+x z
}}\sqrt{\pi } x \big(
+x z
\big)
\biggl[4
\nonumber\\&&
+2 \sqrt{z} \Bigl[
+2 z^{3/2}
-3 \sqrt{z} \big(
\big)
+z \big(6-4 \ln (2)\big)
+8 \ln (2)
\Bigr]
+4 \HA_0(z)
\nonumber\\&&
+4 \HA_1(z)
+2 \Bigl[
-4 \sqrt{z}
+3 z
+2 z^{3/2}
+(-1+z) \HA_0(z)
+(-1+z) \HA_1(z)
+2 \ln (2)
\nonumber\\&&
-2 z \ln (2)
\Bigr]  G_{1}(z)
+(-1+z)  G_{1}(z)^2
-2 (-1+z)  G_{6}(z)
-2 (-1+z)  G_{7}(z)
+6 \zeta_2
\biggr]\biggr\}
\nonumber\\&&
x \Bigl[
+4 x \big(
+96 x
-500 \sqrt{1+x}
\big)
+60 \sqrt{1+x}
\Bigr]
\nonumber\\&&
+2060 \big(
\biggr] \sqrt{\pi }
\Biggr\} \frac{1}{\sqrt{\pi }}
+\frac{3}{8} (-12+35 x) \HA_0(x)
+\frac{1}{8} (14-17 x) \HA_0(x)^2
\nonumber\\&&
+\frac{1}{4} (11-16 x) (1+x)^{3/2} \HA_{-1}(x)
+\frac{5}{8} (1+x)^{3/2} \HA_{-1}(x)^2
+\frac{1}{8} \Bigl[
x \big(
17-2 \sqrt{1
\nonumber\\&&
-2 \big(
\Bigr] \HA_{0,-1}(x)
+8 \sqrt{x}  G_{10}(x)
+8 \sqrt{x}  G_{13}(x)
+\frac{ G_{2}(x)}{2}
\nonumber\\&&
2 (17-4 x) \sqrt{x}
-8 \sqrt{x}  G_{5}(x)
\Bigr]  G_{4}(x)
+2 (-9+4 x) \sqrt{x}  G_{4}(x)
\frac{7}{2}
-\frac{17 x}{4}
\nonumber\\&&
+\frac{1}{2} (-2+x) \HA_0(x)
\biggr]  G_{8}(x)
-\frac{7}{2} (-1+x)
+\frac{1}{2} (-2+x) \HA_0(x)
\biggr]  G_{9}(x)
\Biggr\}
\nonumber\\&&
\end{eqnarray}

\begin{eqnarray}
F_2(y,\varepsilon) &=& 1
-2 y
\frac{1}{4} (-20+y) y \varepsilon
+\varepsilon ^2 \biggl\{
-\frac{1}{48} y \big(
480-765 y-56 y^2+64 y^3+12 y^4\big)
\nonumber\\&&
+\frac{1}{4} (-9+4 y) y^{2/3} (3+y)^{4/3}  G_{26}(y)
+(1-2 y)  G_{30}(y)
\biggr\}
+\varepsilon ^3 \Biggl\{
\frac{1}{192} y \big(
3840-21453 y
\nonumber\\&&
-1672 y^2+1638 y^3+280 y^4-6 y^5\big)
+ G_{26}(y) \biggl[
\frac{1}{16} \big(
243-108 y+2 y^2\big) y^{2/3} (3+y)^{4/3}
\nonumber\\&&
+\frac{1}{6} (9-4 y) y^{2/3} (3+y)^{4/3} \HA_0(y)
+\frac{1}{6} (-9+4 y) y^{2/3} (3+y)^{4/3} \HA_{-3}(y)
\nonumber\\&&
+\frac{1}{270} \big(
-1215+108 y+4 y^2\big) y^{2/3}  G_{24}(y)
+\frac{7}{5} (-1+2 y)  G_{25}(y)
+\frac{2}{3} (-1+2 y)  G_{28}(y)
\nonumber\\&&
-\frac{2}{3} (-1+2 y)  G_{29}(y)
\biggr]
-1215+108 y+4 y^2
\big)
\biggl[-\frac{1}{270} y^{2/3}  G_{27}(y)
-\frac{1}{270} y^{2/3}  G_{33}(y)
\biggr]
\nonumber\\&&
+(-1+2 y) \Bigl(
-\frac{7  G_{34}(y)}{5}
+\frac{2  G_{35}(y)}{3}
\Bigr)
\big(
\frac{2}{3}-\frac{4 y}{3}\big)  G_{31}(y)
+\frac{2}{3} (-1+2 y)  G_{32}(y)
\biggr]  G_{25}(y)
\nonumber\\&&
+\frac{1}{20} \big(
-52+204 y-5 y^2\big)  G_{30}(y)
+\frac{1}{6} (-9+4 y) y^{2/3}
(3+y)^{4/3}  G_{31}(y)
\nonumber\\&&
-\frac{1}{6} (3+y) (-9+4 y) y^{2/3} \sqrt[3]{3+y}  G_{32}(y)
-\frac{2}{3} (-1+2 y)  G_{36}(y)
\Biggr\}+\mathcal{O}(\varepsilon^4).
\end{eqnarray}

By multiplying $F_1(x)$ and $F_2(y)$ one obtains the series expansion of $\mathcal F(x,y)$, with $0<x<1, 0<y<1$.
The functions $G_i$ are iterated integrals [252] defined as in Eq. (<ref>), and they are listed in Appendix <ref>.

§.§ Example 3

In [201] we also considered as an example the system of two differential equations in two variables implied by the following differential operators,

\begin{eqnarray}
&&	-\alpha  \beta -x^6 \partial_x^6-\beta  x^5 \partial_x^5-5 x^5 y \partial_x^5 \partial_y-15 x^5 \partial_x^5-10 \beta  x^4 \partial_x^4-10 x^4 y^2 \partial_x^4 \partial_y^2-5 \beta  x^4 y \partial_x^4 \partial_y
\nonumber\\&&
-60 x^4 y \partial_x^4 \partial_y-66 x^4 \partial_x^4-26 \beta  x^3 \partial_x^3-10 x^3 y^3 \partial_x^3 \partial_y^3-10 \beta  x^3 y^2 \partial_x^3 \partial_y^2-90 x^3 y^2 \partial_x^3 \partial_y^2
\nonumber\\&&
-40 \beta x^3 y \partial_x^3 \partial_y-198 x^3 y \partial_x^3
\partial_y-96 x^3 \partial_x^3-18 \beta  x^2 \partial_x^2-5 x^2 y^4 \partial_x^2 \partial_y^4-10 \beta  x^2 y^3 \partial_x^2 \partial_y^3
\nonumber\\&&
-60 x^2 y^3 \partial_x^2 \partial_y^3-60 \beta  x^2 y^2 \partial_x^2 \partial_y^2-198 x^2 y^2
\partial_x^2 \partial_y^2-78 \beta  x^2 y \partial_x^2 \partial_y-192 x^2 y \partial_x^2 \partial_y
\nonumber\\&&
-38 x^2 \partial_x^2-\alpha  x \partial_x-2 \beta  x \partial_x+\gamma  \partial_x-x y^5 \partial_x \partial_y^5-5 \beta  x y^4 \partial_x
\partial_y^4-15 x y^4 \partial_x \partial_y^4
\nonumber\\&&
-40 \beta  x y^3 \partial_x \partial_y^3-66 x y^3 \partial_x \partial_y^3-78 \beta  x y^2 \partial_x \partial_y^2-96 x y^2 \partial_x \partial_y^2-36 \beta  x y \partial_{x y}-38 x y \partial_{x y}
\nonumber\\&&
\partial_{x y}+x \partial_x^2-2 x \partial_x-\beta  y^5 \partial_y^5-10 \beta  y^4 \partial_y^4-26 \beta  y^3 \partial_y^3-18 \beta  y^2 \partial_y^2-2 \beta  y \partial_y,
\\  [1cm]
-\alpha  \beta_1 -\beta_1  x^5 \partial_x^5-x^5 y \partial_x^5 \partial_y-10 \beta_1  x^4 \partial_x^4-5 x^4 y^2 \partial_x^4 \partial_y^2-5 \beta_1  x^4 y \partial_x^4
\partial_y-15 x^4 y \partial_x^4 \partial_y
\nonumber\\&&
-26 \beta_1  x^3 \partial_x^3-10 x^3 y^3 \partial_x^3 \partial_y^3-10 \beta_1  x^3 y^2 \partial_x^3 \partial_y^2-60 x^3 y^2 \partial_x^3 \partial_y^2-40 {\beta_1} x^3 y \partial_x^3 \partial_y
\nonumber\\&&
-66 x^3 y \partial_x^3 \partial_y-18 \beta_1  x^2 \partial_x^2-10 x^2 y^4 \partial_x^2 \partial_y^4-10 \beta_1  x^2 y^3 \partial_x^2 \partial_y^3-90 x^2 y^3 \partial_x^2
\partial_y^3
\nonumber\\&&
-60 \beta_1  x^2 y^2 \partial_x^2 \partial_y^2-198 x^2 y^2 \partial_x^2 \partial_y^2-78 \beta_1  x^2 y \partial_x^2 \partial_y-96 x^2 y \partial_x^2 \partial_y-2 \beta_1  x \partial_x
\nonumber\\&&
-5 x
y^5 \partial_x \partial_y^5-5 \beta_1  x y^4 \partial_x \partial_y^4-60 x y^4 \partial_x \partial_y^4-40 \beta_1  x y^3 \partial_x \partial_y^3-198 x y^3 \partial_x \partial_y^3
\nonumber\\&&
-78 \beta_1  x y^2 \partial_x
\partial_y^2-192 x y^2 \partial_x \partial_y^2-36 \beta_1  x y \partial_{x y}+x \partial_{x y}-38 x y \partial_{x y}-y^6 \partial_y^6-\beta_1  y^5 \partial_y^5
\nonumber\\&&
-15 y^5 \partial_y^5-10 \beta_1
y^4 \partial_y^4-66 y^4 \partial_y^4-26 \beta_1  y^3 \partial_y^3-96 y^3 \partial_y^3-18 \beta_1  y^2 \partial_y^2-38 y^2 \partial_y^2-\alpha  y \partial_y
\nonumber\\&&
-2 \beta_1  y \partial_y+\gamma
\partial_y+y \partial_y^2-2 y \partial_y.
\end{eqnarray}

Assuming a hypergeometric solution

\begin{equation}
\mathcal F(x,y) = \sum_{m,n\ge 0} f(m,n) x^m y^n,
\end{equation}

the coefficients $f(m,n)$ must obey

\begin{eqnarray}
(m+1) (\gamma +m+n) f(m+1,n)-(\beta +m) \left(\alpha +(m+n)^5+(m+n)^3\right) f(m,n)=0,\\
(n+1) (\gamma +m+n) f(m,n+1)-(\beta_1 +n) \left(\alpha +(m+n)^5+(m+n)^3\right) f(m,n)=0.
\end{eqnarray}

Solving these two equations with the help of Sigma, one obtains

\begin{eqnarray}
f(m,n) &=& \bigg(
\prod_{i_1=1}^n \frac{\big(
\big)
\big(-2
+8 i_1
-13 i_1^2
+11 i_1^3
-5 i_1^4
\big)}{\big(
\big) i_1}\bigg)
\nonumber\\&& \times
\prod_{i_1=1}^m \frac{\big(
\big)
\big) i_1}
\big(-2
+8 n
-13 n^2
+11 n^3
-5 n^4
+8 i_1
-26 n i_1
\nonumber\\&&
+33 n^2 i_1
-20 n^3 i_1
+5 n^4 i_1
-13 i_1^2
+33 n i_1^2
-30 n^2 i_1^2
+10 n^3 i_1^2
+11 i_1^3
-20 n i_1^3
\nonumber\\&&
+10 n^2 i_1^3
-5 i_1^4
+5 n i_1^4
\big) ~.
\end{eqnarray}

This quantity cannot be analytically expressed as a product of Pochhammer symbols due to the high degree of the polynomials appearing.

§.§ A brief descriptions of the commands of HypSeries

Techniques for the solution and classification of the hypergeometric systems (<ref>) are implemented in the Mathematica package HypSeries, which is attached to the paper [201].
The package requires Sigma, EvaluateMultiSums and HarmonicSums [202, 203, 204, 260, 282] to be loaded.

The commands

\[
{\tt solveDE1},~~
{\tt solveDE2},~~
{\tt solveDE3},~~
{\tt solveDE4}
\]

check whether the corresponding set of one to four equations in as many variables has hypergeometric solutions by consulting internal lists of cases, i.e. those discussed in Sec. <ref>. The syntax is e.g.

\[{\tt {solveDE4}[\{eq1==0,eq2==0,eq3==0,eq4==0\},\{x,y,z,t\},\{m,n,p,q\}]}.\]

More general solutions are possible by using the command DEProductSolution. One has to provide the required $n$ differential equations in the list

\[ {\tt sys = \{eq1==0,...,eqn==0\}}.\]

\[
{\tt DEProductSolution[sys,\{x,y,...\},\{m,n,...\}]}
\]

returns the respective expansion coefficient f[m,n,p,q]. Here the tools described in Section <ref> are utilized.

Conversely, from a Pochhammer ratio A = f[m,n,p,q] the command

\[
{\tt {findDE}[A,\{x,y,...\},\{m,n,...\}]}
\]

returns the system of differential equations obeyed by

\[
f(x,y,...) = \sum_{m,n,...\ge =0}^\infty f[m,n,...] \, x^m y^n \cdots
\]

Given a differential equation eq in $n$ variables, the command findRE

\[
{\tt findRE[eq==0,\{x,y,\ldots\}, f[m,n,\ldots]]\} }
\]

returns a corresponding recurrence for f[m,n,...]. The last two commands implement the techniques presented in the beginning of Section <ref>.

The convergence conditions for a number of two- and three-variate hypergeometric functions can be accessed from internal tables using the commands findCond2 and
findCond3. One first has to determine the corresponding function label fcn
via classifier2, classifier3, as e.g.

\[
\tt classifier3[f[m,n,p], \{x,y,z\}, \{m,n,p\}]
\]

returning fcn. Then

\[
{\tt findCond3[fcn,\{x,y,z\}}]
\]

returns the convergence conditions, which are in some cases given in implicit form and are taken from the tabulation in [321].

The function CheckDE[sol,eq] provides a way to check if an expression satisfies a differential equation in a series expansion. It returns a sum which, in that case, should be of higher order in $\ep$ if evaluated.

In [201] a notebook is attached where the usage of the package is presented along with the examples from the previous Sections, together with Mathematica files containing the definitions of the hypergeometric series treated by the package, their convergence conditions, and the translation table to their differential systems.

§ PARTIAL DIFFERENCE EQUATIONS WITH RATIONAL COEFFICIENTS

We examine the problem of linear partial difference equations in several variables, with the solution space being that of rational functions, possibly containing harmonic sums or Pochhammer symbols in the numerator. The corresponding problem in one variable is widely studied, and algorithms exist to find not only rational solutions but also hypergeometric solutions in a wide class of cases [376, 377, 378, 379, 380, 381]; these algorithms are implemented in the package Sigma. It is well-known how univariate difference equations arise frequently in the calculation of Feynman diagrams, for example, in calculating the master integrals in many applications in QCD, one can reduce the problem to a difference equation, which, for OMEs, is in the Mellin variable $N$, see e.g. [206, 207, 160].

In the multivariate case, there are fewer known algorithms that deal with difference equations than in the univariate case. We implemented in a Mathematica package called SolvePartialLDE [201] the algorithms of [382, 383] and describe here how we complemented them with flexible heuristic methods that may be useful, potentially, in future applications to Feynman integrals. As we will see in greater detail in the following sections, the basic idea is to constrain the denominator of the solution or at least parts of the denominator, in a way that we will describe more precisely. Once the denominator is constrained, the problem of finding the numerator can be reduced to that of a linear system of equations, and our program can accept an ansatz provided by a user for what type of object may appear in the numerator, chosen in the space of harmonic sums and/or Pochhammer symbols. The package can also solve the difference equation in a series expansion in one parameter, which in applications would be the dimensional parameter $\ep$, and can factor out a hypergeometric factor chosen by the user as we will elucidate in the following sections. The exposition in this chapter follows [201].

§.§ Description of the basic problem

A partial linear difference equation (PLDE) is an equation for an unknown function $y(n_1,\ldots,n_r)\in \mathbb K(n_1,\ldots,n_r)$, here a rational function in $r$ variables. We define the shift operators
$N_{\mathbf s}$ with respect to the shift $\mathbf s=(s_1,\ldots,s_r)\in \mathbb Z^r$ as:

\begin{equation}
N_\mathbf s y(n_1,\ldots,n_R) = y(n_1+s_1,\ldots,n_r+s_r).
\end{equation}

A PLDE is then an equation of the type

\begin{equation}
\sum_{\mathbf s \in S} a_\mathbf s N_\mathbf s y = f
\label{eq:PLDE}
\end{equation}

with $a_{\mathbf s}$ and $f$ polynomials in $n_1,\ldots,n_r$ and $S$ is a finite subset of $\mathbb Z^r$ called the shift set or structure set. Because (<ref>) is linear, its general solution is the sum of a particular solution and of a linear combination of solutions of the homogeneous equation with $f=0$.

An example of the type of equation under consideration is:

\begin{equation}
-(1 + k + n^2) y(n, k) + (4 + k + 2 n + n^2) y(1 + n, 2 + k) = 0.
\end{equation}

It has the shift set $S=\{(0,0), (1,2)\}$ and its coefficients are

\begin{eqnarray}
a_{(0,0)} &=& -(1+k+n^2),\\
a_{(1,2)} &=& (4 + k + 2 n + n^2).
\end{eqnarray}

For the purpose of the algorithms dealt with in our package, one must distinguish between periodic and aperiodic polynomials [382, 383]. A polynomial $p$ is periodic if there exist infinitely many shifts, mapping $p$ into $p'$, such that $\text{gcd}(p,p')\neq 1$. A polynomial is aperiodic if it is not periodic.
For example, with respect to the variables $\{n,k\}$, the polynomial $(n + k + 2)$ is periodic and the polynomial $(n^2 + k + 6)$ is aperiodic. Any polynomial can be factorized into a periodic and an aperiodic part. Given a PLDE, there are algorithms which constrain the periodic and the aperiodic part of the denominator of the solution. As we will describe further, the aperiodic part can always be constrained, but this is not guaranteed for the periodic part. In the following we describe our implementation choices in our Mathematica package [201].

§.§ Denominator bounds

The algorithms in [382, 383] aim at formulating a denominator bound for the solution of (<ref>). A denominator bound is a polynomial $d$ such that for any solution $y=\frac{n}{p}$ of (<ref>) it must be $p|d$. As observed already in the univariate case [376, 377], obtaining a denominator bound is valuable for the following reason: a naive way to solve the PLDE is by formulating an ansatz for the solution, i.e. a rational function in the variables $n_1,\ldots,n_k$ with undetermined coefficients $c_\mathbf k$,

\begin{equation}
y(n_1,\ldots,n_r)=\frac{\sum\limits_{\mathbf k} c_\mathbf k \prod\limits_i n_i^{k_i} }{\sum\limits_{\mathbf k'} c_{\mathbf k'} \prod\limits_i n_i^{k'_i}}.
\label{eq:rational-ansatz}
\end{equation}

By plugging this ansatz in (<ref>) and clearing denominators one obtains a set of constraints on the $c_\mathbf k$ and $c_{\mathbf k'}$ by imposing the equality of every monomial in the variables $n_i$ on both sides of the
equation. However, these equations are, in general, non-linear.

However, if the denominator bound can be obtained, then only an ansatz for the numerator is required and the  equations for the unknown coefficients are linear. The ansatz for the numerator can then be made to involve harmonic sums [246, 247] and/or Pochhammer symbols.

If we write the solution to a partial linear difference equation as $y=\frac{n}{uv}$ with $u$ aperiodic and $v$ periodic, it is always possible to calculate a bound $d_a$ for the aperiodic part $u$ of the denominator. We refer to [382] for a description of how the aperiodic denominator bound is calculated.

For the periodic part $v$ it is not always possible to obtain a complete denominator bound for a PLDE. This is illustrated for example by the equation

\begin{equation}
y(n+1,k) - y(n,k+1) = 0,
\label{eq:noDenBound}
\end{equation}

which is satisfied by $\frac{1}{(n+k)^\alpha}$ for any $\alpha \in \mathbb N$. Clearly, no polynomial can be a denominator bound for Eq. (<ref>).
This example shows that it is impossible to formulate in all cases a complete denominator bound, because arbitrary powers of periodic factors can appear in the denominator of the solution of some equations. However, it is often possible to formulate a partial bound for $v$ and to deduce some properties of its unknown factors. (A partial bound is a bound for some, but not all, the periodic factors). Specifically, the algorithm in [383] works by successively examining the periodic factors of $a_\mathbf p$ with $\mathbf p$ a “corner point” of the shift set of the equation (see [383] for a definition). It guarantees that $v$ has the property

\begin{equation}
v \mid (d_p \cdot v_{\text{semi-known}} \cdot v_{\text{unknown}})
\end{equation}

where $d_p$ is an explicitly given polynomial obtained through the algorithm, $v_{\text{semi-known}}$ is a polynomial whose factors are to be taken from $\{N_\mathbf s (p^m) | \mathbf s \in \mathbb Z^r, m \in \mathbb Z\}$, where $p\in P$ and $P$ is a subset of the periodic factors of the coefficients of the corner points, which is identified by the algorithm; and $v_{\text{unknown}}$ is a polynomial such that $\text{spread}(v_{\text{unknown}}) = V$ and $V$ is a lattice in $\mathbb Z^r$ identified by the algorithm. The spread of a polynomial $u$ is defined by

\begin{equation}
\text{spread}(u) = \{\mathbf s\in \mathbb Z^r \mid \text{gcd}(u,N_\mathbf s u) \neq 1 \}.
\label{eq:spread}
\end{equation}

In our implementation, if $P\neq\varnothing$ or $V\neq\varnothing$, then corresponding pieces of information are printed in order to guide the user to the formulation of an ansatz for $v_{\text{semi-known}}$ and $v_{\text{unknown}}$. Once the user has decided on an ansatz $d_\text{user}$ for the missing factors in the denominator, it can be included in the search when looking for the numerator of the solution through the option InsertDenFactor $\rightarrow$ $d_{\text{user}}$ of our package.

§.§ Determination of the numerator

Once the denominator has been constrained, possibly including an ansatz chosen by the user, one may look for the numerator of the solution. It has been shown in [384] based on [385] that in general the problem is unsolvable, and no algorithm exists that can find all solutions. Still, one can find some solutions by taking a polynomial ansatz $\text{num}(c_i)$ and substituting the rational function

\begin{equation}\label{Equ:PLDEAnsatz}
y=\frac{\text{num}(c_i)}{d_a d_p d_{\text{user}}}
\end{equation}

into the equation (<ref>). After clearing denominators, one can formulate a linear system for the unknown coefficients $c_i$ such that the equation is satisfied.

In our experience, the determination of the $c_i$ is more computationally demanding than the determination of the denominator bound; for this reason we propose the following strategy, which is implemented in our package as released in [201]. When the PLDE does not contain any symbolic parameters (such as the dimensional regulator $\ep$, or ratios of invariants) other than the shift variables, one may obtain constraints on the undetermined $c_i$ simply by plugging, sufficiently many times, random numerical values for the shift variables. This allows to quickly obtain a linear system for the $c_i$ without resorting to more expensive symbolic comparisons of the coefficients of many monomials and to Gaussian elimination in a potentially large and highly redundant system.

If there are symbols present, instead, one may consider performing a first pass with the symbols replaced by random numbers, with the purpose of identifying and removing redundant constraints. Then, after removing the redundant equations for the $c_i$, the system can be solved in a stepwise manner, i.e. considering one at a time the constraints produced by one monomial, and plugging the result in the rest of the equation. This is what our package does when the function SolvePLDE is called.

It is certainly possible that the use of random numbers to generate constraints can cause the system to generate two (or more) equations for the $c_i$ which are not independent. The probability of such an occurrence can be made arbitrarily small by choosing a sufficiently large range over which the random numbers are chosen. In any event, the consequence of an unfortunate draw of random numbers can only cause the software to output more functions misidentified as solutions when in fact they are not; it cannot cause the software to miss any solutions. By explicitly checking the result, one can guard against this remote possibility, at the expense of additional computation time. The user should be aware of this aspect when calling the function SolvePLDE.

In the following we elaborate further enhancements in order to extend the solution space from the rational function case to more general classes of functions. Besides the examples below, further examples for each aspect can be found in the Mathematica notebooks attached to the paper [201].

§.§.§ Treatment of a hypergeometric prefactor

Given a PLDE (<ref>)

\begin{equation}
\sum_{\mathbf s \in S} a_\mathbf s N_\mathbf s y = f
\label{eq:PLDE1}
\end{equation}

one may want to identify solutions of the form

\begin{equation}
y' = r y
\label{eq:hypFactor}
\end{equation}

with $r=r(n_1,\ldots,n_r)$ a hypergeometric function of its arguments, i.e. a function such that the ratio

\begin{equation}
\frac{N_{\mathbf e_i} r}{r} = \frac{r(n_1,\ldots,n_i+1,\ldots n_r)}{r(n_1,\ldots,n_i,\ldots,n_r)}
\end{equation}

is for all $i$ a rational function of the variables $n_i$. Examples of
hypergeometric functions are Pochhammer symbols, factorials,
$\Gamma$-functions, binomial symbols, and obviously rational functions and polynomials.

The demand (<ref>) implies that $y$ must satisfy another PLDE,

\begin{equation}
\sum_{\mathbf s \in S} a'_\mathbf s N_\mathbf s y' = f',
\label{eq:PLDE2}
\end{equation}

The transformation from (<ref>) to (<ref>) is useful whenever it is possible to formulate an ansatz for $r$. Once some specific form can be postulated for $r$, the equation (<ref>) is obtained by substitution and by exploiting the hypergeometric property.
Consider for example the equation

\begin{eqnarray}
(1 + k) (\varepsilon + k) (1 + k + n^2) \, y(n, k) - 2 k (2 + k + n^2) \, y(n, 1 + k)
\nonumber\\ &&
+ (1 + k) (\varepsilon + k) (2 + k + 2 n + n^2) \, y(1 + n, k) = 0.
\label{eq:example_5}
\end{eqnarray}

We assume that its solution is

\begin{equation}
y(n,k) = (\varepsilon)_k \, y'(n,k)
\end{equation}

with $y'$ a rational function of $n$ and $k$ and $(\varepsilon)_k$ the Pochhammer symbol

\begin{equation}
(\varepsilon)_k = \varepsilon(\varepsilon+1)\cdots(\varepsilon+k-1).
\end{equation}

Then one derives a difference equation for $y'$, namely

\begin{eqnarray}
&& (\varepsilon + k) \big[(1 + k^2 + n^2 + k (2 + n^2)) y'(n, k) -
2 k (2 + k + n^2) y'(n, 1 + k)
\nonumber\\&&
+ (2 + k^2 + 2 n + n^2 + k (3 + 2 n + n^2)) y'(1 + n, k)\big] = 0 .
\end{eqnarray}

We can now solve the new equation, obtaining

\begin{equation}
y'(n,k)=\frac{k}{1 + k + n^2}.
\end{equation}

From this we conclude that the solution of (<ref>) is

\begin{equation}
y(n,k) = (\varepsilon)_k \, \frac{k}{1 + k + n^2} C,
\end{equation}

for some constant $C\in\mathbb K(\varepsilon)$.

§.§.§ Finding solutions in terms of nested sums

In physical applications, it is known that the solutions of difference equations occurring in actual problems contain (cyclotomic) harmonic sums [246, 247, 249] or their generalizations. The algorithm we described can be adapted to look for these objects in the numerator of the solution by modifying the ansatz (<ref>) to be formed by a linear combination, with unknown coefficients, of a polynomial expression in a finite list of harmonic sums having coefficients in $\mathbb K[n_1,\ldots,n_r]$. The list of harmonic sums must be shift-stable, i.e. a shift in any of the variables must not introduce new harmonic sums not already included in the list, and they also should be linearly independent. These properties can be guaranteed by the quasi-shuffle algebra that the harmonic sums satisfy or by difference ring methods [255, 386, 368].

After plugging the ansatz into the equation and clearing denominators, the nested sums at shifted arguments can be rewritten using relations of the type

\begin{equation}
\label{Equ:SumShiftRules}
S_1(n+i) = \frac{1}{n+i} + S_1(n+i-1)
\end{equation}

and similarly for all other nested sums, until only unshifted nested sums appear. Then one applies a coefficient comparison on the power products in the harmonic sums and the shift variables. One notices how the number of unknowns $c_i$ and of equations rises very fast with the degree of the polynomial chosen for the ansatz and with the number of harmonic sums under consideration; the homomorphic image techniques described in the beginning of Section <ref> are instrumental to perform these calculations in reasonable time.

This heuristic method provides in many cases the desired solution. For instance, consider the equation from [201]

\begin{eqnarray}
&&(-k-1) \left(k+n^2+2 n+1\right) f(n,k)
+k \left(k+n^2+2 n+2\right) f(n,k+1)
\nonumber\\&&
+2 (k+1) \left(k+n^2+4 n+4\right) f(n+1,k)
-2 k \left(k+n^2+4 n+5\right) f(n+1,k+1)
\nonumber\\&&
-(k+1) \left(k+n^2+6 n+9\right) f(n+2,k)
+k \left(k+n^2+6 n+10\right) f(n+2,k+1)=0,
\nonumber\\
\end{eqnarray}

Looking for solutions of the form described, with a numerator of degree up to 2, containing the harmonic sums $S_1(n),S_1(k),S_{2,1}(n)$ the algorithm finds the denominator

\begin{equation}
d_p=1 + k + 2 n + n^2
\end{equation}

and the corresponding numerators of the solutions of the homogeneous equation:

\begin{eqnarray}
1, k, k^2, n, k n, S_1(k), k S_1(k), n S_1(k), S_1(k)^2,k S_1(n), k S_{2,1}(n).
\end{eqnarray}

§.§.§ Matching the solution to initial values

Our package allows to look for solutions that conform to a given set of initial conditions. This general solution is found by building a linear combination with undetermined coefficients of the solutions of the homogeneous equation, plus a particular solution of the equation. Next, the initial values are plugged in, and a system of equations is obtained. In the case that the system contains symbolic parameters other than the shift variables, the undetermined coefficients to be searched are not just numbers. In that case, the coefficients of the linear combination are taken to be general rational functions in the parameters up to some chosen degree. The combination of the solutions will be of particular importance for the next subsection.

§.§.§ Finding the solution in a series expansion

In many applications it is desirable to obtain the Laurent series expansion of the solution of a difference equation. This may be easier to achieve than the derivation of a complete solution, because, at each order in the expansion, it is possible to derive a difference equation where the expansion parameter is absent, therefore the linear system to find the coefficients $c_i$ can potentially be solved much more quickly. The procedure, described in the following, generalizes the univariate case given in [387]. It assumes that the initial values of the solution in its $\ep$-expansion are known.

Consider the case where (<ref>) contains a parameter $\varepsilon$ in the coefficients:

\begin{equation}
\sum_{\mathbf s \in S} a_\mathbf s(n_i,\varepsilon) N_\mathbf s y(n_i) = f(n_i,\varepsilon) ,
\label{eq:PLDE-eps}
\end{equation}

where the coefficients $a_{\mathbf s}(n_i,\varepsilon)$ are polynomials in the shift variables and in the parameter $\varepsilon$. We search for a solution of (<ref>) which has, around $\varepsilon=0$, a Laurent expansion starting from the power $\varepsilon^{-\ell}$ of the parameter, with $\ell$ known,

\begin{equation}
y_\varepsilon(n_i) = \varepsilon^{-\ell} y_{-\ell}(n_i) + \cdots + y_0(n_i) + \varepsilon y_1(n_i) + \cdots + \varepsilon^c y_c(n_i) ,
\label{eq:y_eps}
\end{equation}

and we assume that the right-hand side of the equation can be expanded in a series in $\varepsilon$ as

\begin{equation}
f =  \varepsilon^{-\ell} f_{-\ell}(n_i) + \varepsilon^{-\ell+1} f_{-\ell+1} (n_i) + \cdots .
\label{eq:f_eps}
\end{equation}

We assume also that the $a_{\mathbf s}(n_i,\varepsilon=0)$ are not all zero, so that an overall power of $\varepsilon$, if present in the equation, has been factored out.
Then, one may proceed by inserting (<ref>) and (<ref>) into (<ref>) and doing a coefficient comparison of the $\varepsilon^{-\ell}$ terms, obtaining

\begin{equation}
\sum_{\mathbf s \in S} a_\mathbf s(n_i,\varepsilon=0) N_\mathbf s y_{-\ell}(n_i) = f_{-\ell}(n_i).
\label{eq:PLDE_eps-l}
\end{equation}

Equation (<ref>) is now free of $\varepsilon$, which facilitates the task of finding a solution and reduces the computational time required. If (<ref>) can be uniquely solved for $y_{-\ell}$ and the solution matched to initial values, one can move to the next higher power in $\varepsilon$ by plugging the solution into (<ref>). In this new equation one does a coefficient comparison of the next power in $\varepsilon$ and solves for $y_{\ell+1}$. The process is repeated as many times as needed until all the terms of interest in the Laurent expansion are obtained.

For instance, consider the equation

\begin{eqnarray}
&&\big[3 (k+1) (n+1)+4 (n+1)+1\big] (4 k n^2 \varepsilon ^3+5 n \varepsilon +6 \varepsilon ^2+1) f(n+1,k+1)
\nonumber\\&&
-(3 k n+4 n+1) \big[4 (k+1) (n+1)^2 \varepsilon ^3+5 (n+1)
\varepsilon +6 \varepsilon ^2+1\big] f(n,k) = 0 .
\label{eq:example-series}
\end{eqnarray}

Together with a list of 25 initial values, our procedure to compute the expansion encounters at order
$\varepsilon^{-2},\varepsilon^{-1},\varepsilon^0$ the equations

\begin{eqnarray}
(-3 k n-4 n-1) f(n,k)+(3 k n+3 k+7 n+8) f(n+1,k+1) &=& \tau ,
\end{eqnarray}

with $\tau = 0, 5, 0$, respectively,
which are free of $\varepsilon$. The series solution of (<ref>) is found to be

\begin{equation}
f(n,k) = \frac{1}{\varepsilon ^2 (3 k n+4 n+1)}+\frac{5 n}{\varepsilon  (3 k n+4 n+1)}+\frac{6}{3 k n+4 n+1} + \mathcal O(\varepsilon).
\end{equation}

§.§.§ A brief descriptions of the commands of solvePartialLDE

The Mathematica package SolvePartialLDE.m implements the
algorithms for solving partial linear difference equations. It requires Sigma and HarmonicSums to be loaded. Additionally, the software Singular [388] must be installed, and made available by the interface [389] to Mathematica. The installation path of Singular can be set using the command appropriate for the user's system, e.g.

< <Singular.m

SingularCommand = "(path to)/Singular-4.1.3-x86_64-Linux/bin/Singular"

The functions available are

* $\mathtt {spread[p, q, \{n, k, ...\} (, \{eps,...\} ) ]}$: this function calculates the spread, Eq. (<ref>), of the polynomials $p$ and $q$, in the variables $n,k,\ldots$. The symbols in the optional list are treated as an extension to the field over which the polynomials are defined. If the polynomials $p$ and $q$ contain symbolic parameters other than $n,k,\ldots$, such as for instance the dimensional regulator $\varepsilon$, they must be declared in the second list.

* $\mathtt {dispersion[p, q, \{n, k, ...\} (, \{eps,...\} ) ]}$: this function calculates the dispersion (it is the maximum of the spread) of the polynomials $p$ and $q$ in the variables $n,k,\ldots$. The second optional list has the same function as in the function spread.

* $\mathtt {SolvePLDE[eq==rhs,f[n,k,...] ,(options) ] }$. This command solves the linear partial difference equation. It has the following available options:

* UseObject $\rightarrow$ list of Harmonic sums and/or Pochhammer symbols

Allows to define a list of harmonic sums and Pochhammer symbols to be searched in the numerator of the solution.

* PLDEdegBound $\rightarrow$ d

Allows to choose the total degree $d$ of the ansatz for the numerator of the solution. Defaults to 0.

* InsertDenFactor $\rightarrow$ factors

In the case the periodic denominator bound was
complete, the user may force the search to include factors in the denominator.

* PLDESymbols $\rightarrow$ list

Any symbols appearing other than the shift variables must be declared in list.

* InitialValues $\rightarrow$ list

A list of initial values in the form $\{\{var_1\rightarrow val_1, var_2\rightarrow val_2,\ldots, initial value\},\ldots\}$

* SymbolDegree $\rightarrow$ d

When initial conditions are provided, a linear combination of the homogeneous solutions is built, having as coefficients rational functions in the symbols. This option sets the maximum total degree of the numerator and denominator of those rational functions.

* $\mathtt { SolveExpand[eq==rhs,f[n,k,...], PLDEExpandIn\rightarrow\{\varepsilon,\ell_{min},\ell_{max}\} , }$

$ \mathtt{ InitialValues\rightarrow\{\ldots\} ,(options) ] }$ : this command solves the PLDE in a series expansion in the parameter $\varepsilon$, as described in Section <ref>. The options are the same as for SolvePLDE.

* $\mathtt { expandHypergPref[eq==rhs, f[n,k,...], fac] }$. This command derives a new equation whose solution has the hypergeometric factor fac removed, as described in Section <ref>.

§ A NUMERICAL LIBRARY FOR DIS STRUCTURE FUNCTIONS

We present a Fortran library [390] of use for numerically calculating the structure functions $F_2(x,Q^2)$, $g_1(x,Q^2)$ for electromagnetic current exchange and $F_3^{W^+ \pm W^-}(x,Q^2)$ for charged-current exchange in the asymptotic regime of large exchanged momentum $Q^2$, by convolution in $N$-space of parametrized parton distribution functions (PDFs) with Wilson coefficients. The PDFs can be evolved by the library from a parametrization in $N$-space at an initial scale $Q_0^2$ specifiable by the user, and are evolved in the fixed-flavour-number scheme by employing the splitting functions which are known up to $\mathcal O(a_s^3)$ [69, 70, 83]. The library accepts as inputs the combinations $u\pm\bar u$, $d\pm\bar d$ and $s\pm\bar s$ in the unpolarized case and $\Delta u+\Delta\bar u$, $\Delta d+\Delta\bar d$, $\Delta s+\Delta\bar s$ in the polarized case. It also encodes the Wilson coefficients for the above-mentioned deep-inelastic processes and for the Drell-Yan process and Higgs production, which have been calculated in [403]. In the case of $F_2(x,Q^2)$, the known two-mass contributions of the Wilson coefficients are included up to $\mathcal O(a_s^3)$ [121], whereas for $g_1(x,Q^2)$ and $F_3^{W^+ \pm W^-}(x,Q^2)$ single-mass contributions are considered [184, 199, 185, 186, 187, 290].

The applicability of these asymptotic representations is due to the factorization theorems for QCD [59, 112], and any intervening corrections due to heavy quark effects to the Wilson coefficients are suppressed as $\mathcal O(m_{c,b}^2/Q^2)$ in the region $Q^2\gg m_{b,c}^2$. Furthermore, the asymptotic values of the structure functions are used as an ingredient in many definitions of variable flavour number schemes [112, 118].

The $N$-space evaluation is obtained by the asymptotic expansion of the harmonic sums in the limit $|N|\to\infty$, see [247, 256, 257], together with recursion relations. In this way, it is possible to express the physical quantities in the complex $N$ plane. We show what accuracy is attainable in the evaluation of moments of the Wilson coefficients and perform the evolution of a set of test PDFs with our library.

§.§ The structure functions $\boldmath{F_2}$ and $\boldmath{F_L}$

In the fixed flavour number scheme with $N_F=3$, the structure functions $F_{2,L}(x,Q^2)$ are written as the sum of purely massless and massive contributions as follows [121]:

\begin{eqnarray}
F_i(x,Q^2) = F_i^{\rm massless}(x,Q^2) + F_i^{\rm heavy}(x,Q^2)~, \quad i=2,L.
\label{eq:Fi}
\end{eqnarray}

The massless part can be written as

\begin{eqnarray}
\frac{1}{x} F_i^{\rm massless}(x,Q^2)&=& \sum_{q} e_q^2 \Biggl\{\frac{1}{N_F} \Biggl[
\Sigma(x,\mu^2) \otimes C_{i,Q}^{\rm S}\left(x,\frac{Q^2}{\mu^2}\right)
+G\left(x,\mu^2\right) \otimes C_{i,g}\left(x,\frac{Q^2}{\mu^2}\right)\Biggr]
\NN\\&&
+ \Delta_{q}(x,\mu^2) \otimes C_{i,q}^{\rm NS}\left(x,\frac{Q^2}{\mu^2}\right) \Biggr\}~,
\quad i=2,L~,
\label{eq:F2L_Factorization}
\end{eqnarray}

with $\Sigma$ and $\Delta_k$ the flavor singlet and non-singlet distributions given by Eqs. (<ref>) and (<ref>), and $G$ denoting the gluon density.

In our code, for $F_2$, we consider the contributions of both the $c$ and the $b$ quark in the asymptotic region $Q^2\gg m_{c,b}^2$, limited to the OMEs which are fully known in $N$-space, as specified below.
The heavy quark part is then given by

\begin{eqnarray}
\label{eq:F2heavy}
\frac{1}{x}   F_{(2,L)}^{\rm  heavy}(x,N_F\!\!\!&+&\!\!\!2,Q^2,m_1^2,m_2^2) =\NN\\
L_{q,(2,L)}^{\sf NS}\left(x,N_F+2,\frac{Q^2}{\mu^2}
\right)
\otimes
\Bigl[f_k(x,\mu^2,N_F)+f_{\overline{k}}(x,\mu^2,N_F)\Bigr]
\NN\\ &&\hspace{14mm}
+\frac{1}{N_F}L_{q,(2,L)}^{\sf PS}\left(x,N_F+2,\frac{Q^2}{\mu^2}
\right)
\otimes
\Sigma(x,\mu^2,N_F)
\NN\\ &&\hspace{14mm}
+\frac{1}{N_F}L_{g,(2,L)}^{\sf S}\left(x,N_F+2,\frac{Q^2}{\mu^2}
\right)
\otimes
\Biggr\}
\NN\\ && \hspace{7mm}
+\tilde{\tilde{H}}_{q,(2,L)}^{\sf PS}\left(x,N_F+2,\frac{Q^2}{\mu^2}
\right)
\otimes
\Sigma(x,\mu^2,N_F)
\NN\\ &&\hspace{7mm}
+\tilde{\tilde{H}}_{g,(2,L)}^{\sf S}\left(x,N_F+2,\frac{Q^2}{\mu^2}
\right)
\otimes
\end{eqnarray}

The massive Wilson coefficients in their asymptotic form are given in [121] and read:

\begin{eqnarray}
\label{eqWIL12M}
L_{q,(2,L)}^{\sf NS}(N_F+2) &=&
a_s^2 \Bigl[A_{qq,Q}^{(2), {\sf NS}}(N_F+2)~\delta_2 +
\hat{C}^{(2), {\sf NS}}_{q,(2,L)}(N_F)\Bigr]
\NN\\
a_s^3 \Bigl[A_{qq,Q}^{(3), {\sf NS}}(N_F+2)~\delta_2
+  A_{qq,Q}^{(2), {\sf NS}}(N_F+2) C_{q,(2,L)}^{(1), {\sf NS}}(N_F+2)
\NN \\
&& \hspace*{5mm}
+ \hat{C}^{(3), {\sf NS}}_{q,(2,L)}(N_F)\Bigr]~,  \\
\label{eqWIL22M}
L_{q,(2,L)}^{\sf PS}(N_F+2) &=&
a_s^3 \Bigl[~A_{qq,Q}^{(3), {\sf PS}}(N_F+2)~\delta_2
+  A_{gq,Q}^{(2)}(N_F+2) N_F\Ctil_{g,(2,L)}^{(1)}(N_F+2) \NN \\
&& \hspace*{5mm}
+ N_F \hat{\Ctil}^{(3), {\sf PS}}_{q,(2,L)}(N_F)\Bigr]~,
\\
\label{eqWIL32M}
L_{g,(2,L)}^{\sf S}(N_F+2) &=&
a_s^2 A_{gg,Q}^{(1)}(N_F+2)N_F \Ctil_{g,(2,L)}^{(1)}(N_F+2)
\NN\\ &+&
a_s^3 \Bigl[~A_{qg,Q}^{(3)}(N_F+2)~\delta_2
+  A_{gg,Q}^{(1)}(N_F+2) N_F\Ctil_{g,(2,L)}^{(2)}(N_F+2)
\NN\\ && \hspace*{5mm}
+  A_{gg,Q}^{(2)}(N_F+2) N_F\Ctil_{g,(2,L)}^{(1)}(N_F+2)
\NN\\ && \hspace*{5mm}
+  ~A^{(1)}_{Qg}(N_F+2) N_F\Ctil_{q,(2,L)}^{(2), {\sf PS}}(N_F+2)
+ N_F \hat{\Ctil}^{(3)}_{g,(2,L)}(N_F)\Bigr]~,
\\ \NN \\
\label{eq:WILPS2M}
\tilde{\tilde{H}}_{q,(2,L)}^{\sf PS}(N_F+2)
&=& \sum_{i=1}^2 e_{Q_i}^2 a_s^2 \Bigl[~A_{Qq}^{(2), {\sf PS}}(N_F+2,m_i^2)~\delta_2
+~\Ctil_{q,(2,L)}^{(2), {\sf PS}}(N_F+2)\Bigr]
\\
&+& a_s^3 \Bigl[~\tilde{\tilde{A}}_{Qq}^{(3), {\sf PS}}(N_F+2)~\delta_2
+ \sum_{i=1}^2 e_{Q_i}^2 \Bigl[~\Ctil_{q,(2,L)}^{(3), {\sf PS}}(N_F+2) \NN\\
&& %\hspace*{-20mm}
+ A_{gq,Q}^{(2)}(N_F+2)~\Ctil_{g,(2,L)}^{(1)}(N_F+2)
\NN\\&&
+ A_{Qq}^{(2), {\sf PS}}(N_F+2)~C_{q,(2,L)}^{(1), {\sf NS}}(N_F+2)
\Bigr] \Bigr]~,       \label{eqWIL42M}
\NN\\
\label{eq:WILS2M}
\tilde{\tilde{H}}_{g,(2,L)}^{\sf S}(N_F+2) &=&   \sum_{i=1}^2 e_{Q_i}^2\Bigl[
a_s \Bigl[~A_{Qg}^{(1)}(N_F+2)~\delta_2
+~\Ctil^{(1)}_{g,(2,L)}(N_F+2) \Bigr] \NN\\
&+& a_s^2 \Bigl[~A_{Qg}^{(2)}(N_F+2)~\delta_2
+~A_{Qg}^{(1)}(N_F+2)~C^{(1), {\sf NS}}_{q,(2,L)}(N_F+2)\NN\\ &&
\hspace*{5mm}
+~\Ctil^{(2)}_{g,(2,L)}(N_F+2) \Bigr]\Bigr]
\NN\\ &+&
a_s^3 \Bigl[~\tilde{\tilde{A}}_{Qg}^{(3)}(N_F+2)~\delta_2
+ \sum_{i=1}^2 e_{Q_i}^2 \Bigl[A_{Qg}^{(2)}(N_F+2)~C^{(1), {\sf NS}}_{q,(2,L)}(N_F+2)
\NN\\ &&
\hspace*{5mm}
\NN\\ && \hspace*{5mm}
C^{(2), {\sf NS}}_{q,(2,L)}(N_F+2)
+~\Ctil^{(2), {\sf PS}}_{q,(2,L)}(N_F+2)\Bigr\}
\NN\\ && \hspace*{5mm}
+~A_{gg,Q}^{(1)}(N_F+2)~\Ctil^{(2)}_{g,(2,L)}(N_F       +2)
+~\Ctil^{(3)}_{g,(2,L)}(N_F+2) \Bigr]\Bigr]~.
\label{eqWIL52M}
\end{eqnarray}

Here the symbol $\delta_2$ takes the values

\begin{eqnarray}
\delta_2 = \left\{ \begin{array}{lll} 1 & \text{for} & F_2 \\
0 & \text{for} & F_L.
\end{array} \right.
\end{eqnarray}

In this notation, we use

\begin{eqnarray}
\tilde{f}(x)&=&\frac{f(x)}{x}~, \label{tildenotation}
\\
\hat{f}(x)&=&f(x+2)-f(x)~. \label{hatnotation}
\end{eqnarray}

The OMEs $A_{ij}(N_F+2)$ have the structure

\begin{equation}
A_{ij}(N_F+2) = A_{ij}(m_1) + A_{ij}(m_2) + \bar A_{ij}(m_1,m_2) + \bar A_{ij}(m_2,m_1)
\end{equation}

in which $\bar A_{ij}(m_1,m_2)$ denotes the part for which the current couples to the quark of mass $m_1$. The double tilde in $\tilde{\tilde{H}}_{q,(2,L)}$ and $\tilde{\tilde{H}}_{g,(2,L)}^{\sf S}$ refers to the charge weighting in the third-order OMEs:

\begin{equation}
\tilde{\tilde{A}}_{ij}^{(3)} = e_{Q_1}^2 A_{ij}^{(3)}(m_1) + e_{Q_2}^2 A_{ij}^{(3)}(m_2)
+   e_{Q_1}^2 \bar{A}_{ij}^{(3)}(m_1,m_2)
+   e_{Q_2}^2 \bar{A}_{ij}^{(3)}(m_2,m_1)~.
\end{equation}

We implement in a numerical program the evolution of the structure function $F_2(x,Q^2)$ according to Eq. (<ref>), (<ref>), (<ref>), by performing in $N$-space the evolution of the PDFs and pairing them to the Wilson coefficients by multiplication in $N$-space. As a last step, by a numerical contour integral, $F_2(x,Q^2)$ is obtained.

At present, some of the quantities in Eqs. (<ref>)-(<ref>) are unknown in full analytic form. They are: $a_{Qg}^{(3)}$ and the two-mass asymmetrical OMEs $\bar A_{Qg}^{(3)}$ and $\bar A_{Qq}^{(3),{\sf PS}}$.
For this reason, $\bar A_{Qq}^{(3),{\sf PS}}$ and $\tilde{\tilde A}_{Qg}^{(3)}$ are not implemented in our code. For $a_{Qq}^{(3),{\sf PS}}$, which depends on generalized harmonic sums, we employ an approximate representation.

For $a_{Qg}^{(3)}$, approximate representations exist, based on partial information obtained from interpolations of fixed moments and from the known leading-logarithmic contributions to the heavy Wilson coefficients in the small-$x$ limit and their double-log contributions in the large-$x$ limit [183]. The approximation in [183] was further refined in [78]; it is not included at present in our code.

As for the massless Wilson coefficients, pertaining the approximate parametrizations of the three-loop Wilson coefficients, there appear to be discrepancies between Eq. (4.13) in the text of [71] and the attachments to the paper. In particular, an inverse Mellin transform of the formula presented in the attachments would suggest that a term in the 4th line in the text needs to be corrected to read $+932.089L_0$ i.e. a sign flip. A second inconsistency is in the attachment to the paper where the $N$-space approximation is given. There, it appears that the flavour constant $fl_{11}^g$ should be multiplied by $N_F^2$ instead of $N_F$. We suggest these changes so that the formulas become the Mellin transforms of each other and agree numerically with the moments presented in [67]: with the definitions

\begin{equation}
x_1=1-x, \qquad L_0 = \ln(x), \qquad L_1 = \ln(x_1), \qquad fl_{11}^g = \frac{\langle e\rangle^2}{\langle e^2\rangle}
\end{equation}

where $\langle e^i\rangle$ refers to the average of the $i$-th power of the  quark charges, the correct formula appears to be:

\begin{eqnarray}
c_{2,g}^{(3)}(x) &\cong& \textcolor{blue}{N_F} \Big\{\frac{966}{81} L_1^5 - \frac{1871}{18} L_1^4
+ 89.31 L_1^3 + 979.2 L_1^2 - 2405 L_1 +
1372 x_1 L_1^4
\nonumber\\&&
- 15729 - 310510 x + 331570 x^2 - 244150 x L_0^2 - 253.3 x L_0^5
\nonumber\\&&
+ L_0 L_1 (138230 - 237010 L_0) - 11860 L_0 -
700.8 L_0^2 - 1440 L_0^3
\nonumber\\&&
+ \frac{4961}{162} L_0^4 - \frac{134}{9} L_0^5 -
x^{-1} (6362.54 \textcolor{red}{+} 932.089 L_0) + 0.625 \delta(x_1) \Big\}
\nonumber\\&&
\textcolor{blue}{N_F^2} \Big\{\frac{131}{81} L_1^4 - 14.72 L_1^3 + 3.607 L_1^2 - 226.1 L_1 + 4.762 -
190 x - 818.4 x^2
\nonumber\\&&
- 4019 x L_0^2 - L_0 L_1 (791.5 + 4646 L_0) +
739.0 L_0 + 418.0 L_0^2 + 104.3 L_0^3
\nonumber\\&&
+ \frac{809}{81} L_0^4 + \frac{12}{9} L_0^5 +
84.423 x^{-1} \Big\}
\nonumber\\&&
\textcolor{blue}{fl_{11}^g N_F^2 } \Big\{3.211 L_1^2 + 19.04 x L_1 + 0.623 x_1 L_1^3 - 64.47 x +
121.6 x^2 - 45.82 x^3
\nonumber\\&&
- x L_0 L_1 (31.68 + 37.24 L_0) +
11.27 x^2 L_0^3 - 82.40 x L_0 - 16.08 x L_0^2
\nonumber\\&&
+ \frac{520}{81} x L_0^3 +
\frac{20}{27} x L_0^4 \Big\} ~.
\label{eq:approxC2g3}
\end{eqnarray}

The difference to the text of [71] is marked in red. The Mellin transform of (<ref>) is

\begin{eqnarray}
c_{2,g}^{(3)}(N) &\cong&
\textcolor{blue}{N_F} \bigg\{
\frac{1}{27000 (N-1)^2 N^6 (1+N)^6 (2+N)} \big(
96480000+473848000 N+1224036000 N^2
\nonumber\\&&
N^3+1862075200 N^4-699200400 N^5-39517220484 N^6
\nonumber\\&&
N^7+37122825948 N^8+130572396885 N^9+45633145850 N^{10}
\nonumber\\&&
N^{11}-42336233496 N^{12}-9514319157 N^{13}-27750330 N^{14}
\nonumber\\&&
\frac{1}{N^3 (1+N)^4}(474020+2034310 N+3399445 N^2+2702152
\nonumber\\&&
+1041370 N^4+147850 N^5+2405 N^6)
-\frac{3 \big(
8931+566662 N+8931 N^2\big) S_2}{100 N (1+N)^2}
\nonumber\\&&
-\frac{1610 S_2^2}{9 N}
-\frac{4 (-22825+1871 N) S_3}{9 N (1+N)}
-\frac{3220 S_4}{9 N}
\Big] S_1
\frac{48}{5 N (1+N)^3} \big(
102-1409 N
\nonumber\\&&
+306 N^2+102 N^3\big)
+\frac{(22825-1871 N) S_2}{3 N (1+N)}
-\frac{6440 S_3}{27 N}
\Big] S_1^2
\nonumber\\&&
\frac{-8931-566662 N-8931 N^2}{100 N (1+N)^2}
-\frac{3220 S_2}{27 N}
\Big) S_1^3
+\frac{(22825-1871 N) S_1^4}{18 N (1+N)}
-\frac{322 S_1^5}{27 N}
\nonumber\\&&
\frac{2 \big(
1185050+3903173 N+4558059 N^2+2229119
N^3+348023 N^4\big)}{5 N^2 (1+N)^3}
-\frac{6440 S_3}{27 N}
\Big] S_2
\nonumber\\&&
+\frac{(22825-1871 N) S_2^2}{6 N (1+N)}
\frac{\big(
23692069+46835338 N+23692069 N^2\big) S_3}{50 N
\nonumber\\&&
+\frac{(22825-1871 N) S_4}{3 N (1+N)}
-\frac{2576 S_5}{9 N}
-\frac{10 (47402+13823 N) \zeta_2}{N^2}
-\frac{474020 \zeta_3}{N}
\bigg\}
\nonumber\\&&
+\textcolor{blue}{N_F^2} \bigg\{
\frac{1}{27000 (N-1) N^6 (1+N)^3 (2+N)} \big(
8640000+8656000 N+10073200 N^2
\nonumber\\&&
+9115000 N^3-30552800
N^4+8345452 N^5+539404372 N^6-148164201 N^7
\nonumber\\&&
-258366915 N^8-67625199
N^9-24818805 N^{10}\big)
\frac{92920-7915 N+2261 N^2}{10 N^3}
\nonumber\\&&
+\frac{1104 S_2}{25 N}
+\frac{1048 S_3}{81 N}
\Big) S_1
\frac{3607}{1000 N}
+\frac{262 S_2}{27 N}
\Big) S_1^2
+\frac{368 S_1^3}{25 N}
+\frac{131 S_1^4}{81 N}
\nonumber\\&&
+\frac{(9292000-787893 N) S_2}{1000 N^2}
+\frac{131 S_2^2}{27 N}
+\frac{233036 S_3}{25 N}
+\frac{262 S_4}{27 N}
+\frac{(-18584+1583 N) \zeta_2}{2 N^2}
\nonumber\\&&
-\frac{9292 \zeta_3}{N}
\bigg\}
\nonumber\\&&
+\textcolor{blue}{fl_{11}^g} \textcolor{blue}{N_F^{\textcolor{red}{2}}} \bigg\{
\frac{1}{13500 (1+N)^5 (2+N)^4 (3+N)} \big(
97917774+305888240 N
\nonumber\\&&
+418317438 N^2
+369372616 N^3+267827242 N^4+162867261 N^5+71110703 N^6
\nonumber\\&&
+19153530 N^7
+2741175 N^8+152685 N^9\big)
\frac{13749-34880 N-9520 N^2}{500 (1+N)^3}
\nonumber\\&&
-\frac{1869 S_2}{1000 N (1+N)}
\Big) S_1
3211+8291 N+3211 N^2\big) S_1^2}{1000 N (1+N)^2}
-\frac{623 S_1^3}{1000 N (1+N)}
\nonumber\\&&
3211+51091 N-28469 N^2\big) S_2}{1000 N (1+N)^2}
+\frac{7 (-89+5320 N) S_3}{500 N (1+N)}
+\frac{2 (-535+396 N) \zeta_2}{25 (1+N)^2}
\nonumber\\&&
-\frac{1862 \zeta_3}{25 (1+N)}
\bigg\} .
\label{eq:eq:approxC2g3N}
\end{eqnarray}

The difference to the attachment of [71] is marked in red.

For the structure function $F_L$ for photon exchange, the charm-quark contributions to the Wilson coefficients have been calculated to $\mathcal O(a^2)$ in [159, 181] with full mass dependence and in the asymptotic approximation, and to $\mathcal O(a^3)$ in [165]. The massless Wilson coefficients are known from [64, 76, 71]. It is known that the asymptotic approximation to the charm contribution at $\mathcal O(a^2)$ differs significantly from the analytic expression for virtualities below $Q^2\sim 1000\text{~GeV}^2$. For this reason, in our numerical library, only the massless Wilson coefficient for $F_L$ are coded. At $\mathcal O(a^3)$, our program uses the approximate formulas of [76]. The structure function $F_L(x,Q^2)$ is coded according to Eq. (<ref>).

§.§ The structure function $g_1$

In the case of $g_1(x,Q^2)$, we consider the contributions to the Wilson coefficients due to one heavy quark. Their explicit form can be found in [184, 199], and reads:

\begin{eqnarray} \frac{1}{x}g_1(x,Q^2) &=&
\sum_{k=1}^{N_F}e_k^2\Biggl\{
L_{q,g_1}^{\sf NS}\left(x,N_F+1,\frac{Q^2}{\mu^2}
\otimes
\Bigl[\Delta f_k(x,\mu^2,N_F)+\Delta f_{\overline{k}}(x,\mu^2,N_F)\Bigr]
\N\\ &&\hspace{14mm}
+\frac{1}{N_F}L_{q,g_1}^{\sf PS}\left(x,N_F+1,\frac{Q^2}{\mu^2}
\otimes
\Delta \Sigma(x,\mu^2,N_F)
\N\\ &&\hspace{14mm}
+\frac{1}{N_F}L_{g,g_1}^{\sf S}\left(x,N_F+1,\frac{Q^2}{\mu^2}
\otimes
\Delta G(x,\mu^2,N_F)
\Biggr\}
\N\\
H_{q,g_1}^{\sf PS}\left(x,N_F+1,\frac{Q^2}{\mu^2}
\otimes
\Delta \Sigma(x,\mu^2,N_F)
\N\\ &&\hspace{7mm}
+H_{g,g_1}^{\sf S}\left(x,N_F+1,\frac{Q^2}{\mu^2}
\otimes
\Delta G(x,\mu^2,N_F)
\Biggr]~.
\end{eqnarray}

In our numerical program, the polarized splitting functions and Wilson coefficients are coded. They can be obtained from [391, 83] and from [85].

§.§ The structure function $xF_3^{W^+ - W^-}$

For charged-current DIS processes, the structure functions $F_i^{W^\pm}$ are defined by [394]

\begin{eqnarray}
\frac{d \sigma^{\nu(\bar{\nu})}}{dx dy}
&=& \frac{G_F^2 s}{4 \pi} \frac{M_W^4}{(M_W^2 + Q^2)^2}
\\ && \times
\Biggl\{
\left(1 + (1-y)^2\right) F_2^{W^\pm}(x,Q^2)
- y^2 F_L^{W^\pm}(x,Q^2)
\pm \left(1 - (1-y)^2\right) xF_3^{W^\pm}(x,Q^2)\Biggr\}
\nonumber\\
\frac{d \sigma^{l(\bar{l})}}{dx dy}
&=& \frac{G_F^2 s }{4 \pi} \frac{M_W^4}{(M_W^2 + Q^2)^2}
\\ && \times
\Biggl\{
\left(1 + (1-y)^2\right) F_2^{W^\mp}(x,Q^2)
- y^2 F_L^{W^\mp}(x,Q^2)
\pm \left(1 - (1-y)^2\right) xF_3^{W^\mp}(x,Q^2)\Biggr\}
\nonumber
\end{eqnarray}

in terms of the differential cross sections, with $y=q.P/l.P$, $x=Q^2/ys$, $s=(l+P)^2$, $G_F$ the Fermi constant, $M_W$ the $W$-boson mass and $l,P$ are the momenta of the incoming lepton and proton. The combinations

\begin{eqnarray}
xF_3^{W^+ \mp W^-}(x,Q^2) = xF_3^{W^+}(x,Q^2) \pm xF_3^{W^-}(x,Q^2)
\label{eq:F3pm}
\end{eqnarray}

are usually considered because they exist for odd or even moments respectively [40].

Here we consider the combination $xF_3^{W^+ - W^-}(x,Q^2)$, for which the asymptotic charm contributions to the Wilson coefficients have been calculated in [187] to $\mathcal O(a^2)$ correcting the results in [188], and in [290] to $\mathcal O(a^3)$, and the exact $\mathcal O(a)$ in [186].

One can write from the factorization theorems [290]

\begin{eqnarray}
F_3^{W^+ - W^-}(x,Q^2) &=&
\Bigl[
|V_{du}|^2 (d - \overline{d})
+ |V_{su}|^2 (s - \overline{s})
+ V_{u}      (u - \overline{u})\Biggr] \otimes \Biggl[
C_{q,3}^{W^+ - W^-, \rm NS} +
L_{q,3}^{W^+ - W^-, \rm NS} \Biggr] \nonumber\\
+ \Bigl[
|V_{dc}|^2 (d - \overline{d})
+ |V_{sc}|^2 (s - \overline{s}) \Biggr] \otimes
H_{q,3}^{W^+ - W^-, \rm NS}~.
\label{eq:F3WMW}
\end{eqnarray}

Notice how there is no dependence on the gluon density. The coefficients $V_{ij}$ are those of the CKM matrix, whose values are [395]

\begin{eqnarray}
|V_{du}| &=& 0.97370 \\
|V_{su}| &=& 0.2245  \\
|V_{dc}| &=& 0.221   \\
|V_{sc}| &=& 0.987
\end{eqnarray}

\begin{equation}
V_u = |V_{du}|^2 + |V_{su}|^2.
\end{equation}

The asymptotic expressions of the heavy quark Wilson coefficients take the form [290]

\begin{eqnarray}
L_{q,3}^{W^+ - W^-,\rm NS}(N_F+1) &=& a_s^2 A_{qq,Q}^{(2),\rm NS} + \hat{C}_{q,3}^{(2),W^+ - W^-, \rm NS}(N_F) +
a_s^3 \Biggl[A_{qq,Q}^{(3),\rm NS}
\nonumber\\  &&
+ A_{qq,Q}^{(2),\rm NS} C_{q,3}^{(1),W^+ - W^-, \rm NS}(N_F+1) + \hat{C}_{q,3}^{(3),W^+ - W^-, \rm NS}(N_F)\Biggr]~,
\label{eq:WIL1}
\\
H_{q,3}^{W^+ - W^-,\rm NS}(N_F+1) &=&
L_{q,3}^{W^+ - W^-,\rm NS}(N_F+1) + {C}_{q,3}^{W^+ - W^-, \rm
\label{eq:H3}
\end{eqnarray}

with the notation

\begin{equation}
\hat f(N_F) = f(N_F+1)-f(N_F).
\end{equation}

In our program, we calculate $xF_3^{W^+ - W^-}(x,Q^2)$ numerically in $N$-space using Eq. (<ref>) by pairing the valence PDFs obtained by evolving an input parametrization and the logarithmic contributions to the Wilson coefficients computed for $\mu^2=Q^2$ where $\mu^2$ denotes the factorization and renormalization scales.

§.§ The structure function $xF_3^{W^+ + W^-}$

For the structure function $xF_3^{W^+ + W^-}$, defined as in (<ref>), the asymptotic factorization derived in [187] reads:

\begin{align}
\label{eq:3flOddComb3}
F_3^{W^+ + W^-}
= &
\phantom{+}
\left(
\vert V_{du} \vert^2
(d + \bar{d})
\vert V_{su} \vert^2
(s + \bar{s})
(u + \bar{u})
\right)
\nonumber\\ &
\left(
\vert V_{dc} \vert^2
(d + \bar{d})
\vert V_{sc} \vert^2
(s + \bar{s})
\right)
\nonumber\\ &
\left[
\Sigma
\right]
\end{align}

Explicit formulas for the asymptotic Wilson coefficients are given in [187] in $N$-space, obtained from the factorization theorem. These formulas are encoded in our numerical program.

§.§ Drell-Yan process

The Drell-Yan process refers to the inclusive lepton pair production from two hadrons,

\begin{equation}
H_1 + H_2 \to \ell_1+\ell_2+X,
\end{equation}

with the invariant mass of the lepton pair denoted by $Q^2$ and the CM energy of the hadrons by $s$.
From the mass factorization theorems [57, 396, 58], the hadronic structure function can be written as

\begin{equation}
W^{DY}(x,Q^2) = \sum_{i,j=q,\bar q,g} \int_0^1 dx_1\int_0^1 dx_2 \int_0^1 dz \ \delta(x-x_1 x_2 z) \ f_{i/H_1}(x_1,\mu^2) \ f_{j/H_2}(x_2,\mu^2) \ \Delta_{ij}^{DY} \Bigl(z,\frac{Q^2}{\mu^2}\Bigr) ,
\label{eq:sigmaDY}
\end{equation}

where $f_{i,j}$ are the parton distribution functions, and $\mu^2$ is the mass factorization and renormalization scale, here set to be equal, and

\begin{equation}
x=\frac{Q^2}{s} .
\end{equation}

The Drell-Yan structure function is related to the differential cross-section by

\begin{equation}
\frac{1}{x}\frac{d\sigma^{DY}(x,Q^2)}{dQ^2} = \sigma^V(Q^2) \ W^{DY}(x,Q^2) ,
\end{equation}

and $\sigma^V(Q^2)$ is the point-like cross-section, which depends on which of the standard-model bosons $V$ is  exchanged. In the case of photon exchange, one has [397]

\begin{equation}
\sigma^\gamma(Q^2) = \frac{16\pi^2 a_s^2}{3 N_c Q^4} .
\end{equation}

The Wilson coefficients, $\Delta_{ij}^{DY}$, can be calculated in perturbation theory and have been computed at order $a_s$ in [398, 399, 400, 401] and at order $a_s^2$ in [397], see also [193]. Recently, the order $a_s^3$ has been discussed in [404]; the calculation limited to photon exchange only was presented in [405] and for charged-current exchange in [406].

After a Mellin transform, Eq. (<ref>) can be written as:

\begin{equation}
\mathbf M[W^{DY}](N,Q^2) = \mathbf M[f_i](N,\mu^2) \ \mathbf M[f_j](N,\mu^2) \ \mathbf M[\Delta_{ij}^{DY}]\Bigl(N,\frac{Q^2}{\mu^2}\Bigr) .
\end{equation}