# Simulating the Femtouniverse on a Quantum Computer

Nouman<EMAIL_ADDRESS>Patrick<EMAIL_ADDRESS>Jiayu
<EMAIL_ADDRESS>Illinois Quantum Information Science and
Technology Center, Urbana, IL 61801
Illinois Center for Advanced Studies of the Universe, Urbana, IL 61801
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL
61801

###### Abstract

We compute the low-lying spectrum of 4D SU(2) Yang-Mills in a finite volume
using quantum simulations. In contrast to small-volume lattice truncations of
the Hilbert space, we employ toroidal dimensional reduction to the
“femtouniverse” matrix quantum mechanics model. In this limit the theory is
equivalent to the quantum mechanics of three interacting particles moving
inside a 3-ball with certain boundary conditions. We use the variational
quantum eigensolver and quantum subspace expansion techniques to compute the
string tension to glueball mass ratio near the small/large-volume transition
point, finding qualitatively good agreement with large volume Euclidean
lattice simulations.

##  I Introduction

Quantum simulations of quantum field theories (QFTs) in the noisy
intermediate-scale quantum (NISQ) era Preskill2018quantumcomputingin  will be
limited by the gap between the complexity of QFTs and the available quantum
resources. Currently, it is difficult to build a comparison between known
Monte Carlo and quantum simulation results. In particular, gauge theories are
challenging to simulate on quantum devices because of the infinite tower of
states associated with each link on the lattice. To control the rapid growth
in degrees of freedom requires working with small lattices and putting rather
stringent truncations on the subspaces of the Hilbert space associated with
bosonic degrees of freedom. Subsequently one must take care to preserve the
Gauss law and confront various limitations on what can be computed with small
circuits. It is therefore interesting to explore a wide variety of
approximations, truncations, encodings, and simulation techniques to make the
most of available near-term systems. In the process we may hope to develop new
applications and tools that will help lower the threshold for quantum
supremacy in simulating high energy physics phenomena, while creating a large
set of benchmark computations for calibrating quantum simulations against
analytic and classical results.

In this paper we study hybrid quantum simulations of an approximation to 4D
gauge theory known as the femtouniverse Bjorken:1979hv . Complementary to the
small-volume lattice approach, the approximation we use is a matrix quantum
mechanics model obtained by the dimensional reduction of the 4D theory on a
spatial 3-torus. The study of gauge theory on a torus was pioneered by ’t
Hooft thooft  and the effective matrix quantum mechanics theory obtained by
dimensional reduction was systematically developed by Lüscher Luscher:1982ma ,
Lüscher and Münster LUSCHER1984445 , van Baal VANBAAL1986548 ; van1988gauge ;
vanBaal1990 , van Baal and Koller vBK ; van1987qcd , and others. For a
comprehensive review see van2001qcd . We focus on the SU(2) model studied in
detail by van Baal and Koller vBK ; van1987qcd . In our quantum simulations we
use the variational quantum eigensolver (VQE) peruzzo2014variational ;
tilly2022variational  and quantum subspace expansion (QSE) PhysRevA.95.042308
techniques. These are hybrid classical-quantum methods that exploit the
strengths of both types of computation. With these methods we obtain the low-
lying spectrum of the theory, including glueball masses and string tensions,
which can be compared with the large-volume continuum limit results known from
Euclidean lattice Monte Carlo Teper .

In small volumes, the effective theory from integrating out Kaluza-Klein modes
is accurate due to asymptotic freedom. However, it differs markedly from the
large-volume limit. At large torus volumes, the spectrum of the complete 4D
theory is quite insensitive to the volume. There is a “large-small volume
transition region”, roughly around $\Lambda L\sim\mathcal{O}(1)$, where some
strong-coupling physics of the large volume theory is reflected with
reasonable accuracy in the much simpler effective theory. This is the regime
we would like to probe.

Our work is organized as follows. In Section II we begin by reviewing the
effective Hamiltonian described in vBK  and the symmetries of the effective
theory, which allow the decomposition of the Hamiltonian into different
superselection sectors. In Section III we briefly review the QSE algorithm
used in the calculation of excited state energies. Section IV presents our
main results. We compute ground and excited state energies using classical
exact diagonalization (ED), VQE+QSE simulators, and the IBM-Lima quantum
computer. We discuss effects of Hilbert space truncation and different basis
choices and various sources of error. We find that the string tension to
glueball mass ratio of the continuum large-volume theory is reasonably well
captured by VQE+QSE on the real device, working at couplings near the
large/small volume transition.

The study of matrix quantum mechanics in quantum simulations is motivated both
as a simpler approximation to higher-dimensional nonabelian gauge theories,
and by quantum gravity, where matrix models arise in various contexts.
Interesting prior work on this subject includes gharibyan2021toward ;
PRXQuantum.3.010324 . In contrast to these works we use a minimal, gauge-
invariant Hilbert space and focus on the explicit model that arises from
dimensional reduction on a 3-torus, including the loop-induced effective
potential and gauge field topology associated with the parent 4D theory.
Dimensional reduction is also a useful approach for truncating abelian
theories while preserving some structure; previously we have studied the
reduction of the Schwinger model to the quantum mechanical particle on the
circle model in quantum simulation PhysRevD.105.074505 , where dynamics
associated with an ’t Hooft anomaly and the $\theta$ term Hooft1980  are
transparent.

##  II Matrix model for $4D$ $\mathrm{SU}(2)$ gauge theory

We begin by reviewing the physics of the femtouniverse. In this section, we
rely heavily on the pioneering work of van Baal and Koller van1987qcd ; vBK .
We consider pure SU(2) Yang-Mills theory in 4D on a small spatial 3-torus of
length $L$ and work in $A_{0}=0$ gauge. The Hamiltonian is

$H=\int_{[0,L]^{3}}d^{3}\mathbf{x}\bigg{(}\frac{1}{2}g^{2}E^{a}_{k}(\mathbf{x})E^{a}_{k}(\mathbf{x})+\frac{1}{2g^{2}}B^{a}_{k}(\mathbf{x})B^{a}_{k}(\mathbf{x})\bigg{)}$
(1)

where $E^{a}_{k}$, $B^{a}_{k}$ are chromo-electric and chromo-magnetic fields,
and $g$ is the dimensionless strong coupling. Here $k$ denotes the spatial
index and $a$ labels the color index. The gauge connection $A^{a}_{k}$ is
taken to satisfy periodic boundary conditions on the torus.444Working in
sectors of nonzero magnetic flux would be an interesting direction for
generalization.

From lattice analysis we know the theory is gapped near the dynamical scale
$\Lambda=\mu e^{-8\pi^{2}/(bg^{2}(\mu))}$ and becomes exponentially
insensitive to $L$ for $\Lambda L\gtrsim 1$. The effective field theory (EFT)
analysis below, on the other hand, will be valid for $\Lambda L\lesssim 1$,
where $g$ may be thought of as the running coupling at the renormalization
group scale $1/L$. Thus we will be particularly interested in the behavior
near the small volume-large volume transition.

The proper gauge transformations of the SU(2) theory are periodic SU(2)-valued
functions $g(\mathbf{x})$ acting on the connections as

$\displaystyle
A_{k}(\mathbf{x})=g(\mathbf{x})A_{k}(\mathbf{x})g^{-1}(\mathbf{x})-ig(\mathbf{x})\partial_{k}g^{-1}(\mathbf{x}).$
(2)

The classical vacuum manifold, sometimes called the “vacuum valley”, is the
space of flat connections modulo small gauge transformations. The vacua
separate into a union of disjoint sectors characterized by different Chern-
Simons numbers and related by large gauge transformations. We will work in the
sector of fixed vanishing Chern-Simons number, which is acceptable as long as
the states of interest have energies below the “sphaleron” energy barrier
between these sectors.

We may partially fix the gauge so that the classical vacua are given by
spatially constant connections that are aligned in $su(2)$ space. For
convenience of illustration we can go to a gauge where the vacua are of the
Abelian form

$\displaystyle A_{i}=\frac{C_{i}}{L}\frac{\sigma_{3}}{2}.$ (3)

At this point, it appears that the vacuum valley (in a fixed Chern-Simons
sector) is $\mathbb{R}^{3}$. However, there are still residual gauge
transformations, given by

$\displaystyle g(\mathbf{x})$ $\displaystyle=\exp(-4\pi
i\frac{\mathbf{x}\cdot\mathbf{n}}{L}\frac{\sigma_{3}}{2}),$ $\displaystyle g$
$\displaystyle=\sigma_{1}\;.$ (4)

Here $n_{i}$ are integers, and $g=\sigma_{1}$ is the nontrivial element of the
Weyl group. These gauge transformations lead to the following identifications:

$\displaystyle\mathbf{C}\sim\mathbf{C}+4\pi\mathbf{n}$ (5)
$\displaystyle\mathbf{C}\sim-\mathbf{C}.$ (6)

Thus the vacuum valley is the orbifold $T^{3}/\mathbb{Z}_{2}$, where the torus
periodicity is $4\pi$.

At the quantum level, the continuous vacuum valley is lifted, but a discrete
global symmetry is preserved. This is the $\mathbb{Z}_{2}$ electric center
symmetry, which, in the gauge Eq. (3), acts as

$\displaystyle h(\mathbf{x})=\exp(-2\pi
i\frac{\mathbf{x}\cdot\mathbf{n}}{L}\frac{\sigma_{3}}{2})$ (7)

where $n_{i}\in\\{0,1\\}$. These are not gauge transformations because
$h(\mathbf{x})$ and $h(\mathbf{x}+L\mathbf{u}_{i})$, where $\mathbf{u}_{i}$ is
the $i$-th spatial unit vector, differ by a nontrivial element of the
$\mathbb{Z}_{2}$ center of $\mathrm{SU}(2)$. They are global symmetries,
however, because they preserve the action and the periodic boundary conditions
for the gauge field. They act on the classical vacua, which are given by
$\mathbf{C}=2\pi\mathbf{n}$, as

$\displaystyle\mathbf{C}\sim\mathbf{C}+2\pi\mathbf{n}.$ (8)

Therefore we expect 8 minima of the quantum induced potential on the vacuum
valley, lying on the corners of a cube embedded in the 3-torus. We also expect
that the symmetry will be unbroken and the true ground state will be similar
to a symmetric linear combination of the perturbative ground states around
each of these minima. In general the states can be taken to transform in
representations

$\displaystyle\ket{\psi(A^{h})}=(-1)^{\mathbf{k}\cdot\mathbf{e}}\ket{\psi(A)}$
(9)

where $\mathbf{e}$ is $\mathbb{Z}_{2}$-valued electric flux labeling the
representations under center thooft .

It is convenient to work with an effective theory that is partway between the
full theory Eq. (1) and the vacuum valley theory of the $C_{i}$ alone. We
split the gauge field into a linear combination of a spatially constant part
and a spatially varying field:

$\displaystyle A^{a}_{k}(\mathbf{x},t)=c^{a}_{k}(t)+Q^{a}_{k}(\mathbf{x},t).$
(10)

At small $L$ we can integrate out $Q$ to obtain an effective quantum
mechanical theory of the $c^{a}_{k}$. This theory includes, in addition to the
vacuum valley degrees of freedom, the “nonabelian” or “transverse” spatially
constant modes. Although it seems natural to integrate out Kaluza-Klein modes,
it is not obvious that this step is consistent in a gauge theory, where the
energy of a mode depends on the vacuum one is perturbing around. The validity
of this step hinges on being able to stay relatively close to the origin in
field space, with boundaries and suitable boundary conditions at points where
the constant field effective theory breaks down. We discuss this further
below.

It is also convenient to relax our identification of the vacuum valley with
the $\sigma^{3}$ direction. By a constant gauge transformation we can “point”
the vacuum valley in any direction. We will leave this gauge degree of freedom
in the effective theory, then remove it later by averaging over all
directions. In this formulation, the vacuum valley corresponds to $c_{k}$ that
are aligned in $su(2)$ space with any common direction. It may be visualized
as three aligned (or antialigned) vectors in three dimensions, constrained by
the periodicities Eq. (6) to have magnitudes less than $2\pi/L$. The relative
angles between the vectors then encode the transverse degrees of freedom.

The effective theory described thus far is the quantum mechanics of three
interacting particles in three dimensions. In fact, we can also constrain the
particles to move inside a ball,

$\displaystyle r_{i}\equiv\sqrt{\sum_{a}c_{i}^{a}c_{i}^{a}}\leq\frac{\pi}{L},$
(11)

with certain boundary conditions at $\pi/L$:

$\displaystyle\left.\left(\frac{\partial}{\partial
r_{i}}\right)^{1-e_{i}}(r_{i}\psi(c))\right|_{r_{i}=\frac{\pi}{L}}=0.$ (12)

The arguments go as follows vBK . The points $r_{k}\in\\{0,\pi\\}$ in the
vacuum valley are preserved by the center and Weyl transformations. This is
clear in the gauge-fixed version of Eqs. (6), (8) and Eq. (11) generalizes it
to the case where the vacuum valley direction is unfixed. If the wavefunction
and the energy is continuous at these points, then either the wavefunction or
its (covariant) radial derivative must vanish, depending on the electric flux
quantum numbers (cf. Eq. (9)).

Continuity, however, is rather subtle. A global center transformation acts as
a $2\pi$ translation on the vacuum valley, but makes a nontrivial modification
to the spatial dependence of the transverse modes. Since it is a symmetry,
these modes must be reordered. The transverse Kaluza-Klein modes that should
be kept in a consistent effective theory are completely different from the
spatially constant modes valid near the origin, once we have translated by
$2\pi$ away from the origin along the vacuum valley. We should think of the
points $r_{k}\in\\{0,\pi/L\\}$ as boundaries between different patches in
which we have different EFTs. More precisely the “three particles in a ball”
description is valid in any patch, but the map to the underlying degrees of
freedom of the 4D theory changes from patch to patch. For this reason it is
not immediately obvious whether wavefunction continuity (particularly for
derivatives) must hold.

However, Ref. vBK  argued that continuity properties are still expected. In
brief, at weak coupling, all transverse modes are expected to be close to
their ground states at $r_{k}=\pi/L$, and this property was reflected in the
numerical analysis of vBK . At strong coupling, the transverse modes are
excited. However, they also mix, and we do not expect level crossing.
Continuity of the energy density then leads to the conditions Eq. (12). The
interpretation of the wavefunction simply changes discontinuously from patch
to patch.

Thus we are led to the quantum mechanics of three particles in the region Eq.
(11) subject to boundary conditions Eq. (12). The effective Hamiltonian can be
written as Eq. vBK

$H_{\mathrm{eff}}=-\frac{1}{2L}\bigg{(}\frac{1}{g^{2}}+\alpha_{1}\bigg{)}^{-1}\frac{\partial^{2}}{(\partial
c^{a}_{i})^{2}}+V_{T}(c)+V_{l}(c).$ (13)

Here $V_{T}$ is the “transverse” part of the effective potential which
vanishes along the vacuum valley. $V_{l}$ is the effective potential along the
vacuum valley. These potentials take the form:

$\displaystyle V_{T}$
$\displaystyle=\frac{1}{4}\big{(}\frac{1}{g^{2}}+\alpha_{2}\big{)}F^{a}_{ij}F^{a}_{ij}+\alpha_{3}(F^{a}_{ij}F^{a}_{ij})c^{b}_{k}c^{b}_{k}+\alpha_{4}F^{a}_{ij}F^{a}_{ij}c^{b}_{j}c^{b}_{j}+\alpha_{5}(\det
c)^{2}+...$ $\displaystyle V_{l}$
$\displaystyle=\gamma_{1}(g)\sum_{i}r^{2}_{i}+\gamma_{2}(g)\sum_{i}r^{4}_{i}+\gamma_{3}(g)\sum_{i>j}r^{2}_{i}r^{2}_{j}+\gamma_{4}(g)\sum_{i}r^{6}_{i}+\gamma_{5}\sum_{i\neq
j}r^{4}_{i}r^{2}_{j}+\gamma_{6}(g)r^{2}_{1}r^{2}_{2}r^{2}_{3}+\gamma_{7}(g)\sum_{i}r^{8}_{i}+...$
(14)

where $\dots$ indicates terms at higher orders in the fields and in the loop
expansion. We only consider the Hamiltonian with $\alpha_{1},\alpha_{2}\neq 0$
and $\gamma_{i}(g)=\gamma(0)$ for $i\leq 7$ and ignore higher order terms. The
effective theory is truncated at the two derivative order, so it will break
down for energies above the Kaluza-Klein scale $\sim 2\pi/L$ (in addition to
energies above the sphaleron energy, as discussed above.)

In order to extract the low-lying spectrum on a quantum simulator we use
Hamiltonian truncation and compute the matrix elements classically. The
kinetic part can be diagonalized by writing the momentum operator in spherical
coordinates for “particle” $i$ and using a spherical harmonic basis
$Y_{l_{i},m_{i}}(\theta_{i},\phi_{i})$ for each particle. For the radial part
of momentum operator we employ two choices of basis: spherical Bessels
$\chi^{(e)}_{n,l}=j_{l}(k^{e}_{n_{i},l_{i}}r)$ with $V(r)=0$ and a harmonic
oscillator basis $\chi^{(e)}_{n,l}=r^{l}e^{-\omega
r^{2}/2}M(\frac{1}{2}l+\frac{3}{4}-\epsilon^{(e)}_{n,l}/(2\omega),l+\frac{3}{2},\omega
r^{2})$ with $V(r)=\frac{1}{2}\omega^{2}r^{2}$. $M(a,b,z)$ is a confluent
hypergeometric function regular at $z=0$ vBK  and we have set $\omega=1.5$. In
either case the radial basis wavefunctions can be chosen to satisfy the Eq.
(12) in each flux-sector labelled by the flux quantum number $e_{i}$. This
boundary condition determines the “momenta”
$k^{e}_{n,l}(\epsilon^{(e)}_{n,l})$ for each particle in terms of the zeros of
the chosen basis functions.

In general, VQE performs better the more “physics” we can inject into the
choice of basis. At strong coupling, the wavefunctions tend to spread out
more, while for weak coupling they are more concentrated near the origin. For
this reason we use the spherical Bessel basis for strong coupling $g\geq 1.2$
and the harmonic oscillator basis for weak coupling $g<1.2$. (The spherical
bessel basis exactly diagonalizes the tree-level hamiltonian at infinite
coupling.) Of course, in the absence of truncation the spectrum is independent
of the choice of basis, but with truncations a better choice of basis can
minimize the error in eigenvalue computations.

Since we are using gauge-invariant coordinates on the vacuum valley, we can
work in a fully gauge-invariant hilbert space by requiring the angular
wavefunctions to be spin-singlets of $\mathrm{SO}(3)$. Hence the gauge
invariant wavefunctions have no free $m_{i}$ quantum numbers. (Therefore the
approach here differs from other quantum simulation studies of Yang-Mills type
matrix quantum mechanics models where simpler bases are used at the cost of
introducing an extended Hilbert space that includes gauge-non-invariant states
PRXQuantum.3.010324 .)

In sum, the gauge-invariant Rayleigh-Ritz basis for computing the full
Hamiltonian matrix are defined as follows:

$\displaystyle\innerproduct{\\{r_{i},\theta_{i},\phi_{i}\\}_{i=1,2,3}}{l_{1}l_{2}l_{3}n_{1}n_{2}n_{3};\mathbf{e}}=\sum_{m_{1},m_{2},m_{3}}W(l_{1}l_{2}l_{3}m_{1}m_{2}m_{3})\prod^{3}_{i=1}\chi^{e_{i}}_{n_{i},l_{i}}(r_{i})Y_{l_{i}m_{i}}(\theta_{i},\phi_{i})$
$\displaystyle n_{i},l_{i}\in\mathbb{N},\hskip 5.69046pt|l_{1}-l_{2}|\leq
l_{3}\leq l_{1}+l_{2},\hskip 5.69046ptm_{i}\in\\{-l_{i},-l_{i}+1,...,l_{i}\\}$

where $W(l_{1}l_{2}l_{3}m_{1}m_{2}m_{3})$ are the Wigner 3-$j$ symbols.

Now we must address truncation. We use eigenvalues of the operator
$\frac{\partial^{2}}{(\partial c^{a}_{i})^{2}}$ to organize the states in Eq.
(II) in an ascending order. The full Hamiltonian is infinite dimensional so we
truncate it to a finite number of states. We can compute the full matrix in
Eq. (13) for finite number of states by computing angular and radial matrix
elements in the Rayleigh-Ritz basis. The effective hamiltonian can be further
projected onto the irreducible representations of its symmetry group. For
electric flux sectors $\mathbf{e}=0$ or $\mathbf{e}=(1,1,1)$ the Hamiltonian
symmetry group is the full cubic group $O(3,\mathbb{Z})$. The cubic group is a
semidirect product $O(3,\mathbb{Z})=\mathbb{Z}^{3}_{2}\rtimes S_{3}$, where
the $\mathbb{Z}_{2}$ factor corresponds to parity transformations
$P_{i}c^{a}_{k}=-\delta_{ik}c^{a}_{k}$ and $S_{3}$ represents coordinate
permutations $\pi$. These symmetries acts on the gauge-invariant states in the
following manner:

$\displaystyle P_{i}\ket{l_{1}l_{2}l_{3}n_{1}n_{2}n_{3}}$
$\displaystyle=(-1)^{l_{i}}\ket{l_{1}l_{2}l_{3}n_{1}n_{2}n_{3}}$
$\displaystyle\pi\ket{l_{1}l_{2}l_{3}n_{1}n_{2}n_{3}}$
$\displaystyle=\ket{l_{\pi(1)}l_{\pi(2)}l_{\pi(3)}n_{\pi(1)}n_{\pi(2)}n_{\pi(3)}}$
(15)

The cubic group has ten irreps. We focus on the parity even irreps
$A^{+}_{1}$(zero flux) and $e^{+}_{1}$ (unit flux, $\mathbf{e}=(0,0,1)$). For
$\mathbf{e}\neq 0$, $\mathbf{e}\neq(1,1,1)$ the cubic group is broken to
$\mathbb{Z}_{2}\rtimes(\mathbb{Z}^{2}_{2}\rtimes S_{2})$ where $S_{2}$
permutes the directions with equal electric flux. We construct the
Hamiltonians classically in both sectors and numerically diagonalize them for
different numbers $M$ of states. Finally, for the case $M=8$ we use VQE and
QSE to compute the low-lying spectrum and compare with the classical exact
diagonalization.

The ground state in the $A^{+}_{1}$ sector is used as the reference ground
state for the system, identified with the true ground state in the large
volume limit. The gap to the first excitation above the ground state in the
$A^{+}_{1}$ sector is identified as a glueball mass. The difference in ground
state energies in the $e^{+}_{1}$ flux sector and the $A^{1}_{+}$ sector is
identified with the energy of an electric flux string. The spectrum of the
Hamiltonian in these irreps suffices to compute observables, like the square
root of the string tension to the glueball mass ratio, that can be compared
with values from Euclidean lattice simulations.

##  III Quantum subspace expansion for excited states

As mentioned above we use VQE to obtain ground states energies. In order to
obtain excited state energies, we apply the Quantum Subspace Expansion (QSE)
PhysRevA.95.042308  to the femtouniverse model. QSE is an extension to the VQE
algorithm based on quantum measurements of a set of ansatz excitation
operators in the optimal ground state estimate
$\ket{\psi(\boldsymbol{\theta^{*}})}$ obtained from VQE. Here
$\boldsymbol{\theta}=(\theta_{0},\theta_{1},...,\theta_{n_{p}})$ is the set of
parameters in the VQE ansatz and $\boldsymbol{\theta}^{*}$ denotes the optimal
parameters that minimize the energy.

The subspace in QSE is the Hilbert space spanned by a set of ansatz excitation
operators $\\{O_{i}\\}_{i}$ acting on $\ket{\psi(\boldsymbol{\theta}^{*})}$,
i.e., a set of states $\\{O_{i}\ket{\psi(\boldsymbol{\theta}^{*})}\\}_{i}$. To
include the ground state itself, we define $O_{0}=I$, the identity operator.
We denote the dimension of the subspace as $d^{\mathrm{QSE}}$. Then we can
evaluate the Hamiltonian on the subspace as

$H^{\mathrm{QSE}}_{ij}\equiv\matrixelement{\psi(\boldsymbol{\theta}^{*})}{O_{i}^{\dagger}HO_{j}}{\psi(\boldsymbol{\theta}^{*})}.$
(16)

We further require $O_{i}$ to be Hermitian operators. The matrix elements of
$H^{\mathrm{QSE}}$ can be measured on a quantum circuit, since the operator
$O_{i}^{\dagger}HO_{j}$ can be decomposed as a linear combination of Pauli
strings for the original Hilbert space. The basis
$\\{O_{i}\ket{\psi(\boldsymbol{\theta})}\\}_{i}$ is not orthonormal, so
instead of an ordinary eigenvalue problem for $H^{\mathrm{QSE}}$, we will need
to also compute the overlap matrix

$S^{\mathrm{QSE}}_{ij}\equiv\matrixelement{\psi(\boldsymbol{\theta}^{*})}{O_{i}^{\dagger}O_{j}}{\psi(\boldsymbol{\theta}^{*})}$
(17)

and solve a generalized eigenvalue problem (GEVP)

$H^{\mathrm{QSE}}\mathbf{v}=\lambda S^{\mathrm{QSE}}\mathbf{v}.$ (18)

The eigenvector $\mathbf{v}$ found this way does not directly give the
eigenwavefunctions but needs a further transformation explained in Appendix B.
For our purpose, we do not use this transformation and only focus on the
energies given by $\lambda$. Similar to $H^{\mathrm{QSE}}$, the matrix
elements of $S^{\mathrm{QSE}}$ can also be measured on the quantum circuit.
When the ansatz excitation operators are chosen effectively, solving the GEVP
yields a set of eigenstates and eigenenergies that are close to the lowest
$d^{\mathrm{QSE}}$ eigenstates and eigenenergies of $H$ in the original
Hilbert space.

In the example of 3 qubits with the truncation $M=2^{3}=8$, we choose
${O_{i}}$ as $O_{0}=III$, $O_{1}=XII$, $O_{2}=IXI$, $O_{3}=IIX$ in the Pauli-
string notation. This set of ansatz operators becomes approximately accurate
when using the spherical Bessel function basis in the large-$g$ limit, since
the lowest excitations are approximately the single-qubit-flipped states from
the ground state.

Measuring the matrix elements of $H^{\mathrm{QSE}}$ and $S^{\mathrm{QSE}}$
requires evaluations of Pauli strings expectation values at the state
$\ket{\psi(\boldsymbol{\theta}^{*})}$. Since the total number of Pauli strings
in the 3-qubit case is only $64$ and our Hamiltonian is relatively dense, we
measure all the $64$ Pauli strings. It is worth noting that the quantum
computational resource required at QSE in our case is significantly lower than
VQE, because the VQE does quantum measurements at multiple positions in the
space of $\boldsymbol{\theta}$, but in QSE, $\boldsymbol{\theta}$ is fixed to
$\boldsymbol{\theta}^{*}$.

The ground state energy obtained from QSE is not necessarily the same as the
ground state energy from the prior VQE. We include the results for ground
state energy from VQE, and the first excited state energies from the follow-up
QSE.

##  IV Results from numerical simulation and real quantum hardware

We implement our numerical simulation on the measurement-based Aer simulator
in Qiskit Qiskit , a python framework for quantum computation. We work in
units where the volume is $L=1$ and we scan over the running coupling $g(L)$,
which is equivalent to adjusting the dimensionless combination $\Lambda L$.

We also implement the VQE and QSE computation on the IBM quantum computer
Lima. In these computations we use the $M=8$ truncation with 3 qubits and test
three values of the coupling, $g\in\\{1.8,2.2,2.6\\}$. We use measurement
error mitigation https://doi.org/10.48550/arxiv.2010.08520  to reduce the
effect of noise on the real quantum hardware.

The effective Hamiltonian does not exhibit the canonical form of spin or
fermionic Hamiltonian. The simplest digital quantum encoding scheme involves
expanding the Hamiltonian in a Pauli string (tensor product of Paulis) basis.
This naive encoding, for a generic Hamiltonian, leads to number of Pauli
strings that is exponential in the number of qubits rather than polynomial.

We begin by performing multiple runs of noiseless simulated VQE with
$N_{\mathrm{shots}}=10{,}000$ quantum measurements for a given coupling,
starting from fully random initial point for every run. This is a test of
whether the ansatz is able to reliably find the correct result. In Fig. 2 and
Fig. 2 we show histograms of $A^{+}_{1}$ and $e^{+}_{1}$ values obtained from
100 VQE runs with $g=2.6$. The results are clustered around the exact
diagonalization value to within a few percent. The fact that VQE can
occasionally return results below the truth value is a consequence of shot
noise; in our simulations we use 10000 shots, so few-% errors of this type are
expected.

Figure 1: Histogram of 100 noiseless VQE runs vs classical exact
diagonalization result at $g=2.6$ for $A^{+}_{1}$ irrep on Aer simulator with
$N_{\mathrm{shots}}=10{,}000$.

Figure 2: Histogram of 100 noiseless VQE runs vs classical exact
diagonalization result at $g=2.6$ for $e^{+}_{1}$ irrep on Aer simulator with
$N_{\mathrm{shots}}=10{,}000$.

The fidelity between the initial ansatz state and the true ground state, the
flatness of energy landscape around local minima mcclean2018barren , and
limitations arising from the classical optimizer are the dominant factors
which affect the shape of these distributions. In particular the flatness of
energy landscape around a local minima can lead to a VQE result which is far
away from the rest of the values, even in noiseless simulation.

Figure 3: VQE vs classical exact diagonalization result for $A^{+}_{1}$ irrep.

Figure 4: VQE vs classical exact diagonalization result for $e^{+}_{1}$ irrep.

In Fig. 4 and Fig. 4 we show the noiseless VQE simulation results for a range
of $g$ values and compare to the ED values. We also show IBM-Lima results in
Fig. 4 and Fig. 4 with $M=8$ states(3-qubits) for $g=1.8,2.2,2.6$. We used the
measurement error mitigation https://doi.org/10.48550/arxiv.2010.08520  to
reduce the effect of noise on the real quantum hardware. We present both the
VQE minimum and median555We order the energies first. When $n_{\mathrm{runs}}$
of VQE is even, the $(n_{\mathrm{runs}}/2)$-th run (slightly skewed toward the
lower side) is chosen as the median, so that the median is always a single
sample, not an average of 2 samples. When $n_{\mathrm{runs}}$ of VQE is odd,
the usual median is taken. over 100 runs to emphasize the impact of performing
multiple runs.

For most values of $g$, we observe excellent agreement for the ground state
results obtained with high-truncation ED (1024 states), low truncation ED (8
states), and the noiseless VQE simulations with both the minimum and the
median taken over runs. However, a clear discontinuity is apparent at $g=1.2$,
and for $1.2\lesssim g\lesssim 1.5$, the spread in results is pronounced. This
discontinuity arises from the change in basis employed as we go from the weak
to the strong coupling regime. The discontinuity is invisible for $M=1024$,
indicating that with a large truncation cutoff the results are essentially
basis-independent, but the $M=8$ truncation shows a discontinuity. In this
case we also see that taking the minimum over VQE runs generally gets closer
to the exact value than the VQE median. The disagreement between median and
minimum is a consequence of significant numbers of runs in which the VQE
algorithm became stuck in spurious local minima, unable to escape due to the
flatness of the energy landscape around the local minima. In such cases the
minimum is the most reasonable estimate of the ground state energy. However,
the mininum can undershoot the exact result due to shot noise, so in some
circumstances the median may be preferable.

The IBM-Lima results for $g=1.8,2.2,2.6$ are generally around $\sim 20\%$
above the noiseless simulations and exact results, apart from the result in
the $A_{1}^{+}$ sector for $g=2.2$ which is significantly higher. The $g=2.2$
experiment was performed on the same quantum computer but on a different day
from the $g=1.8$ and $g=2.6$ experiments. The larger discrepancy from
noiseless simulations exhibited by the $g=2.2$ experiment may reflect a day-
to-day variation in the quality of the quantum computer in the lab and
calibration of the apparatus therein.

Fig. 6 shows the first excited state energy in the $A_{1}^{+}$ sector,
comparing the QSE+VQE result to exact diagonalization. Fig. 6 shows the
glueball mass $m_{0^{+}}$, which is the gap between this state and the ground
state $m_{0^{+}}=E_{1}(A_{1}^{+})-E_{0}(A_{1}^{+})$. The effect of the
truncation to $M=8$ states substantially overestimates the glueball mass in
the small $g$ region, relative to the $M=1024$ truncation. This situation
improves for larger values of $g$. We repeat the same procedure of multiple
runs using optimal parameters from VQE runs including IBM-Lima runs for
$g=1.8,2.2,2.6$. The minimum of all QSE runs is a fairly good estimate of the
excited state energy, whereas median values from QSE show large deviations
from the exact values. For this reason we show only the minimum results in
Fig. 6. These results converge well to the low-truncation exact
diagonalization values at both large and small $g$, and to the high-truncation
values at large $g$. In addition, the $m_{0^{+}}$ results on real hardware
agree well with the exact results, even in the worst case at $g=2.2$ ($\sim
30\%$), where the upward bias partially cancels in the energy difference.

In Fig. 8 we show the string tension $\sigma$ given by the energy gap between
the ground states of the $e_{1}^{+}$ and $A_{1}^{+}$ sectors, $\sigma
L=E_{0}(e_{1}^{+})-E_{0}(A_{1}^{+})$ (where computationally we work in units
with $L=1$.) The tension is highly suppressed for small couplings, reflecting
the exponential cost of tunneling through the quantum-induced barriers
separating center-conjugate vacua. The suppression disappears at larger
couplings as wavefunctions are more readily able to penetrate the barrier,
signalling the onset of large volume physics.

We emphasize that the largest values of $g$ we consider coincide with the
breakdown of the dimensionally reduced effective quantum mechanics as a good
description of the low-lying states of the full 4D theory. Two neglected
physical effects come into play around the same time: Kaluza-Klein masses
become comparable to the lightest glueball mass and string tension, as does
the sphaleron energy barrier separating sectors of varying Chern-Simons number
van2001qcd . These are large-volume effects that cannot be captured by the
matrix model effective theory.

Figure 5: QSE vs classical exact diagonalization result for $A^{+}_{1}$ irrep
excited state.

Figure 6: Lightest glueball mass in the $A_{1}^{+}$ sector.

Figure 7: String tension, $g\in[0.8,2.6]$.

Figure 8: Square root of string tension to glueball mass ratio vs
$\log(\Lambda L)$.

Nonetheless, as we approach $\Lambda L\sim 1$ from above, some of the large-
volume physics is quantitatively captured by the matrix model. In Fig. 8 we
plot the square root of the string tension-to-glueball mass ratio
$\sqrt{\sigma}/{m_{0^{+}}}$ against $-\log_{10}(\Lambda L)$, using the 2-loop
$\overline{MS}$ estimate for $\Lambda$, $\Lambda
L=e^{-\frac{1}{2\beta_{0}g^{2}}}(\beta_{0}g^{2})^{\frac{-\beta_{1}}{2\beta^{2}_{0}}}$
, where $\beta_{0}=\frac{N}{16\pi^{2}}\frac{11}{3}$ and
$\beta_{1}=(\frac{N}{16\pi^{2}})^{2}\frac{34}{3}$ with $N=2$. We have also
added the continuum, large-volume extrapolation from lattice results Teper  as
a solid horizontal line. We see that the simulations converge to within about
$15\%$ of the continuum, large-volume value for this observable in the largest
volumes (left-hand side of the plot), and agreement is found using IBM-Lima
for the largest value of $g$. In Fig. 9 we zoom in on the large $g$ region of
this plot.

Figure 9: Square root of string tension to glueball mass ratio.

##  V Conclusions and outlook

In the near-term, quantum simulations of the type performed here will be
limited by noise. We find, for example, that IBM Lima does not give stable
results with 4 qubits for our Hamiltonian. This is mainly because of the
increased circuit depth and thus the accumulated error from noise. However,
there are several directions for improvement that could help push the reach of
VQE simulations future, including device-specific choices of gates, optimizing
the complexity of the ansatz, quantum error correction PhysRevA.32.3266 , and
theoretical improvements in the algorithms, such as utilizing commuting
families of Pauli strings Gokhale2019MinimizingSP  to reduce the number of
measurements. Recent developments of Koopman operator learning techniques for
quantum optimization luo2022koopman  can accelerate VQE when using gradient-
based optimizers, saving quantum resources for other use.

Compared to the lattice approaches, the femtouniverse, as an example of
effective theory of the low-momentum modes, represents a different step toward
ultimate target of quantum simulations for high energy physics. Using a
relatively small number of qubits on a quantum computer, we can already
simulate some 4D physics and obtain qualitatively reasonable results in
agreement with 4D lattice simulations on classical computers. The dimensional
reduction approach supplies a different type of regulator from the lattice
regulator, and is relevant to other models of interest in high energy physics
including models of quantum gravity bfss . The manifest gauge invariance in
the quantum basis also helps reduce the number of dimension of the Hilbert
space and thus the number of qubits.

In this paper, using VQE and QSE, we find that the currently available noisy
quantum computers can already produce the glueball mass and string tension for
the 4D SU(2) Yang-Mills theory qualitatively in agreement with other
approaches. Our approach only requires a very small number of qubits, but the
Hamiltonian matrix has a relatively large number of elements, and thus may
require a relatively deep quantum circuit and big number of Pauli strings. The
encoding of the basis states into the qubit computational basis is not fine
tuned, and Hamiltonian matrix elements are all kept without specific ordering
or approximation. In addition, with more delicate choice of the representation
of the femtouniverse on the quantum computer, it will be interesting to study
if the depth of the ansatz can be reduced without loss of essential physics.

The $\mathrm{SU}(2)$ Yang-Mills matrix quantum mechanics model studied here
can be extended in a few different directions. Adding higher-momentum modes
and their interaction with the zero-momentum modes can improve the
approximation to the complete 4D Yang-Mills theory, while adding fermions can
connect to more physical models of interest. Future study on how the
complexity and performance scales with $N$ for the $\mathrm{SU}(N)$ Yang-Mills
theory will likewise inform the scalability and generalizability of the
approach.

### Acknowledgments

We thank Aida El-Khadra and Di Luo for helpful discussions. This work was
supported in part by the U.S. Department of Energy, Office of Science, Office
of High Energy Physics under award number DE-SC0015655 and by its QuantISED
program under an award for the Fermilab Theory Consortium “Intersections of
QIS and Theoretical Particle Physics.” We acknowledge the use of IBM Quantum
services for this work. The views expressed are those of the authors, and do
not reflect the official policy or position of IBM or the IBM Quantum team. In
this paper we used ibmq_lima, which is one of the IBM Quantum Falcon
Processors.

##  Appendix A ED results

Table 1 contains exact diagonalization results for ground state energies of
the irreps $A_{1}^{+}$, $e^{+}_{1}$ and the first excited energy for
$A^{+}_{1}$ denoted by $(A^{+}_{1})^{{}^{\prime}}$.

|  | $A_{1}^{+}$ |  |  | $(A^{+}_{1})^{{}^{\prime}}$ |  |  | $e^{+}_{1}$ |
---|---|---|---|---|---|---|---|---|---
$g$ | 1024 | 16 | 8 | 1024 | 16 | 8 | 1024 | 16 | 8
0.8 | 2.580 | 2.617 | 2.638 | 3.540 | 4.183 | 4.452 | 2.585 | 2.633 | 2.645
0.9 | 2.682 | 2.732 | 2.752 | 3.651 | 4.303 | 4.562 | 2.697 | 2.758 | 2.772
1.0 | 2.752 | 2.821 | 2.844 | 3.763 | 4.486 | 4.722 | 2.786 | 2.865 | 2.919
1.1 | 2.792 | 2.885 | 2.921 | 3.878 | 4.735 | 4.985 | 2.854 | 2.962 | 3.066
1.2 | 2.806 | 2.966 | 3.332 | 3.999 | 4.986 | 5.984 | 2.903 | 3.151 | 3.439
1.3 | 2.799 | 2.923 | 3.222 | 4.134 | 4.844 | 5.743 | 2.941 | 3.107 | 3.352
1.4 | 2.776 | 2.868 | 3.109 | 4.285 | 4.792 | 5.676 | 2.971 | 3.078 | 3.300
1.5 | 2.739 | 2.803 | 2.995 | 4.454 | 4.804 | 5.672 | 2.995 | 3.062 | 3.256
1.6 | 2.690 | 2.731 | 2.883 | 4.641 | 4.872 | 5.679 | 3.018 | 3.058 | 3.216
1.8 | 2.565 | 2.578 | 2.668 | 5.074 | 5.164 | 5.768 | 3.068 | 3.081 | 3.165
2.0 | 2.412 | 2.417 | 2.462 | 5.590 | 5.622 | 6.028 | 3.141 | 3.145 | 3.182
2.2 | 2.251 | 2.253 | 2.273 | 6.194 | 6.206 | 6.463 | 3.256 | 3.258 | 3.273
2.4 | 2.101 | 2.102 | 2.110 | 6.884 | 6.889 | 7.045 | 3.428 | 3.429 | 3.435
2.6 | 1.980 | 1.981 | 1.985 | 7.652 | 7.654 | 7.746 | 3.663 | 3.664 | 3.667
Table 1: Ground state and the first excited state energies of the irrep
$A^{+}_{1}$, and the ground state energies of the irrep $e_{1}^{+}$ from exact
diagonalization with truncations $M=1024,16,8$. As described in the text, we
change the basis used for the truncation at $g=1.2$.

##  Appendix B Details in VQE and QSE

The ansatz wavefunction for the ground state of the Hamiltonian $H$ can be
constructed on a quantum circuit as
$\ket{\psi(\boldsymbol{\theta})}=U(\boldsymbol{\theta})\ket{0}$, where
$U(\boldsymbol{\theta})$ is a unitary operator consisting of quantum gates
parametrized by
$\boldsymbol{\theta}=(\theta_{0},\theta_{1},...,\theta_{n_{p}})$, and
$\ket{0}$ is the default starting state on the quantum circuit. The target
energy function to be minimized is

$E(\boldsymbol{\theta})\equiv\matrixelement{\psi(\boldsymbol{\theta})}{H}{\psi(\boldsymbol{\theta})}.$
(19)

$\boldsymbol{\theta}^{*}\equiv\mathrm{argmin}_{\boldsymbol{\theta}}E(\boldsymbol{\theta})$
gives an approximate ground state of $H$, and the energy evaluated at
$\boldsymbol{\theta}^{*}$ gives the approximate ground state energy. The
Hamiltonian is encoded into the quantum circuit by a linear decomposition in
terms of Pauli strings. We also refer to the review tilly2022variational  for
more details of VQE.

In VQE computations we use the RealAmplitudes ansatz with linear entanglement
and $\mathrm{reps}=2$. There are $9$ parameters in the ansatz. RealAmplitudes
only gives a real-valued wavefunction, and the eigenstates of the Hamiltonian
we consider can always be constrained to be real-valued, so this choice of
ansatz is reasonable for our task with a relatively small number of
parameters. We use the optimizer COBYLA (powell1994direct, ) with
$\mathrm{maxiter}=100$ (maximum number of iterations), $\mathrm{rhobeg}=0.5$
(reasonable initial changes to the variables), $\mathrm{tol}=0.001$ (final
accuracy in the optimization). The number of shots for quantum measurement is
$10{,}000$ so that the statistical error associated with measurement is small.
We also use the measurement error mitigation with $10{,}000$ shots when
constructing the calibration matrix.

When implementing QSE, we make measurements of all $4^{3}=64$ Pauli strings
for the $3$-qubit system, based on the same ansatz as VQE at
$\boldsymbol{\theta}^{*}$. We again use $10{,}000$ shots for each Pauli string
evaluation and also the measurement error mitigation.

We denote the nonorthogonal basis by $\\{\ket{\phi_{i}}\\}$
($i=0,1,...,d^{\mathrm{QSE}}-1$ where $d^{\mathrm{QSE}}=4$ in our case).

The transformation between $\\{\ket{\phi_{i}}\\}$ and an orthonormal basis
$\\{\ket{\psi_{i}}\\}$ is

$\ket{\psi_{i}}=P^{*}_{ij}\ket{\phi_{j}}.$ (20)

Then the following operator identity holds

$PSP^{\dagger}=I,$ (21)

where

$S_{ij}=\innerproduct{\phi_{i}}{\phi_{j}}$ (22)

is the overlap matrix. In the QSE basis, $S=S^{\mathrm{QSE}}$ in the main
text. An operator $O$ when acting on the nonorthonormal basis transforms as

$O^{\prime}=SP^{\dagger}OPS.$ (23)

With the Hamiltonian matrix in the nonorthogonal basis

$H^{\mathrm{QSE}}_{ij}=\matrixelement{\phi_{i}}{H}{\phi_{j}}$ (24)

we need to solve the GEVP

$H^{\mathrm{QSE}}\mathbf{v}=\lambda S\mathbf{v}.$ (25)

We can transform the eigenvector $\mathbf{v}$ back by

$\mathbf{u}=PS\mathbf{v},$ (26)

where $\mathbf{u}$ is the wavefunction of the eigenstate in the original
orthonormal basis $\\{\ket{\psi_{i}}\\}$.

## References

  * [1] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [2] J. D. Bjorken.  Elements of Quantum Chromodynamics.  Prog. Math. Phys., 4:423–561, 12 1979.
  * [3] Gerard ’t Hooft.  A Property of Electric and Magnetic Flux in Nonabelian Gauge Theories.  Nucl. Phys. B, 153:141–160, 1979.
  * [4] M. Luscher.  Some Analytic Results Concerning the Mass Spectrum of Yang-Mills Gauge Theories on a Torus.  Nucl. Phys. B, 219:233–261, 1983.
  * [5] M. Lüscher and G. Münster.  Weak-coupling expansion of the low-lying energy values in the su(2) gauge theory on a torus.  Nuclear Physics B, 232(3):445–472, 1984.
  * [6] Pierre van Baal.  On the ratio of the string tension and the glueball mass squared in the continuum.  Nuclear Physics B, 264:548–556, 1986.
  * [7] Pierre Van Baal.  Gauge theory in a finite volume.  Technical report, 1988.
  * [8] Pierre van Baal.  Beyond the Gribov Horizon in the Femto Universe, pages 131–152.  Springer US, Boston, MA, 1990.
  * [9] Jeffrey Koller and Pierre van Baal.  A non-perturbative analysis in finite volume gauge theory.  Nuclear Physics B, 302(1):1–64, 1988.
  * [10] Pierre van Baal and Jeffrey Koller.  Qcd on a torus, and electric flux energies from tunneling.  Annals of Physics, 174(2):299–371, 1987.
  * [11] Pierre Van Baal.  Qcd in a finite volume.  In At The Frontier of Particle Physics: Handbook of QCD (In 3 Volumes), pages 683–760. World Scientific, 2001.
  * [12] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature communications, 5(1):1–7, 2014.
  * [13] Jules Tilly, Hongxiang Chen, Shuxiang Cao, Dario Picozzi, Kanav Setia, Ying Li, Edward Grant, Leonard Wossnig, Ivan Rungger, George H Booth, et al.  The variational quantum eigensolver: a review of methods and best practices.  Physics Reports, 986:1–128, 2022.
  * [14] Jarrod R. McClean, Mollie E. Kimchi-Schwartz, Jonathan Carter, and Wibe A. de Jong.  Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states.  Phys. Rev. A, 95:042308, Apr 2017.
  * [15] Michael J. Teper.  Glueball masses and other physical properties of SU(N) gauge theories in D = (3+1): A Review of lattice results for theorists.  12 1998.
  * [16] Hrant Gharibyan, Masanori Hanada, Masazumi Honda, and Junyu Liu.  Toward simulating superstring/m-theory on a quantum computer.  Journal of High Energy Physics, 2021(7):1–56, 2021.
  * [17] Enrico Rinaldi, Xizhi Han, Mohammad Hassan, Yuan Feng, Franco Nori, Michael McGuigan, and Masanori Hanada.  Matrix-model simulations using quantum computing, deep learning, and lattice monte carlo.  PRX Quantum, 3:010324, Feb 2022.
  * [18] Jiayu Shen, Di Luo, Chenxi Huang, Bryan K. Clark, Aida X. El-Khadra, Bryce Gadway, and Patrick Draper.  Simulating quantum mechanics with a $\theta$-term and an ’t hooft anomaly on a synthetic dimension.  Phys. Rev. D, 105:074505, Apr 2022.
  * [19] G.’t Hooft.  Naturalness, Chiral Symmetry, and Spontaneous Chiral Symmetry Breaking, pages 135–157.  Springer US, Boston, MA, 1980.
  * [20] MD SAJID ANIS, Abby-Mitchell, Héctor Abrahamand, et al.  Qiskit: An open-source framework for quantum computing, 2021.
  * [21] George S. Barron and Christopher J. Wood.  Measurement error mitigation for variational quantum algorithms, 2020\.
  * [22] Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  Nature communications, 9(1):1–6, 2018.
  * [23] Asher Peres.  Reversible logic and quantum computers.  Phys. Rev. A, 32:3266–3276, Dec 1985.
  * [24] Pranav Gokhale et. al.  Minimizing state preparations in variational quantum eigensolver by partitioning into commuting families.  arXiv: Quantum Physics, 2019.
  * [25] Di Luo, Jiayu Shen, Rumen Dangovski, and Marin Soljačić.  Koopman operator learning for accelerating quantum optimization and machine learning.  arXiv preprint arXiv:2211.01365, 2022.
  * [26] T. Banks, W. Fischler, S. H. Shenker, and L. Susskind.  $m$ theory as a matrix model: A conjecture.  Phys. Rev. D, 55:5112–5128, Apr 1997.
  * [27] Michael JD Powell.  A direct search optimization method that models the objective and constraint functions by linear interpolation.  In Advances in optimization and numerical analysis, pages 51–67. Springer, 1994.


disposition section title

DESY 18-043

Five-Particle Phase-Space Integrals in QCD

O. Gituliar

<EMAIL_ADDRESS>
V. Magerya

<EMAIL_ADDRESS>
A. Pikelner111On leave of absence from Joint Institute for Nuclear Research,
141980 Dubna, Russia.

<EMAIL_ADDRESS>
II. Institut für Theoretische Physik, Universität Hamburg,
Luruper Chaussee 149, D-22761 Hamburg, Germany

###### Abstract

We present analytical expressions for the 31 five-particle phase-space master
integrals in massless QCD as an $\epsilon$-series with coefficients being
multiple zeta values of weight up to 12. In addition, we provide computer code
for the Monte-Carlo integration in higher dimensions, based on the RAMBO
algorithm, that has been used to numerically cross-check the obtained results
in 4, 6, and 8 dimensions.

##  1 Introduction

Nowadays, perturbative calculations play the key role in describing data from
high-energy particle colliders, such as the LHC, as well as in improving the
precision of numerical parameters in Standard Model and other models. It is
clear now that higher-order calculations will play an even more crucial role
in processing data from future high-luminosity colliders, like the FCC or the
ILC, where theoretical errors will dominate over experimental statistical
errors. These arguments motivate us to make one step forward beyond available
fully- inclusive phase-space integrals for a four-particle decay [GGH04] and
calculate a set of yet unknown integrals that corresponds to a five-particle
decay of a color-neutral off-shell particle in Quantum Chromodynamics. These
results, among others, are necessary ingredients in calculating various
exclusive quantities with the method of differential equations where they are
needed to determine integration constants, as for example discussed in [Git16,
GM15] for the three-loop time-like splitting functions [AMV12].

In this paper, we focus on the calculation of master integrals that can be
used to express any other integral of the corresponding topology provided a
set of integration-by-parts rules (IBP) [CT81] is known. Our approach is based
on techniques for solving dimensional recurrence relations (DRR) [Tar96]
described in [LM17a, LM17b]. In particular, we use DREAM package [LM17a] to
obtain numerical results for desired integrals with 2000-digit precision and
restore analytical form in terms of multiple zeta values (MZV) [Fur03, BBV10,
LM16] up to weight 12 using the PSLQ method [FBA99] as implemented in
Mathematica. We also present a Monte-Carlo code, based on the RAMBO algorithm
[KSE86], for numerical integration of the phase-space integrals in arbitrary
(integer) number of dimensions that has been used to check consistency of the
obtained results.

This paper is organized as following. In Section 2 we introduce our notation
and describe our calculational method in more details. In Section 3 we provide
complete results for four-particle integrals and discuss numerical cross-
checks using Monte-Carlo integration. In Section 4 we make final remarks. In
Appendix A we provide the complete list of master integrals.

Additionally we provide auxiliary files on arXiv222https://arxiv.org
containing the complete master integrals with MZV weight up to 12 as well as
the Monte-Carlo integration routines with the corresponding results.

##  2 The Method

& missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{1}$ $F_{2}$: 12 14 23 34 $F_{3}$: 01 05 $F_{4}$: 12 13 02 03 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{5}$: 12 34 02 03 $F_{6}$: 12 13 02 04 45 35 $F_{7}$: 123 124 $F_{8}$: 13
14 123 124 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{9}$: 12 23 023 012 $F_{10}$: 01 123 $F_{11}$: 01 123 014 05 $F_{12}$: 12
023 05 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{13}$: 12 13 124 05 02 $F_{14}$: 02 12 023 05 $F_{15}$: 12 13 124 134 02 03
$F_{16}$: 12 34 02 03 124 134 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{17}$: 12 34 013 $F_{18}$: 12 13 25 124 05 04 $F_{19}$: 12 13 25 124 013 05
$F_{20}$: 12 13 124 134 35 25 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{21}$: 12 13 24 012 05 $F_{22}$: 12 13 24 012 45 05 $F_{23}$: 12 13 25 45
134 012 $F_{24}$: 12 13 45 134 05 03 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{25}$: 12 13 45 134 012 05 $F_{26}$: 012 123 $F_{27}$: 12 01 123 013
$F_{28}$: 123 012 05 01 missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

missing

\makebox(0.0,0.0)[c]{\footnotesize$$}

$F_{29}$: 12 123 25 013 05 01 $F_{30}$: 12 13 034 134 05 03 $F_{31}$: 12 13
034 134 05 04

Table 1: Cut diagrams for five-particle phase-space master integrals in QCD.
Dashed lines represent cut propagators and carry final-state momenta
$p_{1},\dots,p_{5}$. Labels represent propagators, so that "123" corresponds
to $p_{1}+p_{2}+p_{3}$ and "012" to $q-p_{1}-p_{2}$ (where $q$ is the initial-
state momentum, i.e., $q=p_{1}+\dots+p_{5}$).

We start by identifying a set of five-particle phase-space master integrals by
constructing auxiliary topologies containing the four-loop massless
propagators from [BC10], taking all physical five-particle cuts of those, and
performing the Laporta-style IBP reduction [Lap00, Lap17] implemented in FIRE5
[Smi15]. This gives us 31 master integrals listed in Table 2 with up to 6
unique propagators each. Our notation is

$F_{i}=S_{\Gamma}\int\textrm{d}\textrm{PS}_{5}~{}\frac{1}{D^{(i)}_{1}~{}\dots~{}D^{(i)}_{n}},$
(1)

where $D^{(i)}_{j}$ are propagators that take the form of invariant scalar
products

$s_{kl\dots q}=\left(p_{k}+p_{l}+\dots+p_{q}\right)^{2},$ (2)

$\textrm{d}\textrm{PS}_{5}$ is a five-particle phase-space element in $D$
dimensions

$\textrm{d}\textrm{PS}_{N}=\left(\prod_{i=1}^{N}\textrm{d}^{D}p_{i}\,\delta^{+}\big{(}p_{i}^{2}\big{)}\right)\delta^{(D)}\Big{(}q-p_{1}-\ldots-
p_{N}\Big{)},$ (3)

and $S_{\Gamma}$ is a common normalization factor chosen for convenience to be

$S_{\Gamma}=\left(q^{2}\right)^{5-2D}\frac{\big{(}2\pi\big{)}^{4}}{\pi^{2D}}\Gamma\left(\frac{D}{2}-1\right)\Gamma\left(3\frac{D}{2}-3\right).$
(4)

With this normalization and knowing the volume of the complete $N$-particle
phase space333The dependence on $q^{2}$ here is trivial, and can be restored
by power counting. We will omit it from now on, setting $q^{2}$ to $1$.

$\int\textrm{d}\textrm{PS}_{N}=\left(q^{2}\right)^{\frac{D}{2}(N-1)-N}\,\frac{\pi^{\frac{D}{2}(N-1)}}{\big{(}2\pi\big{)}^{N-1}}\frac{\Gamma\left(\frac{D}{2}-1\right)^{N}}{\Gamma\Big{(}\left(\frac{D}{2}-1\right)\big{(}N-1\big{)}\Big{)}\,\Gamma\Big{(}\left(\frac{D}{2}-1\right)N\Big{)}}$
(5)

we can already fix the value of $F_{1}$ as

$F_{1}=S_{\Gamma}\int\textrm{d}\textrm{PS}_{5}=\frac{\Gamma\big{(}\frac{D}{2}-1\big{)}^{6}\,\Gamma\big{(}3\frac{D}{2}-3\big{)}}{\Gamma\big{(}4\frac{D}{2}-4\big{)}\,\Gamma\big{(}5\frac{D}{2}-5\big{)}}.$
(6)

Next, with the help of LiteRed [Lee14] and FIRE5 [Smi15] we derive a set of
lowering dimensional recurrence relations which express master integrals in
$D+2$ dimensions in terms of master integrals in $D$ dimensions:

$F_{i}(D+2)=M_{ij}(D)\,F_{j}(D).$ (7)

In our case $M$ can be shuffled into triangular form, with each $F_{i}$ only
depending on itself and master integrals from lower sectors $S_{i}$:

$F_{i}(D+2)=M_{ii}(D)\,F_{i}(D)\,+\sum_{k\in S_{i}}M_{ik}(D)F_{k}(D)$ (8)

and the general solution being

$F_{i}(D)=\omega_{i}(D)\,H_{i}(D)+R_{i}(D),$ (9)

where $H_{i}$ is a homogeneous solution of eq. (8), $R_{i}$ is a partial
solution that can be constructed numerically with DREAM [LM17a] provided
$F_{1}$ is known, and $\omega_{i}(D)$ is an arbitrary periodic function that
needs to be determined from separate considerations.

We argue that all $\omega_{i>1}$ are zero. To see this, first let us look at
the asymptotic behavior of $F_{i}$ at large $D$. Rewriting eq. (1) as an
integral over invariants $s_{ij}$ gives

$F_{i}=S_{\Gamma}\left(\prod_{k=1}^{N-1}\Omega_{D-k}\right)\int\left(\prod_{l<m}\textrm{d}s_{lm}\right)\big{(}{\Delta_{N}}\big{)}^{\frac{D-N-1}{2}}\Theta\big{(}\Delta_{N}\big{)}\,\delta\left(1-s_{1\dots
N}\right)\frac{1}{D^{(i)}_{1}~{}\dots~{}D^{(i)}_{n}},$ (10)

where $\Delta_{N}$ is the Gram determinant defined as

$\Delta_{N}=\frac{(-1)^{N+1}}{2^{N}}\begin{vmatrix}s_{11}&s_{12}&\cdots&s_{1N}\\\
s_{12}&s_{22}&\cdots&s_{2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
s_{1N}&s_{2N}&\cdots&s_{NN}\end{vmatrix},$ (11)

and $\Omega_{k}$ is the surface of a unit hypersphere in $k$-dimensional space

$\Omega_{k}=2\pi^{\frac{k}{2}}\,\Gamma\left(\frac{k}{2}\right)^{-1}.$ (12)

If $\Delta_{N}(s_{ij})$ has a unique global maximum inside the integration
region, we can apply Laplace’s method to eq. (10) and find its asymptotic as

$F_{i}(D\to\infty)=S_{\Gamma}\left(\prod_{k=1}^{N-1}\Omega_{D-k}\right)\left(\Delta^{max}_{N}\right)^{\frac{D}{2}}\left(\frac{2\pi}{D}\right)^{\frac{1}{2}\left(\frac{N(N-1)}{2}-1\right)}\left(\mathcal{C}_{i}+\mathcal{O}\left(D^{-1}\right)\right),$
(13)

where $\mathcal{C}_{i}$ is a constant that depends on the location of the
maximum and the denominators $D^{(i)}_{j}$, but not on $D$.

The global maximum of $\Delta_{N}$ is reached when all $s_{ij}$ ($i\neq j$)
are identical and equal to $\frac{2}{N(N-1)}$. Geometrically this
configuration corresponds to the vectors $\vec{p}_{i}$ pointing to the
vertices of a regular $N$-hedron embedded into Euclidean space of $(N-1)$
dimensions. The maximum value is then

$\Delta^{max}_{N}=\frac{1}{N^{N}(N-1)^{N-1}}$ (14)

and explicitly we get

$F_{i}(D\to\infty)=\pi^{\frac{7}{2}}2^{\frac{25}{2}}\frac{\left(4^{4}5^{5}\right)^{-\frac{D}{2}}\Gamma\left(\frac{3D}{2}-3\right)}{D^{\frac{9}{2}}\Gamma\left(\frac{D-4}{2}\right)\Gamma\left(\frac{D-3}{2}\right)\Gamma\left(\frac{D-1}{2}\right)}\left(\mathcal{C}_{i}+\mathcal{O}\left(D^{-1}\right)\right).$
(15)

It follows that all $F_{i}$ have identical asymptotic behavior up to a
constant $\mathcal{C}_{i}$. As a confirmation, it can be shown that eq. (15)
is asymptotically the same expression as we had for $F_{1}$ in eq. (6).

Next, we can find the asymptotics of the homogeneous parts of eq. (9),
$H_{i}(D)$, using e.g. the routine FindAsymptotics from DREAM. Comparing these
to eq. (15), we determine that all $H_{i}(D)$ for $i>1$ are growing
exponentially faster than $F_{i}(D)$, which can only happen if the
corresponding periodic functions $\omega_{i}(D)$ are zero.

Thus, to find $F_{i}$ we only need to find $R_{i}$, the inhomogeneous
solutions to eq. (8). We compute them as a series in $\epsilon=(4-D)/2$ using
DREAM with 2000-digit accuracy and then restore the analytical form of the
series coefficients in terms of MZVs using PSLQ method [FBA99]. In this way we
obtain the analytical result for all master integrals up to MZVs of weight 12
using the corresponding bases from [Fur03] and the SummerTime package [LM16]
for their numerical evaluation. Corresponding expressions are presented in
Appendix A as well as in the auxiliary files on arXiv.

##  3 Crosschecks

###  3.1 Four-Particle Integrals

As the first consistency check of our method we reproduce results for four-
particle phase-space integrals reported in [GGH04]. We perform all the steps
described in Section 2. Generating the IBP rules with the help of LiteRed and
then proceeding with DREAM we obtain the final result with 2000-digit accuracy
and MZVs up to weight 12. The series reconstructed with PSLQ (using the
original notation, and omitting $S_{\Gamma}$ and $q^{2}$ factors) are:

$\displaystyle R_{6}$
$\displaystyle=-1+\zeta_{2}+\epsilon\Bigg{(}-12+5\zeta_{2}+9\zeta_{3}\Bigg{)}+\epsilon^{2}\Bigg{(}-91+27\zeta_{2}+45\zeta_{3}+\frac{61}{5}\zeta_{2}^{2}\Bigg{)}$
(16)
$\displaystyle+\epsilon^{3}\Bigg{(}-558+161\zeta_{2}+197\zeta_{3}+61\zeta_{2}^{2}-80\zeta_{3}\zeta_{2}+207\zeta_{5}\Bigg{)}$
$\displaystyle+\epsilon^{4}\Bigg{(}-3025+939\zeta_{2}+897\zeta_{3}+\frac{1157}{5}\zeta_{2}^{2}-400\zeta_{3}\zeta_{2}+1035\zeta_{5}+\frac{288}{5}\zeta_{2}^{3}-153\zeta_{3}^{2}\Bigg{)},$
$\displaystyle R_{8,a}$
$\displaystyle=\frac{5}{\epsilon^{4}}-\frac{40\zeta_{2}}{\epsilon^{2}}-\frac{126\zeta_{3}}{\epsilon}+14\zeta_{2}^{2}+\epsilon\Bigg{(}1008\zeta_{2}\zeta_{3}-1086\zeta_{5}\Bigg{)}+\epsilon^{2}\Bigg{(}-\frac{272}{7}\zeta_{2}^{3}+1602\zeta_{3}^{2}\Bigg{)},$
(17) $\displaystyle R_{8,b}$
$\displaystyle=\frac{3}{4\epsilon^{4}}-\frac{17\zeta_{2}}{2\epsilon^{2}}-\frac{44\zeta_{3}}{\epsilon}-\frac{183}{5}\zeta_{2}^{2}+\epsilon\Bigg{(}376\zeta_{2}\zeta_{3}-790\zeta_{5}\Bigg{)}+\epsilon^{2}\Bigg{(}-\frac{19088}{105}\zeta_{2}^{3}+698\zeta_{3}^{2}\Bigg{)}.$
(18)

###  3.2 Numerical Verification

$i$ | Numerical results | Analytic results
---|---|---
| $D=4$ | $D=6$ | $D=8$ | $D=4$ | $D=6$ | $D=8$
2 | – | 1708(2)00 | 4699(1)0 | – | 171085.62 | 47000.531
3 | 3.7823(4) | 3.1704(2) | 3.0221(1) | 3.7823736 | 3.1704486 | 3.0221118
4 | – | 1504.7(8) | 725.3(1) | – | 1504.4507 | 725.26806
5 | – | 1007.4(5) | 580.80(9) | – | 1007.5235 | 580.76347
6 | – | 6191(5)00 | 14496(6)0 | – | 619633.25 | 144975.32
7 | 46.46(4) | 18.533(2) | 15.205(1) | 46.435253 | 18.532303 | 15.205538
8 | – | 2031(2)0 | 5357(2) | – | 20297.189 | 5355.3611
9 | – | 4313(3) | 2406.7(4) | – | 4312.8823 | 2406.7943
10 | 10.436(2) | 7.1508(5) | 6.5093(3) | 10.435253 | 7.1507477 | 6.5092878
11 | 228.8(1) | 62.67(1) | 47.663(4) | 229.11836 | 62.667046 | 47.663194
12 | – | 157.34(4) | 102.26(1) | – | 157.33521 | 102.26408
13 | – | 13729(8) | 4000(1) | – | 13732.166 | 4000.2779
14 | – | 268.46(8) | 172.80(2) | – | 268.45969 | 172.79805
15 | – | 6322(6)0 | 16048(5) | – | 63316.356 | 16049.857
16 | – | 4414(3)0 | 12952(4) | – | 44117.898 | 12951.443
17 | – | 1243.4(7) | 709.6(1) | – | 1243.1369 | 709.52840
18 | – | 3002(2)00 | 5899(3)0 | – | 300402.99 | 58965.517
19 | – | 4982(4)00 | 10637(4)0 | – | 498329.79 | 106357.81
20 | – | 2360(2)000 | 5777(2)00 | – | 2362594.9 | 577686.64
21 | – | 6312(5)0 | 20642(5) | – | 63147.876 | 20642.071
22 | – | 8402(7)00 | 24407(8)0 | – | 840453.94 | 244075.75
23 | – | 1443(1)000 | 4556(1)00 | – | 1443198.3 | 455543.43
24 | – | 1391(1)00 | 3997(1)0 | – | 139263.92 | 39966.878
25 | – | 3347(3)00 | 8526(3)0 | – | 335128.10 | 85254.217
26 | 25.563(6) | 15.376(1) | 13.6042(7) | 25.564747 | 15.376404 | 13.604247
27 | – | 697.3(3) | 397.84(6) | – | 697.18948 | 397.83514
28 | 143.9(1) | 52.855(7) | 42.917(3) | 143.97886 | 52.853837 | 42.917424
29 | – | 4409(3)0 | 13702(3) | – | 44117.898 | 13700.597
30 | – | 6327(6)0 | 16181(5) | – | 63316.356 | 16178.566
31 | – | 8955(8)0 | 19055(8) | – | 89611.062 | 19051.115
Table 2: Numerical results for the ratio $F_{i}/F_{1}$ with the corresponding
uncertainties (standard deviations) indicated in the parenthesis. Missing
entries correspond to divergent integrals.

As another cross-check we have calculated the leading terms of the
$\epsilon$-expansion of $F_{i}$ numerically using the direct way: through
Monte-Carlo integration of eq. (1) over the phase-space. While such a
technique can not be easily applied to divergent integrals, we can sidestep
such an issue by noting that our master integrals only suffer from IR
divergences that disappear already at $D=6$. In this way we can check several
leading terms of the expansion at $D=4-2\epsilon$ by calculating the
corresponding integrals in $D\geq 6$ since both are connected by dimensional
recurrence relations.

To calculate a finite integral of the form eq. (1) we choose a uniform mapping
from a hypercube into momentum coordinates using an algorithm similar to RAMBO
[KSE86] but extended into arbitrary $D$. Then we calculate the integrand from
scalar products of the momenta, and finally we integrate over the hypercube
using the Vegas [Lep78] implementation from Cuba [Hah05].

Note that although the integrals we are calculating are finite, the integrands
are not. Exposing an integration algorithm like Vegas to such infinities may
lead to unpredictable behaviour, so as a precaution we choose to regulate
these infinities by adding a small parameter $\alpha$ to the denominator of
the integrand, and then to calculate the integral with progressively smaller
values of $\alpha$ (from $2^{-30}$ to $2^{-100}$), checking if convergence was
reached afterwards.

The results of this method are summarized in Table 2, and show good agreement
between numerical and analytic results. Our integration program is written in
C using the GNU Scientific Library [Gou09] and Cuba. Its source code can be
found at https://hg.tx97.net/rambo, and also in the auxiliary files on arXiv.
With a requested accuracy of $0.1\%$ the complete integration takes less than
two days on a 12-core machine, with each integration taking between a minute
and two hours.

##  4 Conclusions

In this paper we present analytical expressions for five-particle phase-space
integrals expressed in terms of multiple zeta values up to weight 12. The
results are calculated using dimensional recurrence relation method with a
2000-digit accuracy using the DREAM package. We also present computer code for
the numerical integration of phase-space integrals in a higher-number of
dimensions that has been used to cross-check the obtained results with an
accuracy of $0.1\%$. The approach presented here shows excellent performance
for calculating single-scale integrals without ultra-violet divergences and
can be easily applied to other problems of this kind.

##  5 Acknowledgments

We are thankful to Sven Moch for numerous discussions and helpful suggestions
concerning this work, and for proofreading this paper. We were pleased to use
Axodraw2 [CV16] to draw diagrams for this paper.

This work was supported in part by the German Research Foundation DFG through
the Collaborative Research Centre No. SFB 676 Particles, Strings and the Early
Universe: the Structure of Matter and Space-Time.

##  Appendix A Results

The main results of our work are listed below. For brevity, we truncate them
up to MZVs of weight 6. Complete results with weight up to 12 are available in
auxiliary files on arXiv.

$\displaystyle F_{1}$
$\displaystyle=\frac{1}{72}+\epsilon\frac{53}{288}+\epsilon^{2}\Bigg{(}\frac{15961}{10368}-\frac{13}{72}\zeta_{2}\Bigg{)}+\epsilon^{3}\Bigg{(}\frac{436013}{41472}-\frac{689}{288}\zeta_{2}-\frac{13}{18}\zeta_{3}\Bigg{)}$
(19)
$\displaystyle+\epsilon^{4}\Bigg{(}\frac{96102601}{1492992}-\frac{207493}{10368}\zeta_{2}-\frac{689}{72}\zeta_{3}+\frac{17}{240}\zeta_{2}^{2}\Bigg{)}+\epsilon^{5}\Bigg{(}\frac{2206279853}{5971968}-\frac{5668169}{41472}\zeta_{2}$
$\displaystyle-\frac{207493}{2592}\zeta_{3}+\frac{901}{960}\zeta_{2}^{2}+\frac{169}{18}\zeta_{3}\zeta_{2}-\frac{65}{6}\zeta_{5}\Bigg{)}+\epsilon^{6}\Bigg{(}\frac{437728233961}{214990848}-\frac{1249333813}{1492992}\zeta_{2}$
$\displaystyle-\frac{5668169}{10368}\zeta_{3}+\frac{271337}{34560}\zeta_{2}^{2}+\frac{8957}{72}\zeta_{3}\zeta_{2}-\frac{3445}{24}\zeta_{5}-\frac{4007}{5040}\zeta_{2}^{3}+\frac{169}{9}\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{2}$
$\displaystyle=\frac{10}{3\epsilon^{5}}-\frac{55}{3\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{160}{3}-\frac{130}{3}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{490}{3}+\frac{715}{3}\zeta_{2}-\frac{512}{3}\zeta_{3}\Bigg{)}$
(20)
$\displaystyle+\frac{1}{\epsilon}\Bigg{(}\frac{1450}{3}-\frac{2080}{3}\zeta_{2}+\frac{2816}{3}\zeta_{3}+25\zeta_{2}^{2}\Bigg{)}-\frac{4390}{3}+\frac{6370}{3}\zeta_{2}-\frac{8192}{3}\zeta_{3}-\frac{275}{2}\zeta_{2}^{2}$
$\displaystyle+\frac{6656}{3}\zeta_{3}\zeta_{2}-2504\zeta_{5}+\epsilon\Bigg{(}\frac{13090}{3}-\frac{18850}{3}\zeta_{2}+\frac{25088}{3}\zeta_{3}+400\zeta_{2}^{2}-\frac{36608}{3}\zeta_{3}\zeta_{2}$
$\displaystyle+13772\zeta_{5}-211\zeta_{2}^{3}+\frac{13136}{3}\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{3}$
$\displaystyle=\frac{7}{8}-\frac{\zeta_{2}}{2}+\epsilon\Bigg{(}\frac{219}{16}-\frac{7}{2}\zeta_{2}-6\zeta_{3}\Bigg{)}+\epsilon^{2}\Bigg{(}\frac{4231}{32}-\frac{227}{8}\zeta_{2}-42\zeta_{3}-\frac{111}{10}\zeta_{2}^{2}\Bigg{)}$
(21)
$\displaystyle+\epsilon^{3}\Bigg{(}\frac{65347}{64}-\frac{3999}{16}\zeta_{2}-\frac{499}{2}\zeta_{3}-\frac{777}{10}\zeta_{2}^{2}+86\zeta_{3}\zeta_{2}-237\zeta_{5}\Bigg{)}+\epsilon^{4}\Bigg{(}\frac{887695}{128}$
$\displaystyle-\frac{64219}{32}\zeta_{2}-\frac{6303}{4}\zeta_{3}-\frac{5967}{16}\zeta_{2}^{2}+602\zeta_{3}\zeta_{2}-1659\zeta_{5}-\frac{1827}{20}\zeta_{2}^{3}+204\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{4}$
$\displaystyle=-\frac{\zeta_{2}}{\epsilon^{3}}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{\zeta_{2}}{2}-13\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{3}{2}\zeta_{2}+\frac{13}{2}\zeta_{3}-\frac{267}{10}\zeta_{2}^{2}\Bigg{)}+4\zeta_{2}+\frac{39}{2}\zeta_{3}+\frac{267}{20}\zeta_{2}^{2}$
(22)
$\displaystyle+191\zeta_{3}\zeta_{2}-\frac{1129}{2}\zeta_{5}+\epsilon\Bigg{(}10\zeta_{2}+52\zeta_{3}+\frac{801}{20}\zeta_{2}^{2}-\frac{191}{2}\zeta_{3}\zeta_{2}+\frac{1129}{4}\zeta_{5}-\frac{16547}{70}\zeta_{2}^{3}$
$\displaystyle+487\zeta_{3}^{2}\Bigg{)}$ $\displaystyle F_{5}$
$\displaystyle=\frac{7\zeta_{2}^{2}}{5\epsilon^{2}}+\frac{1}{\epsilon}\Bigg{(}-\frac{7}{2}\zeta_{2}^{2}-18\zeta_{3}\zeta_{2}+87\zeta_{5}\Bigg{)}-\frac{7}{10}\zeta_{2}^{2}+45\zeta_{3}\zeta_{2}-\frac{435}{2}\zeta_{5}+\frac{647}{5}\zeta_{2}^{3}-108\zeta_{3}^{2}$
(23) $\displaystyle F_{6}$
$\displaystyle=-\frac{35}{9\epsilon^{5}}-\frac{19}{6\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{278}{3}-\frac{191}{3}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{1697}{3}+\frac{1687}{18}\zeta_{2}-\frac{3236}{9}\zeta_{3}\Bigg{)}$
(24)
$\displaystyle+\frac{1}{\epsilon}\Bigg{(}-\frac{7793}{3}+\frac{10450}{9}\zeta_{2}+\frac{2386}{3}\zeta_{3}-\frac{4927}{18}\zeta_{2}^{2}\Bigg{)}+\frac{31259}{3}-\frac{67705}{9}\zeta_{2}+4296\zeta_{3}$
$\displaystyle+\frac{209681}{180}\zeta_{2}^{2}+\frac{14708}{3}\zeta_{3}\zeta_{2}-\frac{28148}{3}\zeta_{5}+\epsilon\Bigg{(}-\frac{117869}{3}+\frac{317353}{9}\zeta_{2}-31444\zeta_{3}$
$\displaystyle-\frac{66809}{45}\zeta_{2}^{2}-\frac{101858}{9}\zeta_{3}\zeta_{2}+\frac{84602}{3}\zeta_{5}-\frac{357871}{126}\zeta_{2}^{3}+\frac{96472}{9}\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{7}$
$\displaystyle=-1+\zeta_{2}+\epsilon\Bigg{(}-17+10\zeta_{2}+9\zeta_{3}\Bigg{)}+\epsilon^{2}\Bigg{(}-\frac{351}{2}+\frac{163}{2}\zeta_{2}+90\zeta_{3}+\frac{36}{5}\zeta_{2}^{2}\Bigg{)}$
(25)
$\displaystyle+\epsilon^{3}\Bigg{(}-\frac{5709}{4}+\frac{2495}{4}\zeta_{2}+\frac{1337}{2}\zeta_{3}+72\zeta_{2}^{2}-151\zeta_{3}\zeta_{2}+207\zeta_{5}\Bigg{)}+\epsilon^{4}\Bigg{(}-\frac{80649}{8}$
$\displaystyle+\frac{35823}{8}\zeta_{2}+\frac{18035}{4}\zeta_{3}+\frac{4881}{10}\zeta_{2}^{2}-1510\zeta_{3}\zeta_{2}+2070\zeta_{5}-\frac{387}{10}\zeta_{2}^{3}-387\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{8}$
$\displaystyle=\frac{1}{2\,\epsilon^{5}}-\frac{11}{4\,\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}8-\frac{49}{6}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{49}{2}+\frac{539}{12}\zeta_{2}-\frac{127}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{145}{2}-\frac{392}{3}\zeta_{2}$
(26)
$\displaystyle+\frac{1397}{6}\zeta_{3}-\frac{823}{60}\zeta_{2}^{2}\Bigg{)}-\frac{439}{2}+\frac{2401}{6}\zeta_{2}-\frac{2032}{3}\zeta_{3}+\frac{9053}{120}\zeta_{2}^{2}+\frac{1829}{3}\zeta_{3}\zeta_{2}$
$\displaystyle-\frac{2381}{3}\zeta_{5}+\epsilon\Bigg{(}\frac{1309}{2}-\frac{7105}{6}\zeta_{2}+\frac{6223}{3}\zeta_{3}-\frac{3292}{15}\zeta_{2}^{2}-\frac{20119}{6}\zeta_{3}\zeta_{2}+\frac{26191}{6}\zeta_{5}$
$\displaystyle+\frac{48983}{1260}\zeta_{2}^{3}+1397\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{9}$
$\displaystyle=-\frac{\zeta_{2}}{\epsilon^{3}}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{\zeta_{2}}{2}-7\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{3}{2}\zeta_{2}+\frac{7}{2}\zeta_{3}+\frac{17}{10}\zeta_{2}^{2}\Bigg{)}+4\zeta_{2}+\frac{21}{2}\zeta_{3}-\frac{17}{20}\zeta_{2}^{2}+131\zeta_{3}\zeta_{2}$
(27)
$\displaystyle-\frac{127}{2}\zeta_{5}+\epsilon\Bigg{(}10\zeta_{2}+28\zeta_{3}-\frac{51}{20}\zeta_{2}^{2}-\frac{131}{2}\zeta_{3}\zeta_{2}+\frac{127}{4}\zeta_{5}+\frac{7677}{70}\zeta_{2}^{3}+328\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{10}$
$\displaystyle=\zeta_{2}-\frac{3}{2}+\epsilon\Bigg{(}-\frac{47}{2}+\frac{15}{2}\zeta_{2}+11\zeta_{3}\Bigg{)}+\epsilon^{2}\Bigg{(}-\frac{903}{4}+58\zeta_{2}+\frac{165}{2}\zeta_{3}+\frac{171}{10}\zeta_{2}^{2}\Bigg{)}$
(28)
$\displaystyle+\epsilon^{3}\Bigg{(}-\frac{13795}{8}+\frac{951}{2}\zeta_{2}+\frac{1003}{2}\zeta_{3}+\frac{513}{4}\zeta_{2}^{2}-159\zeta_{3}\zeta_{2}+\frac{739}{2}\zeta_{5}\Bigg{)}+\epsilon^{4}\Bigg{(}-\frac{184655}{16}$
$\displaystyle+\frac{14539}{4}\zeta_{2}+3092\zeta_{3}+\frac{6507}{10}\zeta_{2}^{2}-\frac{2385}{2}\zeta_{3}\zeta_{2}+\frac{11085}{4}\zeta_{5}+\frac{1293}{14}\zeta_{2}^{3}-374\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{11}$
$\displaystyle=10\zeta_{2}\zeta_{3}-16\zeta_{5}+\epsilon\Bigg{(}-5\zeta_{3}\zeta_{2}+8\zeta_{5}-\frac{48}{35}\zeta_{2}^{3}+43\zeta_{3}^{2}\Bigg{)}$
(29) $\displaystyle F_{12}$
$\displaystyle=-\frac{2\zeta_{3}}{\epsilon}-13\zeta_{3}-\frac{37}{5}\zeta_{2}^{2}+\epsilon\Bigg{(}-76\zeta_{3}-\frac{481}{10}\zeta_{2}^{2}+30\zeta_{3}\zeta_{2}-115\zeta_{5}\Bigg{)}+\epsilon^{2}\Bigg{(}-422\zeta_{3}$
(30)
$\displaystyle-\frac{1406}{5}\zeta_{2}^{2}+195\zeta_{3}\zeta_{2}-\frac{1495}{2}\zeta_{5}-\frac{158}{7}\zeta_{2}^{3}+112\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{13}$
$\displaystyle=\frac{1}{6\epsilon^{5}}-\frac{11}{12\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{8}{3}+\frac{\zeta_{2}}{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{49}{6}-\frac{11}{4}\zeta_{2}+\frac{86}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{145}{6}+8\zeta_{2}-\frac{473}{3}\zeta_{3}$
(31)
$\displaystyle+\frac{5251}{60}\zeta_{2}^{2}\Bigg{)}-\frac{439}{6}-\frac{49}{2}\zeta_{2}+\frac{1376}{3}\zeta_{3}-\frac{57761}{120}\zeta_{2}^{2}-466\zeta_{3}\zeta_{2}+1794\zeta_{5}+\epsilon\Bigg{(}\frac{1309}{6}$
$\displaystyle+\frac{145}{2}\zeta_{2}-\frac{4214}{3}\zeta_{3}+\frac{21004}{15}\zeta_{2}^{2}+2563\zeta_{3}\zeta_{2}+-9867\zeta_{5}+\frac{1349771}{1260}\zeta_{2}^{3}-1372\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{14}$
$\displaystyle=-\frac{8\zeta_{2}^{2}}{5\epsilon}-\frac{28}{5}\zeta_{2}^{2}-48\zeta_{5}+\epsilon\Bigg{(}-\frac{164}{5}\zeta_{2}^{2}-168\zeta_{5}-\frac{6248}{105}\zeta_{2}^{3}+6\zeta_{3}^{2}\Bigg{)}$
(32) $\displaystyle F_{15}$
$\displaystyle=\frac{1}{6\epsilon^{6}}+\frac{41}{36\,\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{311}{36}-\frac{73}{18}\zeta_{2}\Bigg{)}+\Bigg{(}\frac{445}{18}-\frac{563}{36}\zeta_{2}-\frac{281}{9}\zeta_{3}\Bigg{)}\frac{1}{\epsilon^{3}}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{689}{9}$
(33)
$\displaystyle+\frac{5273}{36}\zeta_{2}-\frac{907}{18}\zeta_{3}-\frac{7103}{180}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{2024}{9}-\frac{7759}{18}\zeta_{2}+\frac{13933}{18}\zeta_{3}+\frac{10553}{120}\zeta_{2}^{2}$
$\displaystyle+\frac{1489}{3}\zeta_{3}\zeta_{2}-\frac{3257}{3}\zeta_{5}\Bigg{)}-\frac{6158}{9}+\frac{12437}{9}\zeta_{2}-\frac{22193}{9}\zeta_{3}+\frac{28621}{120}\zeta_{2}^{2}+\frac{14065}{18}\zeta_{3}\zeta_{2}$
$\displaystyle+\frac{3631}{6}\zeta_{5}-\frac{134489}{420}\zeta_{2}^{3}+1189\zeta_{3}^{2}$
$\displaystyle F_{16}$
$\displaystyle=\frac{1}{6\epsilon^{6}}+\frac{7}{12\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{185}{36}-\frac{65}{18}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{209}{18}-\frac{157}{12}\zeta_{2}-\frac{289}{9}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{76}{3}$
(34)
$\displaystyle+\frac{3563}{36}\zeta_{2}-\frac{215}{2}\zeta_{3}-\frac{10927}{180}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{239}{9}-\frac{2123}{18}\zeta_{2}+\frac{14431}{18}\zeta_{3}-\frac{7041}{40}\zeta_{2}^{2}$
$\displaystyle+463\zeta_{3}\zeta_{2}-1411\zeta_{5}\Bigg{)}+\frac{857}{9}-\frac{688}{3}\zeta_{2}-\frac{5623}{9}\zeta_{3}+\frac{170591}{120}\zeta_{2}^{2}+\frac{8963}{6}\zeta_{3}\zeta_{2}$
$\displaystyle-\frac{8349}{2}\zeta_{5}-\frac{134141}{180}\zeta_{2}^{3}+\frac{4169}{3}\zeta_{3}^{2}$
$\displaystyle F_{17}$
$\displaystyle=\frac{4\zeta_{3}}{\epsilon^{2}}+\frac{1}{\epsilon}\Bigg{(}10\zeta_{3}+\frac{84}{5}\zeta_{2}^{2}\Bigg{)}+48\zeta_{3}+42\zeta_{2}^{2}-52\zeta_{3}\zeta_{2}+280\zeta_{5}+\epsilon\Bigg{(}236\zeta_{3}+\frac{1008}{5}\zeta_{2}^{2}$
(35)
$\displaystyle-130\zeta_{3}\zeta_{2}+700\zeta_{5}+\frac{588}{5}\zeta_{2}^{3}-160\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{18}$
$\displaystyle=-\frac{5}{12\epsilon^{6}}+\frac{1}{8\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}\frac{205}{12}+\frac{323}{36}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{1937}{18}-\frac{71}{24}\zeta_{2}+\frac{223}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{4019}{9}$
(36)
$\displaystyle-\frac{3137}{12}\zeta_{2}-\frac{443}{18}\zeta_{3}+\frac{46931}{360}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}-\frac{3047}{2}+\frac{26681}{18}\zeta_{2}-1465\zeta_{3}-\frac{5473}{144}\zeta_{2}^{2}$
$\displaystyle-\frac{9113}{9}\zeta_{3}\zeta_{2}+\frac{27953}{9}\zeta_{5}\Bigg{)}+\frac{39460}{9}-\frac{54149}{9}\zeta_{2}+\frac{60682}{9}\zeta_{3}-\frac{33061}{24}\zeta_{2}^{2}$
$\displaystyle+\frac{6415}{18}\zeta_{3}\zeta_{2}-\frac{15967}{18}\zeta_{5}+\frac{12313867}{7560}\zeta_{2}^{3}-2302\zeta_{3}^{2}$
$\displaystyle F_{19}$
$\displaystyle=-\frac{1}{4\epsilon^{6}}+\frac{19}{8\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}\frac{64}{9}+\frac{173}{36}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{1895}{18}-\frac{2767}{72}\zeta_{2}+\frac{371}{9}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{9157}{18}$
(37)
$\displaystyle-\frac{1451}{18}\zeta_{2}-\frac{3439}{18}\zeta_{3}+\frac{32357}{360}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}-\frac{5798}{3}+\frac{24635}{18}\zeta_{2}-\frac{7103}{18}\zeta_{3}-\frac{17023}{720}\zeta_{2}^{2}$
$\displaystyle-\frac{4127}{9}\zeta_{3}\zeta_{2}+\frac{18073}{9}\zeta_{5}\Bigg{)}+\frac{111845}{18}-\frac{120679}{18}\zeta_{2}+\frac{49330}{9}\zeta_{3}-\frac{3664}{5}\zeta_{2}^{2}+\frac{50743}{18}\zeta_{3}\zeta_{2}$
$\displaystyle-\frac{53357}{18}\zeta_{5}+\frac{10750309}{7560}\zeta_{2}^{3}-\frac{7078}{9}\zeta_{3}^{2}$
$\displaystyle F_{20}$
$\displaystyle=\frac{20}{9\epsilon^{6}}+\frac{85}{6\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{775}{12}-\frac{325}{9}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{985}{18}-\frac{1705}{9}\zeta_{2}-\frac{1595}{9}\zeta_{3}\Bigg{)}$
(38)
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{6445}{9}-\frac{88}{9}\zeta_{2}^{2}+\frac{3815}{4}\zeta_{2}-\frac{4775}{6}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}-3635+\frac{8935}{18}\zeta_{2}+\frac{26455}{6}\zeta_{3}$
$\displaystyle-\frac{743}{36}\zeta_{2}^{2}+\frac{23485}{9}\zeta_{3}\zeta_{2}-\frac{7535}{3}\zeta_{5}\Bigg{)}+\frac{121310}{9}-\frac{80005}{9}\zeta_{2}+\frac{7415}{9}\zeta_{3}+\frac{9299}{24}\zeta_{2}^{2}$
$\displaystyle+\frac{189155}{18}\zeta_{3}\zeta_{2}-\frac{26065}{2}\zeta_{5}+\frac{134581}{126}\zeta_{2}^{3}+\frac{54985}{9}\zeta_{3}^{2}$
$\displaystyle F_{21}$
$\displaystyle=\frac{1}{2\epsilon^{5}}-\frac{11}{4\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}8+\frac{3}{2}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{49}{2}-\frac{33}{4}\zeta_{2}+\frac{158}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{145}{2}+24\zeta_{2}$
(39)
$\displaystyle-\frac{869}{3}\zeta_{3}+\frac{1205}{12}\zeta_{2}^{2}\Bigg{)}-\frac{439}{2}-\frac{147}{2}\zeta_{2}+\frac{2528}{3}\zeta_{3}-\frac{13255}{24}\zeta_{2}^{2}-\frac{2326}{3}\zeta_{3}\zeta_{2}+\frac{5614}{3}\zeta_{5}$
$\displaystyle+\epsilon\Bigg{(}\frac{1309}{2}+\frac{435}{2}\zeta_{2}-\frac{7742}{3}\zeta_{3}+\frac{4820}{3}\zeta_{2}^{2}+\frac{12793}{3}\zeta_{3}\zeta_{2}-\frac{30877}{3}\zeta_{5}+\frac{565729}{1260}\zeta_{2}^{3}$
$\displaystyle-\frac{5108}{3}\zeta_{3}^{2}\Bigg{)}$ $\displaystyle F_{22}$
$\displaystyle=\frac{23}{9\epsilon^{6}}+\frac{55}{9\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{2521}{36}-\frac{127}{3}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{1777}{9}-\frac{190}{3}\zeta_{2}-\frac{692}{3}\zeta_{3}\Bigg{)}$
(40)
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{11117}{18}+\frac{11419}{12}\zeta_{2}-\frac{400}{3}\zeta_{3}-\frac{11507}{90}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{32441}{18}-\frac{7993}{3}\zeta_{2}$
$\displaystyle+\frac{12061}{3}\zeta_{3}+\frac{6109}{18}\zeta_{2}^{2}+\frac{28472}{9}\zeta_{3}\zeta_{2}-\frac{15326}{3}\zeta_{5}\Bigg{)}-\frac{99143}{18}+\frac{50243}{6}\zeta_{2}-\frac{33388}{3}\zeta_{3}$
$\displaystyle-\frac{16091}{360}\zeta_{2}^{2}+\frac{13690}{9}\zeta_{3}\zeta_{2}+\frac{5165}{3}\zeta_{5}-\frac{539759}{630}\zeta_{2}^{3}+\frac{61504}{9}\zeta_{3}^{2}$
$\displaystyle F_{23}$
$\displaystyle=\frac{115}{18\epsilon^{6}}-\frac{115}{36\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{2645}{36}-\frac{1345}{18}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{3565}{18}+\frac{1345}{36}\zeta_{2}-\frac{2255}{9}\zeta_{3}\Bigg{)}$
(41)
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{5750}{9}+\frac{30935}{36}\zeta_{2}+\frac{2255}{18}\zeta_{3}+\frac{1235}{12}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{16445}{9}-\frac{41695}{18}\zeta_{2}$
$\displaystyle+\frac{51865}{18}\zeta_{3}-\frac{1235}{24}\zeta_{2}^{2}+\frac{26705}{9}\zeta_{3}\zeta_{2}-\frac{9910}{3}\zeta_{5}\Bigg{)}-\frac{50945}{9}+\frac{67250}{9}\zeta_{2}-\frac{69905}{9}\zeta_{3}$
$\displaystyle-\frac{28405}{24}\zeta_{2}^{2}-\frac{26705}{18}\zeta_{3}\zeta_{2}+\frac{4955}{3}\zeta_{5}-\frac{280151}{252}\zeta_{2}^{3}+\frac{45385}{9}\zeta_{3}^{2}$
$\displaystyle F_{24}$
$\displaystyle=\frac{1}{3\epsilon^{6}}+\frac{17}{9\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{179}{12}-\frac{67}{9}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{370}{9}-\frac{89}{3}\zeta_{2}-\frac{568}{9}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{2197}{18}$
(42)
$\displaystyle+\frac{9253}{36}\zeta_{2}-\frac{466}{3}\zeta_{3}-\frac{9659}{90}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{669}{2}-\frac{5840}{9}\zeta_{2}+\frac{14671}{9}\zeta_{3}-\frac{6923}{90}\zeta_{2}^{2}$
$\displaystyle+\frac{2836}{3}\zeta_{3}\zeta_{2}-2610\zeta_{5}\Bigg{)}-\frac{16739}{18}+\frac{3441}{2}\zeta_{2}-\frac{33956}{9}\zeta_{3}+\frac{628781}{360}\zeta_{2}^{2}+\frac{18284}{9}\zeta_{3}\zeta_{2}$
$\displaystyle-3109\zeta_{5}-\frac{746701}{630}\zeta_{2}^{3}+\frac{7664}{3}\zeta_{3}^{2}$
$\displaystyle F_{25}$
$\displaystyle=\frac{14}{9\epsilon^{5}}-\frac{1}{\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{358}{9}-\frac{203}{9}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{2182}{9}+\frac{187}{6}\zeta_{2}-\frac{1159}{9}\zeta_{3}\Bigg{)}$
(43)
$\displaystyle+\frac{1}{\epsilon}\Bigg{(}-\frac{10072}{9}+\frac{3946}{9}\zeta_{2}+\frac{5333}{18}\zeta_{3}-\frac{6451}{45}\zeta_{2}^{2}\Bigg{)}+\frac{40714}{9}-\frac{25759}{9}\zeta_{2}+1442\zeta_{3}$
$\displaystyle+\frac{48637}{90}\zeta_{2}^{2}+\frac{13613}{9}\zeta_{3}\zeta_{2}-3923\zeta_{5}+\epsilon\Bigg{(}-\frac{154768}{9}+\frac{122299}{9}\zeta_{2}-10793\zeta_{3}$
$\displaystyle-\frac{19313}{45}\zeta_{2}^{2}-\frac{21757}{6}\zeta_{3}\zeta_{2}+\frac{22801}{2}\zeta_{5}-\frac{3443939}{1890}\zeta_{2}^{3}+\frac{25937}{9}\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{26}$
$\displaystyle=2-\zeta_{2}+\epsilon\Bigg{(}32-9\zeta_{2}-10\zeta_{3}\Bigg{)}+\epsilon^{2}\Bigg{(}313-\frac{165}{2}\zeta_{2}-90\zeta_{3}-\frac{111}{10}\zeta_{2}^{2}\Bigg{)}+\epsilon^{3}\Bigg{(}\frac{4855}{2}$
(44)
$\displaystyle-\frac{2905}{4}\zeta_{2}-669\zeta_{3}-\frac{999}{10}\zeta_{2}^{2}+164\zeta_{3}\zeta_{2}-268\zeta_{5}\Bigg{)}+\epsilon^{4}\Bigg{(}\frac{65815}{4}-\frac{45421}{8}\zeta_{2}$
$\displaystyle-\frac{9533}{2}\zeta_{3}-\frac{12339}{20}\zeta_{2}^{2}+1476\zeta_{3}\zeta_{2}-2412\zeta_{5}+\frac{1383}{70}\zeta_{2}^{3}+430\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{27}$
$\displaystyle=\frac{2\zeta_{3}}{\epsilon^{2}}+\frac{1}{\epsilon}\Bigg{(}-\zeta_{3}+\frac{52}{5}\zeta_{2}^{2}\Bigg{)}-3\zeta_{3}-\frac{26}{5}\zeta_{2}^{2}-8\zeta_{3}\zeta_{2}+178\zeta_{5}+\epsilon\Bigg{(}-8\zeta_{3}-\frac{78}{5}\zeta_{2}^{2}$
(45)
$\displaystyle+4\zeta_{3}\zeta_{2}-89\zeta_{5}+\frac{5342}{35}\zeta_{2}^{3}+10\zeta_{3}^{2}\Bigg{)}$
$\displaystyle F_{28}$ $\displaystyle=-\frac{\zeta_{2}^{3}}{5}+2\zeta_{3}^{2}$
(46) $\displaystyle F_{29}$
$\displaystyle=\frac{5}{9\epsilon^{6}}-\frac{25}{18\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{5}{18}-\frac{55}{9}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{5}{9}+\frac{371}{18}\zeta_{2}-\frac{10}{9}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{10}{9}$
(47)
$\displaystyle-\frac{425}{18}\zeta_{2}+\frac{949}{9}\zeta_{3}+\frac{9167}{90}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}-\frac{20}{9}+\frac{595}{9}\zeta_{2}-\frac{4597}{9}\zeta_{3}+\frac{11021}{180}\zeta_{2}^{2}+\frac{578}{9}\zeta_{3}\zeta_{2}$
$\displaystyle+\frac{4357}{3}\zeta_{5}\Bigg{)}-\frac{40}{9}-\frac{946}{9}\zeta_{2}+\frac{10234}{9}\zeta_{3}-\frac{292403}{180}\zeta_{2}^{2}-\frac{10673}{9}\zeta_{3}\zeta_{2}+\frac{4225}{2}\zeta_{5}$
$\displaystyle+\frac{797803}{630}\zeta_{2}^{3}+\frac{430}{9}\zeta_{3}^{2}$
$\displaystyle F_{30}$
$\displaystyle=\frac{1}{6\epsilon^{6}}+\frac{13}{12\epsilon^{5}}+\frac{1}{\epsilon^{4}}\Bigg{(}-\frac{25}{3}-\frac{25}{6}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{3}}\Bigg{(}\frac{143}{6}-\frac{181}{12}\zeta_{2}-\frac{98}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}-\frac{443}{6}$
(48)
$\displaystyle+\frac{431}{3}\zeta_{2}-\frac{169}{3}\zeta_{3}-\frac{163}{4}\zeta_{2}^{2}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}\frac{1301}{6}-\frac{833}{2}\zeta_{2}+802\zeta_{3}+\frac{5887}{120}\zeta_{2}^{2}+\frac{1606}{3}\zeta_{3}\zeta_{2}$
$\displaystyle-\frac{3344}{3}\zeta_{5}\Bigg{)}-\frac{3959}{6}+\frac{7979}{6}\zeta_{2}-\frac{7286}{3}\zeta_{3}+\frac{5891}{15}\zeta_{2}^{2}+797\zeta_{3}\zeta_{2}+\frac{8}{3}\zeta_{5}-\frac{348377}{1260}\zeta_{2}^{3}$
$\displaystyle+\frac{4240}{3}\zeta_{3}^{2}$ $\displaystyle F_{31}$
$\displaystyle=\frac{7}{9\epsilon^{5}}-\frac{17}{18\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{(}-\frac{143}{9}-\frac{125}{9}\zeta_{2}\Bigg{)}+\frac{1}{\epsilon^{2}}\Bigg{(}\frac{902}{9}+\frac{133}{6}\zeta_{2}-\frac{236}{3}\zeta_{3}\Bigg{)}+\frac{1}{\epsilon}\Bigg{(}-\frac{4190}{9}$
(49)
$\displaystyle+\frac{716}{3}\zeta_{2}+\frac{1418}{9}\zeta_{3}-\frac{265}{6}\zeta_{2}^{2}\Bigg{)}+\frac{16892}{9}-\frac{4709}{3}\zeta_{2}+\frac{9718}{9}\zeta_{3}+\frac{3373}{20}\zeta_{2}^{2}+1228\zeta_{3}\zeta_{2}$
$\displaystyle-\frac{17612}{9}\zeta_{5}+\epsilon\Bigg{(}-\frac{63902}{9}+\frac{22181}{3}\zeta_{2}-\frac{68062}{9}\zeta_{3}-\frac{377}{5}\zeta_{2}^{2}-\frac{23666}{9}\zeta_{3}\zeta_{2}+\frac{48610}{9}\zeta_{5}$
$\displaystyle-\frac{688249}{1890}\zeta_{2}^{3}+\frac{27128}{9}\zeta_{3}^{2}\Bigg{)}$

## References

  * [AMV12] A. A. Almasy, S. Moch, and A. Vogt.  On the Next-to-Next-to-Leading Order Evolution of Flavour-Singlet Fragmentation Functions.  Nucl. Phys., B854:133–152, 2012.  arXiv:1107.2263, doi:10.1016/j.nuclphysb.2011.08.028.
  * [BBV10] J. Blumlein, D. J. Broadhurst, and J. A. M. Vermaseren.  The Multiple Zeta Value Data Mine.  Comput. Phys. Commun., 181:582–625, 2010.  arXiv:0907.2557, doi:10.1016/j.cpc.2009.11.007.
  * [BC10] P. A. Baikov and K. G. Chetyrkin.  Four Loop Massless Propagators: An Algebraic Evaluation of All Master Integrals.  Nucl. Phys., B837:186–220, 2010.  arXiv:1004.1153, doi:10.1016/j.nuclphysb.2010.05.004.
  * [CT81] K. G. Chetyrkin and F. V. Tkachov.  Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops.  Nucl. Phys., B192:159–204, 1981.  doi:10.1016/0550-3213(81)90199-1.
  * [CV16] J. C. Collins and J. A. M. Vermaseren.  Axodraw Version 2.  2016\.  arXiv:1606.01177.
  * [FBA99] H. Ferguson, D. Bailey, and S. Arno.  Analysis of PSLQ, an integer relation finding algorithm.  Mathematics of Computation of the American Mathematical Society, 68(225):351–369, 1999.  doi:10.1090/S0025-5718-99-00995-3.
  * [Fur03] H. Furusho.  The multiple zeta value algebra and the stable derivation algebra.  Publications of the Research Institute for Mathematical Sciences, 39(4):695–720, 2003.  arXiv:math/0011261, doi:10.2977/prims/1145476044.
  * [GGH04] A. Gehrmann-De Ridder, T. Gehrmann, and G. Heinrich.  Four particle phase space integrals in massless QCD.  Nucl.Phys., B682:265–288, 2004.  arXiv:hep-ph/0311276, doi:10.1016/j.nuclphysb.2004.01.023.
  * [Git16] O. Gituliar.  Master integrals for splitting functions from differential equations in QCD.  JHEP, 02:017, 2016.  arXiv:1512.02045, doi:10.1007/JHEP02(2016)017.
  * [GM15] O. Gituliar and S. Moch.  Towards three-loop QCD corrections to the time-like splitting functions.  Acta Phys. Polon., B46(7):1279–1289, 2015.  arXiv:1505.02901, doi:10.5506/APhysPolB.46.1279.
  * [Gou09] B. Gough.  GNU Scientific Library Reference Manual.  Network Theory Ltd., 3rd edition, 2009.
  * [Hah05] T. Hahn.  CUBA: A Library for multidimensional numerical integration.  Comput. Phys. Commun., 168:78–95, 2005.  arXiv:hep-ph/0404043, doi:10.1016/j.cpc.2005.01.010.
  * [KSE86] R. Kleiss, W. J. Stirling, and S. D. Ellis.  A New Monte Carlo Treatment of Multiparticle Phase Space at High-energies.  Comput. Phys. Commun., 40:359, 1986.  doi:10.1016/0010-4655(86)90119-0.
  * [Lap00] S. Laporta.  High precision calculation of multiloop Feynman integrals by difference equations.  Int. J. Mod. Phys., A15:5087–5159, 2000.  arXiv:hep-ph/0102033, doi:10.1016/S0217-751X(00)00215-7,10.1142/S0217751X00002157.
  * [Lap17] S. Laporta.  High-precision calculation of the 4-loop contribution to the electron g-2 in QED.  Phys. Lett., B772:232–238, 2017.  arXiv:1704.06996, doi:10.1016/j.physletb.2017.06.056.
  * [Lee14] R. N. Lee.  LiteRed 1.4: a powerful tool for reduction of multiloop integrals.  J. Phys. Conf. Ser., 523:012059, 2014.  arXiv:1310.1145, doi:10.1088/1742-6596/523/1/012059.
  * [Lep78] G. P. Lepage.  A New Algorithm for Adaptive Multidimensional Integration.  J. Comput. Phys., 27:192, 1978.  doi:10.1016/0021-9991(78)90004-9.
  * [LM16] R. N. Lee and K. T. Mingulov.  Introducing SummerTime: a package for high-precision computation of sums appearing in DRA method.  Comput. Phys. Commun., 203:255–267, 2016.  arXiv:1507.04256, doi:10.1016/j.cpc.2016.02.018.
  * [LM17a] R. N. Lee and K. T. Mingulov.  DREAM, a program for arbitrary-precision computation of dimensional recurrence relations solutions, and its applications.  2017\.  arXiv:1712.05173.
  * [LM17b] R. N. Lee and K. T. Mingulov.  Meromorphic solutions of recurrence relations and DRA method for multicomponent master integrals.  2017\.  arXiv:1712.05166.
  * [Smi15] A. V. Smirnov.  FIRE5: a C++ implementation of Feynman Integral REduction.  Comput. Phys. Commun., 189:182–191, 2015.  arXiv:1408.2372, doi:10.1016/j.cpc.2014.11.024.
  * [Tar96] O. V. Tarasov.  Connection between Feynman integrals having different values of the space-time dimension.  Phys. Rev., D54:6479–6490, 1996.  arXiv:hep-th/9606018, doi:10.1103/PhysRevD.54.6479.