# On the importance of data encoding in quantum Boltzmann methods

Merel A. Schalkers
Matthias Möller

(January 2023)

###### Abstract

In recent years, quantum Boltzmann methods have gained more and more interest
as they might provide a viable path towards solving fluid dynamics problems on
quantum computers once this emerging compute technology has matured and fault-
tolerant many-qubit systems become available. The major challenge in
developing a start-to-end quantum algorithm for the Boltzmann equation
consists in encoding relevant data efficiently in quantum bits (qubits) and
formulating the streaming, collision and reflection steps as one comprehensive
unitary operation. The current literature on quantum Boltzmann methods mostly
proposes data encodings and quantum primitives for individual phases of the
pipeline assuming that they can be combined to a full algorithm.

In this paper we disprove this assumption by showing that for encodings
commonly discussed in literature either the collision or the streaming step
cannot be unitary. Building on this landmark result we propose a novel
encoding in which the number of qubits used to encode the velocity depends on
the number of time steps one wishes to simulate, with the upper bound
depending on the total number of grid points.

In light of the non-unitarity result established for existing encodings, our
encoding method is to the best of our knowledge the only one currently known
that can be used for a start-to-end quantum Boltzmann solver where both the
collision and the streaming step are implemented as a unitary operation.
Furthermore our theoretical unitarity result can serve as a guideline on which
types of encodings to consider or whether a ‘stop-and-go’ method with repeated
measurements and re-initializations is the method of choice.

##  1 Introduction

Since the first quantum computing boom in the 1990s, quantum computational
fluid dynamics (QCFD) has been a field of interest to researchers worldwide.
Due to the high computational demands of classic CFD the exponential potential
of quantum computers in combination with quantum parallelism and quantum
indeterminacy has caused interest in the application. The first QCFD
algorithms were proposed by Yepez and his co-workers around the turn of the
century [Yep98, Yep01, YB01, Yep02, Pra+03]. These algorithms are based on a
quantum distributed computing approach assuming that many small-scale quantum
computers are more realistic than one large many-qubit system. The core idea
of the so-called quantum lattice-gas model is that each grid point of
position-space gets its own 6-qubit quantum computer associated to it (which
can also be groups of 6 qubits of a future many-qubit quantum computer). The
benefit of this approach is that the possible quantum circuit depth and stable
entanglement required remains very low, making it a realistic and relatively
near-term approach given the capabilities of current quantum devices. Its
downside is that to encode a grid of size $N$ a total of $6N$ qubits are
required, which means that the amount of qubits required grows linearly with
the size of the grid. Given the limited amount of quantum devices available
and the large amount of grid points required for solving practical problems
with modern Boltzmann methods, this distributed approach proves a significant
drawback. Furthermore, as we will show below, the computational basis state
encoding of the velocity vector adopted in the aforementioned papers does not
allow for implementing the streaming step as a unitary operator so that
measurement and state re-initialization is mandatory after each time step.

After these early results by Yepez et al., the QCFD field became stagnant for
about a decade until its recent resurgence, in particular, in the form of
quantum Boltzmann methods. Most recent are the methods presented in [TS20,
Bud20, Bud21, MVS22, SM22, Ste23], that all have their own strengths and
weaknesses. Some papers include a streaming and specular reflection mechanism,
but no collision methods yet [TS20, Bud20, SM22]. Other approaches have
implemented a collision method using the linear combination of unitary
approach [CW12], causing the algorithm to require a measurement-and-restart
strategy after each time step [Bud21]. Due to the high costs of quantum state
preparation and the chance of measurement errors this ‘stop-and-go’ strategy
is hardly usable in practice. Other algorithms have managed to create a
unitary collision operator, but have not yet been able to combine this with a
streaming step into one start-to-end algorithm [MVS22, Ste23].

What remained an open problem is the development of a full-fledged quantum
Boltzmann method (QBM) that implements both the streaming and the collision
step as unitary operations. In this paper we present the first-of-its-kind
full-fledged QBM building on a novel encoding scheme of the velocity vector
that scales with the number of time steps. Furthermore, we prove rigorously
that for the encoding schemes considered for universal quantum computers in
all previous publications it is impossible to implement both streaming and
collision as a unitary, downgrading them as candidates for any practical QBM.
Taking both contributions of this paper together, our new encoding and the
theoretical (negative) result on existing encodings, we hope to stimulate a
paradigm shift in QBM research from focusing on encodings and algorithms for
individual steps of the pipeline to developing full-fledged QBM algorithms.

##  2 Lattice Boltzmann method

In the Boltzmann method the macroscopic behavior of a fluid is simulated by
considering the microscopic behavior of the fluid particles as they move
through space and deriving the macroscopic quantities via averaging-based
post-processing, instead of encoding the macroscopic variables directly, as is
commonly done in other CFD methods like the finite volume method.

In this paper we consider the discrete lattice Boltzmann method, where a
particle can only move with specific velocities taken from a finite set of
discrete velocities. We define the structure of the method using the D$n$Q$m$
system, where $n$ represents the amount of spatial dimensions and $m$ the
amount of discrete velocities considered. Figure 1 gives examples of the
commonly used D1Q2, D1Q3, D2Q5 and D2Q9 systems, respectively, in standard
Boltzmann convention. For an in depth review of the lattice Boltzmann method
we refer to the book [Krü+17].

Figure 1: Four examples of different types of D$n$Q$m$ possible. The top
picture on the left portrays the D1Q2 setting and the top picture on the right
portrays the D1Q3 setting (where a stationary particle can be included). The
picture on the left bottom portrays the D2Q5 setting and the picture on the
right below shows the D2Q9 setting.

Boltzmann methods simulate the macroscopic behavior of a fluid or gas by
implementing a streaming step followed by particle collision on the
microscopic level in each time step. When obstacles are present an additional
reflection step is performed in each time step. For brevity we omit a detailed
description of the latter and refer the interested reader to our recent work
[SM22] on this topic.

The streaming step is implemented by letting the particles move by one grid
point per time step in the direction they are traveling currently. Figure 2
illustrates how the particles travel in one time step from the point $x$ to
$x\pm 1$ respectively for the D1Q2 case. Similar illustrations can be
constructed for two- and three-dimensional cases but are omitted here for
brevity reasons.

To implement the collision step we define so-called equivalence classes of
streaming patterns which have the same total mass and momentum and are thus
considered to be equivalent. A combination of colliding particles can
therefore be transformed into any combination from the same equivalence class
upon collision without changing the total mass and momentum. Figure 3 shows an
example of two equivalent velocity combinations for the D2Q5 (and D2Q4) case.

Figure 2: Illustration of the streaming step for the D1Q2 case. Figure (a)
shows the velocity vectors at position $x$ at time $t$. Figure (b) shows the
same after configuration at time $t+1$ after particles have moved to positions
$x-1$ and $x+1$, respectively. Red and blue colors identify the different
streaming directions and their propagation pattern.

In Section 4, we present a lattice Boltzmann encoding for which both the
collision and the streaming step can be performed through unitary operations
and thus admit a straightforward implementation on a sufficiently large fault-
tolerance quantum computer. Before that, in Section 3, we provide rigorous
mathematical proofs that show that such a unitary treatment of both streaming
_and_ collision is impossible with the encodings adopted in current
literature, thereby underpinning the uniqueness of our proposed encoding.

Figure 3: Illustration of two velocity combinations of the D2Q5 (and D2Q4)
velocity spectrum that belong to the same equivalence class with total
momentum 0 and mass 2: (a) particles streaming in the $q_{1}$ and $q_{3}$
direction, and (b) particles streaming in the $q_{2}$ and $q_{4}$ direction.

##  3 Data encoding

As in any computational field, data encoding is pivotal for reaching a good
result. More than five decades of classical CFD research and application have
established ‘good practices’ for storing field data such as densities and
velocities at, e.g., the grid points or cell centers as floating-point numbers
following the IEEE-754 standard. Every now and then new hardware developments
stimulate research into non-standard formats, like reduced or mixed-precision
[Fre+22], but, in general, data encoding is not considered to be an open
problem.

Not so in QCFD and, in particular, quantum Boltzmann methods. In this section
we will review the main data encodings currently used for QBM and show that in
all of them either the streaming step or the collision step cannot be unitary.
This result, though discouraging at first sight, should be interpreted as
wake-up call that novel quantum encodings for CFD states are imperative for
devising full-fledged QCFD applications in the future. We propose one such
novel encoding in Section 4 and discuss its potential and limitations.

The two mainstream encodings of the velocity vector are the amplitude based
encoding and the computational basis state encoding. In what follows, we will
consider both approaches separately and show how they both lead to a
contradiction in the unitarity of either the collision or the streaming
operation.

###  3.1 Amplitude based encoding

The first type of encoding we consider is the so-called amplitude based
encoding, used for several quantum Boltzmann methods [TS20, Bud20, Bud21,
SM22]. The amplitude based encoding of the velocity vector is such that at
each location $\ket{x}$ there can be multiple particles with different
velocities, for instance $\ket{v_{0}}$, $\ket{v_{1}}$, $\ket{v_{2}}$ and
$\ket{v_{3}}$ for D2Q4. Here and below, $\ket{i}$ denotes the representation
of $i$ as bit string. The state of the system at this point $x$ can then be
encoded as111Note that we distinguish in our notation between the grid point
$x$ and its representation $\ket{x}$ as part of the quantum register.

$\ket{x}\left(\alpha_{0}\ket{v_{0}}+\alpha_{1}\ket{v_{1}}+\alpha_{2}\ket{v_{2}}+\alpha_{3}\ket{v_{3}}\right),$
(1)

where $\alpha_{0}$, $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are complex
numbers that simply represent the relative weight or amount of particles
traveling at the given velocity at grid point $x$. For simplicity we will
assume that
$|\alpha_{0}|^{2}+|\alpha_{1}|^{2}+|\alpha_{2}|^{2}+|\alpha_{3}|^{2}=1$, and
so in this example there are only particles at grid point $x$ but the proof
extends trivially to the general case with particles spread around the grid.

In order to show that this encoding of the velocity vector inevitably leads to
non-unitary collision operators, let us assume a system in a specific state
$\ket{\psi_{1}}$, with only particles traveling with velocities $\ket{v_{0}}$
and $\ket{v_{1}}$, meaning that $|\alpha_{0}|,|\alpha_{1}|>0$ and
$\alpha_{2}=\alpha_{3}=0$. Then we can write the state of the system as

$\ket{\psi_{1}}=\ket{x}\left(\alpha_{0}\ket{v_{0}}+\alpha_{1}\ket{v_{1}}\right).$
(2)

Now assume that an equivalent velocity combination exists consisting of
particles traveling with velocities
$\beta_{2}\ket{v_{2}}+\beta_{3}\ket{v_{3}}$, where we have
$|\beta_{2}|,|\beta_{3}|>0$. To realize this potential outcome of a collision
as a quantum algorithm, we need to implement the transformation between both
equivalent states as a unitary operation $U_{\text{col}}$ which changes the
states of the velocity encodings as follows

$\begin{split}\ket{\psi_{1}^{\prime}}&=I\otimes
U_{\text{col}}\ket{\psi_{1}}\\\ &=\ket{x}\otimes
U_{\text{col}}\left(\alpha_{0}\ket{v_{0}}+\alpha_{1}\ket{v_{1}}\right)\\\
&=\ket{x}\left(\gamma_{0}(\alpha_{0}\ket{v_{0}}+\alpha_{1}\ket{v_{1}})+\gamma_{1}(\beta_{2}\ket{v_{2}}+\beta_{3}\ket{v_{3}})\right).\end{split}$
(3)

Here, if $\gamma_{0}=1$ and $\gamma_{1}=0$ no collision is taking place (and
we simply implement an identity operation) and if $\gamma_{1}=1$ we fully
change from the original velocities to its alternative representative from the
same equivalence class.222Here $\gamma_{0}$, $\gamma_{1}$ are chosen to
reflect the fact that a collision operation should switch weight of a
combination of velocities in an equivalent class to another combination of
velocities in the same equivalence class. This equation could be written in a
less restrictive way by splitting $\gamma_{0}$ and $\gamma_{1}$ up into
separate amplitudes $\gamma_{i}$ for all the basis states $\ket{v_{i}}$, the
same contradiction of unitarity as presented below however could be reached.
Note that to preserve unitarity $|\gamma_{0}|^{2}+|\gamma_{1}|^{2}=1$ must
hold.

Let us now consider another system in state
$\ket{\psi_{2}}=\ket{x}\ket{v_{2}}$. Applying the unitary operation
$U_{\text{col}}$ should not effect the state at all as a single speed is only
in an equivalence class with itself, and so the required behavior for
$U_{\text{col}}$ is

$\begin{split}\ket{\psi_{2}^{\prime}}&=I\otimes
U_{\text{col}}\ket{\psi_{2}}\\\ &=\ket{x}U_{\text{col}}\ket{v_{2}}\\\
&=e^{i\theta}\ket{x}\ket{v_{2}},\end{split}$ (4)

with $\theta\in(0,2\pi]$. That is, the collision operator must preserve the
single-velocity state except for changes in the phase factor $e^{i\theta}$
that can be neglected.

Now that we have identified the required behavior for $U_{\text{col}}$ to
implement a collision operation, we can prove that any $U_{\text{col}}$ that
meets both requirements simultaneously cannot be unitary. Here, we resort to
the characterization $U_{\text{col}}^{\dagger}U_{\text{col}}=I$ of unitary
operators, with superscript $\dagger$ denoting the adjoint operator.

###### Proof.

To reach a contradiction, assume that $U_{\text{col}}$ is a unitary operator.
Then it must preserve the inner product for all possible states
$\ket{\phi_{1}}$ and $\ket{\phi_{2}}$

$\braket{\phi_{1}}{\phi_{2}}=\bra{\phi_{1}}U_{\text{col}}^{\dagger}U_{\text{col}}\ket{\phi_{2}}.$
(5)

However, for a collision operation $U_{\text{col}}$ that behaves as expected
on the system states described in Equations (2) to (4), it follows that

$\begin{split}0&=\braket{\psi_{1}}{\psi_{2}}\\\ &=\bra{\psi_{1}}\left(I\otimes
U_{\text{col}}\right)^{\dagger}\left(I\otimes
U_{\text{col}}\right)\ket{\psi_{2}}\\\
&=e^{i\theta}\left(\gamma_{0}(\alpha_{0}\bra{v_{0}}+\alpha_{1}\bra{v_{1}})+\gamma_{1}(\beta_{2}\bra{v_{2}}+\beta_{3}\bra{v_{3}})\right)\bra{x}\ket{x}\ket{v_{2}}\\\
&=e^{i\theta}\gamma_{1}\beta_{2}.\end{split}$ (6)

The first equality follows from the fact that $\ket{\psi_{1}}$ and
$\ket{\psi_{2}}$ are orthogonal by construction. The second one holds under
the assumption of $U_{\text{col}}$ being unitary, which is disproved by the
fact that the entire equality chain only holds for the trivial case
$\gamma_{1}=0$ (as $|\beta_{2}|>0$ by definition of the state
$\ket{\psi_{1}}$), that is, when $U_{\text{col}}$ does not implement the
collision operation. From this we can conclude that an amplitude based
encoding of the velocity does not allow for a unitary implementation of the
collision operation. ∎

Notice that this proof works for any amplitude based encoding of $v$ where the
different possible velocities at a position are all represented by their own
basis state as there will always be a case with only a single incoming
velocity, for which an identity operation up to a phase shift should take
place, while at the same time there will be combinations of velocities for
which we want some weight of the system to change from one combination of
velocities to another combination of velocities in the same equivalence class.
These two antagonizing requirements will always lead to the same contradiction
of unitarity proven above and we further expand on this intuition in Section
3.3.

###  3.2 Computational basis state encoding

The second type of encoding of a quantum state considered is the computational
basis encoding, used in several quantum lattice Boltzmann papers such as
[Yep98, Yep01, YB01, Yep02, Pra+03, MVS22, Ste23]. Using this encoding the
contradiction of unitarity in the collision operation can be avoided by
encoding the velocity of the qubits at a position $\ket{x}$ in space by
identifying each direction particles could be streamed from with its own
qubit, which will be set to one if and only if there is a particle streaming
from that direction.

As an example consider the D2Q4 lattice depicted in Figure 4. In this case the
velocity can be encoded using four qubits $q_{0}$, $q_{1}$, $q_{2}$ and
$q_{3}$ where the state

$\ket{x}\ket{v}=\ket{x}\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{x}\ket{0110}$ (7)

is such that from the center point $(1,1)$, there is a particle streaming to
$(1,2)$ and a particle streaming to $(0,1)$ but not to $(2,1)$ or $(1,0)$.

Figure 4: Illustration of the computational basis state encoding for the D2Q4
lattice. For each grid point $x$ we set the respective qubit $q_{j}$ to one if
and only if there is a particle streaming in that direction, i.e.
$\ket{v}=\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{0110}$.

Using this encoding the collision step can be defined quite naturally as
unitary operation. However, we run into trouble when attempting to define a
unitary streaming step $U_{\text{str}}$ as we demonstrate in what follows.

To simplify notation let us restrict ourselves to the D1Q2 lattice and
consider the two settings at time $t$ from Figures 5 and 6, which can be
encoded as

$\displaystyle\ket{\psi_{1}}$ $\displaystyle=\sum_{x=0}^{3}\ket{x}\ket{v}$ (8)
$\displaystyle=\frac{1}{2}\left(\ket{00}\ket{00}+\ket{01}\ket{11}+\ket{10}\ket{10}+\ket{11}\ket{10}\right),$
(9)

and

$\ket{\psi_{2}}=\frac{1}{2}\left(\ket{00}\ket{01}+\ket{01}\ket{01}+\ket{10}\ket{00}+\ket{11}\ket{11}\right),$
(10)

respectively. It then follows directly that

$\braket{\psi_{1}}{\psi_{2}}=0.$ (11)

Upon streaming, the systems from Figures 5 and 6 change from their state at
time $t$ (top lattice) to that at time $t+1$ (bottom lattice), i.e.

$\ket{\psi_{1}^{\prime}}=\frac{1}{2}\left(\ket{00}\ket{11}+\ket{01}\ket{00}+\ket{10}\ket{10}+\ket{11}\ket{10}\right),$
(12)

and

$\ket{\psi_{2}^{\prime}}=\frac{1}{2}\left(\ket{00}\ket{11}+\ket{01}\ket{00}+\ket{10}\ket{01}+\ket{11}\ket{01}\right),$
(13)

respectively. As in the previous section, we will show by contradiction that
any operation $U_{\text{str}}$ for which
$U_{\text{str}}\ket{\psi_{1}}=\ket{\psi_{1}^{\prime}}$ and
$U_{\text{str}}\ket{\psi_{2}}=\ket{\psi_{2}^{\prime}}$ cannot be unitary.

###### Proof.

Let us assume that $U_{\text{str}}$ is unitary, i.e. it preserves the inner
product

$\braket{\phi_{1}}{\phi_{2}}=\bra{\phi_{1}}U^{\dagger}_{\text{str}}U_{\text{str}}\ket{\phi_{2}},$
(14)

for all states $\ket{\phi_{1}},\ket{\phi_{2}}$. Substituting the states (9)
and (10) on the left side, and (12) and (13) into the right inner product we
arrive at the contradiction

$0=\braket{\psi_{1}}{\psi_{2}}=\bra{\psi_{1}}U_{\text{str}}^{\dagger}U_{\text{str}}\ket{\psi_{2}}=\braket{\psi_{1}^{\prime}}{\psi_{2}^{\prime}}=\frac{1}{2}.$
(15)

The first equality follows from the orthogonality property (11), and the
second one from the assumption that $U_{\text{str}}$ is a unitary operator,
which we just disproved. ∎

Figure 5: D1Q2 example setting 1. The binary encoding above the arrows
indicate whether or not a particle is flowing there in that time step. 1
indicates that there is a particle streaming there and 0 indicates that there
is no particle. In the example setting we consider periodic boundary
conditions. The top figure shows the state of the system at time $t$. The
figure below shows the state of the system at time $t+1$.

Figure 6: D1Q2 example setting 2, the binary encoding above the arrows
indicate whether or not a particle is flowing there in that time step. 1
indicates that there is a particle streaming there and 0 indicates that there
is no particle. In the example setting we consider periodic boundary
conditions. The top figure shows the state of the system at time $t$. The
figure below shows the state of the system at time $t+1$.

As in Section 3.1 this proof extends to any computational basis encoding where
each possible combination of velocities at a specific lattice point is encoded
using its own basis state, as one can always construct two situations with no
overlap at time $t$ that will have non-zero overlap after streaming at time
$t+1$. This proof also extends trivially to any other D$n$Q$m$ setting as the
streaming possibilities of D1Q2 are essentially a subset of any other system
and thus the same example can be used by setting the other streaming
directions to 0.

###  3.3 Intuition and extension of non-unitarity proofs

In this section we expand on our non-unitarity proofs by providing physical
intuition behind the proofs presented above. It is intended to give insight
into what types of encodings our non-unitarity proof extends to, and what
physical features of the system necessarily lead to the non-unitarity for
these encodings.

Consider the proof from Section 3.1 that shows that the amplitude based
encoding, where each velocity direction is identified through its own basis
state leaving the total velocity at a position $x$ to be a superposition of
such basis states, prevents the collision operator $U_{\text{col}}$ from being
unitary. Since it encodes each streaming direction as a different basis state,
the quantum encodings of the velocity directions are all orthogonal to one
another. This is also necessary, since if the basis states of the possible
streaming directions are not orthogonal, we cannot fully distinguish between
them. However, this orthogonality of the different velocity directions leads
directly to the non-unitarity of $U_{\text{col}}$. Since a collision operator
that will rotate a given linear combination of basis states into a linear
combination of other basis states in such a way that the represented streaming
patterns belong to the same equivalence class, it will also rotate ‘pure’
velocities represented by a single basis state into another basis state
leading to a nonphysical and undesired change of velocities.

Following this line of argumentation it can be seen that the non-unitarity of
$U_{\text{col}}$ is not so much a result of a specific choice of encoding but
an inherent non-unitarity of the collision step itself that directly leads to
the idea of computational basis state encoding, where each velocity pattern
(i.e. the combination of velocities) at a grid point is encoded as its own
basis state, and not as a unitary combination of all the basis states
representing a non-zero contribution.

When encoding the velocity pattern at each grid point as a basis state,
naturally, the non-unitarity of collision falls away and we can find a
straightforward unitary operator to implement the collision step. However,
such an encoding will always lead to non-unitarity of streaming due to the
non-local nature of a streaming operation. Consider an arbitrary point in
space $x$ and imagine two different scenarios with two different combinations
of speeds $\ket{v_{1}}$ and $\ket{v_{2}}$ at this point. Then the inner
product between $\ket{x}\ket{v_{1}}$ and $\ket{x}\ket{v_{2}}$ must be 0, as
these are different basis states. However, the velocity states of the systems
at position $x$ in the next time step do not depend on the current velocity
states in the lattice point. In fact, they only depend on the velocity states
of the neighboring lattice points. Since the inner product of the states at
the point $x$ at the next time step does not depend on the current states at
the point $x$, in the next time step the velocity at the point $x$ of the two
systems could be identical, and hence, the inner product could be one. There
is no way of ensuring that this can only happen when the inner product at some
other point $x^{\prime}$ of the systems was non-zero before as each grid point
has velocity vectors in multiple directions determining its associated
velocity basis state.

This shows that any quantum encoding that successfully implements both
streaming and collision as a unitary operation must belong to one of the
following three types. The first type is an amplitude based type encoding,
where the different velocities are not orthogonal and thus not entirely
distinguishable. The second type is a computational basis state encoding where
the non-locality of streaming is somehow avoided. The last type is a
completely novel encoding method that avoids both non-unitarity problems
entirely. In the next section we will present precisely one such idea.

##  4 Space-time data encoding

In this section we propose a novel space-time data encoding that enables
unitary collision _and_ streaming at the same time. To the best of our
knowledge, this is the first-of-its-kind start-to-end quantum Boltzmann
algorithm that does not require measurement and quantum-state re-
initialization after each time step.

In what follows, we adopt an extended computational basis state encoding,
where at each location $x$ we take into account the velocities at all grid
points in the vicinity of $x$. Here, ‘in the vicinity of $x$’ means that a
particle can theoretically reach the grid point $x$ within the number of time
steps still to be performed before measurement. Mathematically speaking being
‘in the vicinity of $x$’ means being, respectively, in the so-called extended
von Neumann, Moore or hexagonal neighborhood of the point $x$, depending on
the lattice structure.333The von Neumann neighborhood of extent $r$ defines
the diamond-shaped set of points at a Manhattan distance of up to $r$ from the
point $x$. It applies to, e.g., D2Q4, D2Q5, D3Q6, and D3Q7. The Moore distance
extends the former one by diagonal directions and applies to, e.g., D2Q8,
D2Q9, D3Q26, and D3Q27. As its name suggests, the hexagonal neighborhood
applies to D2Q6 and its extension to its three-dimensional counterpart.

This leads to a trade-off between the number of time steps that can be
performed between measurements and the number of qubits required to encode the
velocity at each grid point $x$. The more time steps one wishes to take
between measurement-and-re-initialization cycles, the more qubits are required
for our space-time encoding. Obviously the maximum number of qubits required
to implement the velocity without any in-between measurements must be such
that the entire grid is spanned. For a D$n$Q$m$ lattice this will be $mN_{g}$,
where $N_{g}$ is the total number of grid points. When encoding the proposed
method on a classical computer $mN_{g}$ bits would also be required, so when
encoding the full domain there is no quantum benefit in terms of (qu)bit
numbers. The quantum improvement comes from exploiting quantum parallelism,
which is done as long as we do not encode the whole space.

In what follows, let $N_{t}$ denote the number of streaming steps to be
performed between (re-)initialization and measurement. We extend the
computational basis state encoding of velocity directions from Section 3.2 to
take into account all the speed states from grid points in the neighborhood of
$x$ that can (at least theoretically) reach $x$ within $N_{t}$ streaming
steps. This takes away the non-locality of the streaming operator, which led
to the non-unitarity of $U_{\text{str}}$ for the ‘regular’ computational basis
state encoding at the cost of increasing the number of qubits required to
encode all required velocity data.

We will give a detailed description of this encoding for the D2Q4 lattice, but
want to note that it can be extended naturally to any other choice of
D$n$Q$m$. Consider the D2Q4 lattice given in Figure 7 with qubit $q_{j}$ set
to one if and only if there is a particle traveling with velocity direction
$j$ from grid point $x$ into a neighboring grid point in the current time
step. We now extend this encoding to include _all_ possible velocities at
positions ‘in the vicinity of $x$’ for the total of $N_{t}$ time steps in
order to obtain a unitarily streamable encoding. This is illustrated in Figure
8 for a single time step, i.e. $N_{t}=1$ yielding the encoding

$\ket{x}\ket{q_{19}q_{18}\dots q_{0}}.$ (16)

For D2Q4, the number of qubits encoding the possible velocity states per grid
location $x$ grows with the number of time steps (still) to be taken as

$n_{v}=4+\sum_{i=1}^{N_{t}}16i=8N_{t}^{2}+8N_{t}+4,$ (17)

where the maximum number of qubits required to encode all velocity directions
over the entire grid equals $4N_{g}$ as stated before.444Note that the growth
rate of qubit numbers per time step depends on the choice of D$n$Q$m$. The
number of qubits required is equal to the number of points in the extended Von
Neumann, Moore or hexagonal neighborhood, depending on which choice of $n$ and
$m$ considered. Similarly it can be shown that for $d$ dimensions the growth
rate is of the order $\mathcal{O}\left(N_{t}^{d}\right)$.

We can now encode the collision step by first identifying the equivalence
class for the D2Q4 lattice. We note that at each grid point $x$ as represented
in Figure 7 the states $\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{1010}$ and
$\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{0101}$ belong to the same equivalence class
(cf. Figure 3), as they have the same total mass and momentum.555The other
equivalence classes are $\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{1000}$ and
$\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{1100}$ and all cyclic shifts of these
patterns, and $\ket{q_{0}q_{1}q_{2}q_{3}}=\ket{1111}$. However, they all have
just a single representative so that we define the collision operator based on
the ambiguous case. We implement the collision step by defining a unitary
operator $U_{\text{col}}$ which performs the following mappings

$\displaystyle U_{\text{col}}\ket{1010}$
$\displaystyle=\phantom{-}\alpha\ket{1010}+\beta\ket{0101},$ (18)
$\displaystyle U_{\text{col}}\ket{0101}$
$\displaystyle=-\beta\ket{1010}+\alpha\ket{0101},$ (19)

with $\alpha,\beta\in\mathbb{C}$ and $|\alpha|^{2}+|\beta|^{2}=1$, while
acting as the identity operation on any other basis state. It can easily be
verified that this operation is unitary. With the so-defined
$U_{\text{col}}\in\mathbb{C}^{2^{4}\otimes 2^{4}}$, we can write the total
collision operation for an encoding of the velocities states $v$ consisting of
$n_{v}=4k$ qubits as $k$-fold Kronecker products of $U_{\text{col}}$
operations, i.e.
$U_{\text{col}}^{\text{tot}}=U_{\text{col}}\otimes\dots\otimes
U_{\text{col}}$. Since each $U_{\text{col}}$ requires a few CNOT and a single
triple controlled rotation gate the total collision operator can be
efficiently implemented even on near-term devices.666We can implement the
described collision operator by first applying three CNOT operations to the
system turning the states into $\ket{1010}\mapsto\ket{1110}$ and
$\ket{0101}\mapsto\ket{1111}$. Subsequently a triple controlled rotation
operation of choice is applied to the right-most qubit (controlled on the
three left-most qubits). Finally the initial three CNOT operations are applied
in reverse order to reset all velocity states correctly.

In practice the total collision operator $U_{\text{col}}^{\text{tot}}$ differs
per time step, since its local counterpart $U_{\text{col}}$ only needs to be
applied to velocity states ‘in the vicinity of $x$’. In the first out of the
$N_{t}$ time steps it is important for all qubits representing velocity states
‘in the vicinity of $x$’ to be updated correctly. In the very last time step,
however, it is only important for the qubits $q_{0}$, $q_{1}$, $q_{2}$ and
$q_{3}$ to end up in the correct state. The more time steps $t$ have been
taken, the less time steps $N_{t}-t$ are still to be taken and so the 4-qubit
local collision operator $U_{\text{col}}$ only needs to be applied to the
remaining qubits relevant for encoding the ‘directly connected’ velocity
states as given in Equation (17).

With this logic we can define a collision operator per time step $t$ as

$U_{\text{col},t}^{\text{tot}}=\underbrace{U_{\text{col}}\otimes\dots\otimes
U_{\text{col}}}_{\text{$c$ collision
operations}}\otimes\underbrace{I\otimes\dots\otimes I}_{\text{identity
operations}},$ (20)

where $c=2(N_{t}-t)^{2}+2(N_{t}-t)+1$ and the identity operations are added to
avoid dimensionality issues. In practice no operation will be applied on the
qubits encoding velocity states not ‘in the vicinity of $x$’ within $N_{t}-t$
time steps.

Our space-time encoding enables different manners of implementing the
streaming step. It can easily be seen that the way the streaming method should
be implemented differs per time step $t$ depending on which positions will be
‘in the vicinity of $x$’ in the next time step as well. At the first time step
it is important for (almost) all qubits to be streamed to a very specific
position, whereas in the last time step it is only important for the qubits
$q_{0}$, $q_{1}$, $q_{2}$ and $q_{3}$ to end up in the correct state. For the
example shown below we are only considering a total of one step to be taken
(i.e. $N_{t}=1$) and so we only need to consider the speeds that will stream
to location $x$ in one time step. In this case that means that streaming
consists of performing a swap operation between the following qubit pairs
$q_{0}$ and $q_{12}$, $q_{1}$ and $q_{17}$, $q_{2}$ and $q_{6}$ as well as
$q_{3}$ and $q_{11}$.

Also in general (i.e. $N_{t}>1$), the streaming step can be implemented by a
combination of swap gates. Following the same in-the-vicinity-of-$x$ argument
as was used for the collision step, a total of

$n_{\text{swap}}(t)=4+\sum_{i=1}^{N_{t}-t}16i=8\left(N_{t}-t\right)^{2}+8\left(N_{t}-t\right)+4$
(21)

swap gates are required to update as many velocity-encoding qubits in time
step $t$, whereby these swap operations can be performed largely in
parallel.777In each time step the swap operations in the 4 (or generally
speaking $m$) different directions can be performed in parallel. Furthermore
the swap operations for the velocities in the same direction but not in the
same ‘line of streaming’ can all be performed in parallel. Therefore we only
need to take into account the velocities in the same line of streaming and the
depth of the circuit is determined by the longest ‘line of streaming’, which
is equal to $T-t$. In each layer of the swap operations at least half of the
$T-t$ velocities can be swapped to the correct position. Therefore a total of
$\log_{2}\left(T-t\right)$ swap operations needs to be performed in the $t$-th
time step. The depth of the streaming circuit at time $t$ will amount to

$d_{\text{str}}(t)=\log_{2}\left(T-t\right)$ (22)

swap operations at time $t$.

Figure 7: Illustration of the computational basis state encoding for D2Q4.
Figure 8: Illustration of the space-time encoding for D2Q4 for a single time
step.

##  5 Conclusion

In this paper we have shown that current data encoding methods considered for
quantum Boltzmann methods do not allow for treating both streaming and
collision as unitary quantum operations. We have provided both a mathematical
proof of its impossibility, and insight into the physical properties of the
system and encodings that lead to this behavior. Using this insight we
subsequently developed a new space-time data encoding method that does allow
for both streaming and collision to be implemented as a unitary operation.
This paper should serve as a guideline on where (not) to look for successful
quantum encodings of the lattice Boltzmann and other QCFD methods.

## References

  * [Yep98] Jeffrey Yepez  “Quantum Computation of Fluid Dynamics”  In _Quantum Computing and Quantum Communications: lecture notes in computer science_ , 1998  URL: https://www.phys.hawaii.edu/~yepez/papers/publications/pdf/1999LectNotesCompSciVol1509Pg35.pdf
  * [Yep01] Jeffrey Yepez  “Quantum Lattice-Gas Model for computational fluid dynamics”  In _Physical Review E_ , 2001  URL: https://journals.aps.org/pre/abstract/10.1103/PhysRevE.63.046702
  * [YB01] Jeffrey Yepez and Bruce Boghosian  “An efficient and accurate quantum lattice-gas model for the many-body Schrödinger wave equation”  In _Computer Physics Communications_ , 2001  URL: https://www.phys.hawaii.edu/~yepez/static/papers/pdf/2002CompPhysCommVol146No3%2C15Pg280.pdf
  * [Yep02] Jeffrey Yepez  “Quantum Lattice-Gas Model for the Burgers Equation”  In _Journal of statistical physics_ , 2002  URL: https://link.springer.com/article/10.1023/A:1014514805610
  * [Pra+03] Marco A. Pravia, Zhiying chen, Jeffrey Yepez and David G. Cory  “Experimental demonstration of Quantum Lattice Gas Computation”  In _Quantum Information Processing_ , 2003  URL: https://link.springer.com/article/10.1023/A:1025835216975
  * [CW12] Andrew M. Childs and Nathan Wiebe  “Hamiltonian Simulation Using Linear Combinations of Unitary Operations”  In _Quantum Information & Computation_, 2012  URL: https://arxiv.org/pdf/1202.5822.pdf
  * [Krü+17] Timm Krüger et al.  “The lattice Boltzmann method”  Springer, 2017  URL: https://link.springer.com/book/10.1007/978-3-319-44649-3
  * [Bud20] Ljubomir Budinski  “Quantum algorithm for the advection-diffusion equation simulated with the lattice Boltzmann method”  In _Quantum Information Processing 2021_ , 2020  URL: https://link.springer.com/article/10.1007/s11128-021-02996-3
  * [TS20] B.. Todorova and R. Steijl  “Quantum algorithm for the collisionless Boltzmann equation”  In _Journal of Computational Physics, 409, 109347_ , 2020  DOI: http://dx.doi.org/10.1016/j.jcp.2020.109347
  * [Bud21] Ljubomir Budinski  “Quantum algorithm for the Navier-Stokes equations by using the streamfunction-vorticity formulation and the lattice Boltzmann method”  In _International Journal of Quantum information_ , 2021  URL: https://arxiv.org/abs/2103.03804
  * [Fre+22] Gabriel Freytag, João V.. Lima, Paolo Rech and Philippe O.. Navaux  “Impact of Reduced and Mixed-Precision on the Efficiency of a Multi-GPU Platform on CFD Applications”  In _Lecture Notes in Computer Science_ , 2022  URL: https://link.springer.com/chapter/10.1007/978-3-031-10542-5_39
  * [MVS22] Y. Moawad, W. Vanderbauwhede and R. Steijl  “Investigating hardware acceleration for simulation of CFD quantum circuits”  In _Front. Mech. Eng_ , 2022  DOI: 10.3389/fmech.2022.925637
  * [SM22] Merel A. Schalkers and Matthias Möller  “Efficient and fail-safe collisionless quantum Boltzmann method”  In _Arxiv preprint_ , 2022  DOI: https://doi.org/10.48550/arXiv.2211.14269
  * [Ste23] René Steijl  “Quantum Circuit Implementation of Multi-Dimensional Non-Linear Lattice Models”  In _MDPI: Applied Sciences_ , 2023  DOI: https://doi.org/10.3390/app13010529


Towards DIS at N4LO

A. Basdew-Sharma1$\star$

1 Nikhef, Amsterdam, The Netherlands

*<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

We will report on an ongoing effort towards calculating the N4LO perturbative
QCD corrections to the DIS total inclusive cross-section. We are developing a
method based on differential equations and series expansion in the inverse
Bjorken parameter. As a byproduct our calculation should also deliver analytic
or at least precise numerical approximations for the four-loop splitting
functions.

###### Contents

  1. 1 Introduction
  2. 2 Setup
  3. 3 Method
    1. 3.1 Mellin moments
    2. 3.2 Recursive algorithm
      1. 3.2.1 One-loop
      2. 3.2.2 Multiple loops
  4. 4 Results
  5. 5 Conclusion

##  1 Introduction

With the forthcoming high luminosity upgrade at the Large Hadron Collider
(LHC), the experimental uncertainty of the production of Higgs particles is
expected to reduce from $\Delta_{\sigma_{H}}\sim 8\%$ (run 2) to
$\Delta_{\sigma_{H}}\sim 3\%$ [1]. To achieve the desired theoretical
uncertainty for the predictions, one needs to perform calculations at N3LO
precision in QCD perturbation theory. These calculations, in turn, formally
require the four-loop splitting functions, which can be extracted from deep
inelastic scattering (DIS) calculations at N4LO. However, due to the
nonexistence of the four-loop splitting functions, predictions for a number of
key observables have been calculated with the three-loop splitting functions
[2, 3], resulting in a higher theoretical uncertainty.

Secondly, one can imagine that it will also be desirable to obtain more
precise parton distribution functions (PDFs). Most of the data used to
describe PDFs comes from DIS [4]. Therefore, to get a better idea of the PDFs,
we need a better extraction of DIS data, for which N4LO DIS calculations can,
at least, be useful. Moreover, the scale-evolution of the PDFs is governed by
the splitting functions [4, 5], which is another incentive to extract
splitting functions at N4LO, if one wants to reduce the theoretical
uncertainty of the PDFs.

We are developing a new method to calculate the N4LO perturbative QCD
corrections to the DIS total inclusive cross-section. As a byproduct we will
be able to extract the four-loop splitting functions. The method is based on
differential equations and series expansion in the inverse Bjorken parameter
($\omega$). In particular, we will calculate the Mellin moments in a recursive
manner, to compute coefficient functions from which we can build the structure
functions and extract splitting functions.

In the past the full $\omega$ dependance of the coefficient functions was
obtained from a finite number of Mellin moments. To mention some results: in
1999 Moch and Vermaseren calculated the DIS structure functions at two loops
[6]. In 2004 Moch, Vermaseren and Vogt calculated the three-loop splitting
functions in the non-singlet [7] and singlet case [8]. After that it took Moch
et al thirteen years to calculate the four-loop non-singlet splitting
functions using an operator product expansion method [9]. In 2018 numerical
approximations for the singlets were obtained using moments with $n<10$ [10].
In the same year Herzog et al calculated the $n=2,3$ moments of the non-
singlets at five-loops using the R* method [11].

In our method we exploit the analytic behaviour of the moments at $\omega\to
0$ to calculate expansions of the master integrals in $\omega$, by using the
differential equations obtained from the Integration-By-Parts (IBP) reduction
tables (with Kira [12] and FIRE [13], including LiteRed [14, 15, 16]). We use
FORM [17, 18, 19] to perform the symbolic manipulations and FORCER [20] to
obtain our boundary conditions of the differential equations. Our ultimate
goal is to calculate the N4LO DIS QCD corrections and to obtain the four-loop
splitting functions. However, in order to do so, we first need to validate our
algorithm. To that end our short term goal is to first check known results at
three-loop which have been checked before [21, 22, 23, 24], as well as
unchecked results like the polarised splitting functions.

This report is structured in the following way: in section 2 we will sketch
the basic setup. In section 3.1 we will reconstruct the coefficient function
from a finite number of Mellin moments. After that we will describe our
recursive algorithm in section 3.2 and we wil list some intermediate results
in section 4. We will conclude in section 5.

##  2 Setup

This section is largely based on [25] and [26]. Our starting point will be the
hadronic tensor $W_{\mu\nu}$ as sketched in figure 1.

Figure 1: Depicted we see the one loop DIS hadronic tensor, with $p$ and $q$
the momenta of the hadron and photon. The black dots are the hadrons and the
red line is a cut.

By using the optical theorem, one can relate the hadronic tensor to the
forward scattering amplitude $T_{\mu\nu}$ as

$W_{\mu\nu}(p,q)=2\text{ Im }T_{\mu\nu}(p,q),$ (1)

effectively removing the cut. In the following we will focus on the partonic
level of the forward scattering amplitude to calculate the coefficient
functions, as illustrated in figure 2. It is known that the structure
functions can be constructed from the coefficient functions [27].

Figure 2: To build the coefficient functions, we focus on the partonic part of
the forward scattering diagram, highlighted in green.

##  3 Method

Having sketched the setup, we can now proceed to our method. First we will
show how to relate the Mellin moments to the forward scattering amplitude, as
explained in [25] and [26]. After that we will sketch our recursive algorithm
to obtain the Mellin moments.

###  3.1 Mellin moments

The forward scattering amplitude can be written as

$T_{\mu\nu}(p,q)=-\left(\eta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}}\right)T_{1}(\omega,Q^{2})+\left(p_{\mu}-\frac{q_{\mu}p\cdot
q}{q^{2}}\right)\left(p_{\nu}-\frac{q_{\nu}p\cdot
q}{q^{2}}\right)T_{2}(\omega,Q^{2}),$ (2)

where $\omega\equiv x^{-1}$ is the inverse Bjorken parameter and
$Q^{2}=-q^{2}$. The functions $T_{1}(\omega,Q^{2})$ and $T_{2}(\omega,Q^{2})$
have branch cuts for $\omega\leq-1$ and $\omega\geq 1$ (figure 3).

Figure 3: The functions $T_{1}(\omega,Q^{2})$ and $T_{2}(\omega,Q^{2})$ have
branch cuts for $\omega\leq-1$ and $\omega\geq 1$.

Using Cauchy’s theorem, one can relate the derivatives at $\omega=0$ to a
contour integral around it ($C_{0}$ in figure 4) as

$T^{(n)}_{i}(Q^{2})=\frac{1}{n!}\frac{d^{n}T_{i}(\omega,Q^{2})}{d\omega^{n}}\bigg{|}_{w=0}=\oint_{C_{0}}\frac{d\omega}{2\pi
i}\frac{T_{i}(\omega,Q^{2})}{\omega^{n+1}}.$ (3)

However, the contour around the origin can be deformed into two contours
($C_{1}$ in figure 4) around the branch cuts.

Figure 4: The contour $C_{0}$ around the origin can be deformed into two
contours $C_{1}$ around the branch cuts.

Now we can use the fact that the functions $T_{1}$ and $T_{2}$ are symmetric
in $\omega$, such that we obtain

$T_{i}^{(n)}(Q^{2})=\frac{(1+(-1)^{n})}{2\pi
i}\int_{1}^{\infty}d\omega\frac{Disc_{\omega}T_{i}(\omega,Q^{2})}{\omega^{n+1}}.$
(4)

Notice that the discontinuity across the branch cuts is the imaginary part of
the forward scattering amplitude. Therefore, by again using the optical
theorem for even $n$, we can go back to $x$-space to obtain

$T_{i}^{(n)}(Q^{2})=\frac{1}{\pi}\int_{0}^{1}dx\
x^{n-1}W_{i}(x,Q^{2})=\frac{1}{\pi}\mathcal{M}_{n}[W_{i}(Q^{2})].$ (5)

Thus we see that from a finite number of Mellin moments, one can obtain the
full $\omega$ dependance of the coefficient functions. This has been used
extensively in the past to calculate splitting functions [6, 7, 8]. In the
next section we will exhaust the analytic behaviour at $\omega\to 0$ to
calculate the solutions to the loop integrals.

###  3.2 Recursive algorithm

In this section we will explain our algorithm to recursively obtain solutions
to the loop integrals involved in the calculations. We will first tackle the
one-loop problem and then generalise to higher orders.

####  3.2.1 One-loop

One can scalerize the loop integrals to sets of scalar integrals which belong
different topologies. At one-loop order (see figure 5) there is only one
topology,

$\text{Topo}(A1,A2,A3)=\int
d^{D}K\left(\frac{1}{(K^{2})^{A1}}\right)\left(\frac{1}{((P-K)^{2})^{A2}}\right)\left(\frac{1}{((Q+K)^{2})^{A3}}\right).$
(6)

Figure 5: The topology of topo1 is illustrated. Here we used the fact that one
can exploit the symmetry $Q\leftrightarrow-Q$ to reduce the diagrams to a
single topology.

Subsequently we used the program Kira [12] to reduce all integrals in the
topology to a set of master integrals, $M_{1}=\text{Topo}(1,0,1)$ and
$M_{2}=\text{Topo}(0,1,1)$. Using the reduction tables, we can obtain a
differential equation for $M_{2}$

$\frac{\partial M_{2}}{\partial\omega}=-\frac{\epsilon}{1+\omega}M_{2}.$ (7)

Notice that $M_{1}$ does not depend on $\omega$. Our boundary condition for
$M_{2}$ is $M_{2}(\omega\to 0)=M_{1}$. Next we use our series solution ansatz
for the master integral

$M_{2}(\omega,\epsilon)=\sum_{n=0}^{\infty}\omega^{n}a_{n}(\epsilon)$ (8)

and notice that the left hand side of eq. (7) always is of a higher order in
$n$ for fixed $\omega$ than the right hand side. This allows us to derive a
recursive relation for the coefficients

$a_{n}=-\frac{\left(\epsilon+n-1\right)a_{n-1}}{n}.$ (9)

Thus, knowing the boundary conditions, one can build the full
$\omega$-dependent solutions to the master integrals.

####  3.2.2 Multiple loops

There arise three problems when one continues to higher order loops: 1. the
boundary conditions should be calculated (which is hard in general), 2. one
needs to deal with higher order poles in the differential matrix, and 3. the
reductions will require a lot of computing power and runtime.

The boundary conditions are obtained by setting $\omega\to 0$. This is
equivalent to taking $P\to 0$. However, notice that the latter limit reduces
all the diagrams to massless self energy diagrams. Luckily for us these have
all been calculated (see [15, 16, 28, 29, 30, 31, 32]) and stored by our group
members using FORCER [20]. Thus, all boundary conditions are known.

For the second problem, let us consider a toy topology with three masters,

$\frac{d}{d\omega}M_{i}=A_{ij}\ M_{j},$ (10)

with

$\displaystyle A$ $\displaystyle=\begin{pmatrix}1&0&0\\\
\frac{a}{\omega^{4}}&1&0\\\
\frac{b}{\omega^{3}}&0&\frac{c}{\omega}\end{pmatrix}.$ (11)

Notice the higher order poles in the differential matrix $A$. We will
transform these away order by order by using a scaling transformation $T$

$\displaystyle\begin{split}\vec{M}^{\prime}&=T\vec{M}\\\
A^{\prime}&=TAT^{-1}-T\partial_{\omega}T^{-1}.\end{split}$ (12)

$T_{ij}$ is a diagonal matrix in which the $\text{i}^{th}$ entry equals the
(negative) exponent in $\omega$ of the highest order pole in the
$\text{i}^{th}$ row in $A_{ij}$. Thus, to scale away the pole in the term
$\frac{a}{\omega^{4}}$ in equation eq. (11), $T$ is chosen as

$\displaystyle T_{4}$ $\displaystyle=\begin{pmatrix}1&0&0\\\ 0&\omega^{4}&0\\\
0&0&1\end{pmatrix},$ (13)

such that

$\displaystyle A^{\prime}$ $\displaystyle=\left(\begin{array}[]{ccc}1&0&0\\\
a&\frac{\omega+4}{\omega}&0\\\ \frac{b}{\omega^{3}}&0&\frac{c}{\omega}\\\
\end{array}\right).$ (17)

Subsequently we have to scale away the $O(\omega^{-}{3})$ term. We do this, by
transforming with

$\displaystyle T_{3}$ $\displaystyle=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&\omega^{3}\end{pmatrix},$ (18)

resulting in

$\displaystyle A^{\prime\prime}$
$\displaystyle=\left(\begin{array}[]{ccc}1&0&0\\\
a&\frac{\omega+4}{\omega}&0\\\ b&0&\frac{c+3}{\omega}\\\ \end{array}\right).$
(22)

Therefore we are left with $O(\omega^{-1})$ poles. As a remark we would like
to mention that it is not possible to scale away the $O(\omega^{-1})$ poles in
the same manner. The final transformation matrix $T=T_{3}\cdot T_{4}$ becomes

$\displaystyle T$ $\displaystyle=\left(\begin{array}[]{ccc}1&0&0\\\
0&\omega^{4}&0\\\ 0&0&\omega^{3}\\\ \end{array}\right).$ (26)

As a consequence, according to eq. (12), $M^{\prime\prime}_{2}$ and
$M^{\prime\prime}_{3}$ start at order $O(\omega^{4})$ and $O(\omega^{3})$
respectively in their series solution. We conjecture that a diagonal
transformation matrix is sufficient to map away deeper poles from the
differential matrix. So far this has always worked, but a rigorous proof is
lacking.

Lastly we need to deal with the $O(\omega^{-1})$ terms left in the
differential matrix. To that end it is convenient to pull out a factor of
$\omega^{-1}$ from the differential matrix, such that is it not singular
anymore and obtain

$\frac{d}{d\omega}M_{i}=\frac{A_{ij}}{\omega}M_{j}.$ (27)

By expanding the masters and differential matrix as

$\displaystyle\begin{split}M_{i}=\sum_{k=0}^{\infty}c_{i}^{(k)}\omega^{k},\\\
A_{ij}=\sum_{k=0}^{\infty}A_{ij}^{(k)}\omega^{k},\end{split}$ (28)

we can, by separating out the coefficient, build a new recursive relation for
the coefficients $c_{l}$ of the master integrals as

$c_{l}^{(k)}=\left(B_{li}^{(k)}\right)^{-1}\sum_{0\leq
k{{}^{\prime}}\leq(k-1)}A_{ij}^{(k-k{{}^{\prime}})}c_{j}^{(k{{}^{\prime}})},$
(29)

in which

$B_{ij}^{(k)}\equiv\left(k\ \delta_{ij}-A_{ij}^{(0)}\right).$ (30)

We would like to stress that this is a recurrence relation: for a fixed $l$
the order $k$ on the left hand side of eq. (29) is always larger than
$k^{\prime}$ on the right hand side. We have implemented this relation in a
highly parallelized code in FORM [17, 18, 19].

##  4 Results

In this section we will list some intermediate results to display the
performance of the algorithm.

#loops | #topos
---|---
1 | 1
2 | 3
3 | 32
4 | 364

The number of topologies at each loop ramps up quite fast from three loops and
onward.

#loops | time
---|---
1 | 3.7s (8 cores)
2 | 90s (24 cores)
3 |  week (64 cores)
4 | ?

The reduction time of the topologies at each loop is small at one and two
loop, however it increases to a week at three loops. We have not run the four-
loop reduction yet.

#loops | time
---|---
1 | 2
2 | 31
3 | 2989
4 | ?

For the number of masters at each loop we see a similar pattern. The number of
masters increases by two orders of magnitude at three loops.

Our toughest calculation for the coefficients so far is for topology 31 at
three loops with 156 masters. It took ten hours to calculate 100 coefficients
on a computer using 58 cores. It is uncertain whether we will be able to run
the full four-loop reduction on our machine (64 cores, 1TB RAM memory). At
least feasible should be $n_{f}^{2}$ pieces. We are in the process of
discussing possible improvements of the reduction with one of the Kira authors
J. Usovitsch.

##  5 Conclusion

We have developed an algorithm to calculate the DIS cross-section at four
loops. To introduce our new algorithm we first explained our setup to
calculate the coefficient functions. We related the Mellin moments to the
forward scattering amplitude. We showed how to obtain differential equations
for the resulting master integrals, using IBP reduction tables. At one loop
this was sufficient to derive a recursive relation. At higher orders in
perturbation theory, we needed to cope with extra poles in $\omega$ in the
differential matrix. To this end we scaled away the poles, by transforming the
master integrals and we modified the recurrence relation.

To conclude, our algorithm has been checked to work at one- and two-loop
orders. We are finishing up the calculation at three loops. We will start the
four-loop reductions for the $n_{f}^{2}$ pieces soon.

## Acknowledgements

I would like to thank Andreas Vogt, Jos Vermaseren, Andreas Pelloni and Franz
Herzog for their collaboration on this project.

##### Funding information

This work is supported by the NWO Vidi grant 680-47-551.

## References

  * [1] S. Amoroso _et al._ ,  _Les Houches 2019: Physics at TeV Colliders: Standard Model Working Group Report_ ,  In _11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches_ (2020), 2003.01700.
  * [2] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger,  _Higgs Boson Gluon-Fusion Production in QCD at Three Loops_ ,  Phys. Rev. Lett. 114, 212001 (2015),  10.1103/PhysRevLett.114.212001,  1503.06056.
  * [3] F. A. Dreyer and A. Karlberg,  _Vector-Boson Fusion Higgs Production at Three Loops in QCD_ ,  Phys. Rev. Lett. 117(7), 072001 (2016),  10.1103/PhysRevLett.117.072001,  1606.00840.
  * [4] A. Accardi _et al._ ,  _A Critical Appraisal and Evaluation of Modern PDFs_ ,  Eur. Phys. J. C 76(8), 471 (2016),  10.1140/epjc/s10052-016-4285-4,  1603.08906.
  * [5] S. Alekhin, J. Blümlein, S. Moch and R. Placakyte,  _Parton distribution functions, $\alpha_{s}$, and heavy-quark masses for LHC Run II_,  Phys. Rev. D 96(1), 014011 (2017),  10.1103/PhysRevD.96.014011,  1701.05838.
  * [6] S. Moch and J. A. M. Vermaseren,  _Deep inelastic structure functions at two loops_ ,  Nucl. Phys. B 573, 853 (2000),  10.1016/S0550-3213(00)00045-6,  hep-ph/9912355.
  * [7] S. Moch, J. A. M. Vermaseren and A. Vogt,  _The Three loop splitting functions in QCD: The Nonsinglet case_ ,  Nucl. Phys. B 688, 101 (2004),  10.1016/j.nuclphysb.2004.03.030,  hep-ph/0403192.
  * [8] A. Vogt, S. Moch and J. A. M. Vermaseren,  _The Three-loop splitting functions in QCD: The Singlet case_ ,  Nucl. Phys. B 691, 129 (2004),  10.1016/j.nuclphysb.2004.04.024,  hep-ph/0404111.
  * [9] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,  _Four-Loop Non-Singlet Splitting Functions in the Planar Limit and Beyond_ ,  JHEP 10, 041 (2017),  10.1007/JHEP10(2017)041,  1707.08315.
  * [10] A. Vogt, F. Herzog, S. Moch, B. Ruijl, T. Ueda and J. A. M. Vermaseren,  _Anomalous dimensions and splitting functions beyond the next-to-next-to-leading order_ ,  PoS LL2018, 050 (2018),  10.22323/1.303.0050,  1808.08981.
  * [11] F. Herzog, S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,  _Five-loop contributions to low-N non-singlet anomalous dimensions in QCD_ ,  Phys. Lett. B 790, 436 (2019),  10.1016/j.physletb.2019.01.060,  1812.11818.
  * [12] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch,  _Integral reduction with Kira 2.0 and finite field methods_ ,  Comput. Phys. Commun. 266, 108024 (2021),  10.1016/j.cpc.2021.108024,  2008.06494.
  * [13] A. V. Smirnov and F. S. Chuharev,  _FIRE6: Feynman Integral REduction with Modular Arithmetic_ ,  Comput. Phys. Commun. 247, 106877 (2020),  10.1016/j.cpc.2019.106877,  1901.07808.
  * [14] R. N. Lee,  _Group structure of the integration-by-part identities and its application to the reduction of multiloop integrals_ ,  JHEP 07, 031 (2008),  10.1088/1126-6708/2008/07/031,  0804.3008.
  * [15] R. N. Lee,  _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ (2012),  1212.2685.
  * [16] R. N. Lee,  _LiteRed 1.4: a powerful tool for reduction of multiloop integrals_ ,  J. Phys. Conf. Ser. 523, 012059 (2014),  10.1088/1742-6596/523/1/012059,  1310.1145.
  * [17] J. A. M. Vermaseren,  _New features of FORM_ (2000),  math-ph/0010025.
  * [18] J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga,  _FORM version 4.0_ ,  Comput. Phys. Commun. 184, 1453 (2013),  10.1016/j.cpc.2012.12.028,  1203.6543.
  * [19] M. Tentyukov and J. A. M. Vermaseren,  _The Multithreaded version of FORM_ ,  Comput. Phys. Commun. 181, 1419 (2010),  10.1016/j.cpc.2010.04.009,  hep-ph/0702279.
  * [20] B. Ruijl, T. Ueda and J. A. M. Vermaseren,  _Forcer, a FORM program for the parametric reduction of four-loop massless propagator diagrams_ ,  Comput. Phys. Commun. 253, 107198 (2020),  10.1016/j.cpc.2020.107198,  1704.06650.
  * [21] J. Ablinger, J. Blumlein, S. Klein, C. Schneider and F. Wissbrock,  _The $O(\alpha_{s}^{3})$ Massive Operator Matrix Elements of $O(n_{f})$ for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_,  Nucl. Phys. B 844, 26 (2011),  10.1016/j.nuclphysb.2010.10.021,  1008.3347.
  * [22] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider and F. Wißbrock,  _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_,  Nucl. Phys. B 886, 733 (2014),  10.1016/j.nuclphysb.2014.07.010,  1406.4654.
  * [23] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider,  _The 3-loop pure singlet heavy flavor contributions to the structure function $F_{2}(x,Q^{2})$ and the anomalous dimension_,  Nucl. Phys. B 890, 48 (2014),  10.1016/j.nuclphysb.2014.10.008,  1409.1135.
  * [24] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider,  _The three-loop splitting functions $P_{qg}^{(2)}$ and $P_{gg}^{(2,N_{F})}$_,  Nucl. Phys. B 922, 1 (2017),  10.1016/j.nuclphysb.2017.06.004,  1705.01508.
  * [25] G. F. Sterman,  _An Introduction to quantum field theory_ ,  Cambridge University Press,  ISBN 978-0-521-31132-8 (1993).
  * [26] D. Bonocore, E. Laenen and R. Rietkerk,  _Unitarity methods for Mellin moments of Drell-Yan cross sections_ ,  JHEP 05, 079 (2016),  10.1007/JHEP05(2016)079,  1603.05252.
  * [27] J. A. M. Vermaseren, A. Vogt and S. Moch,  _The Third-order QCD corrections to deep-inelastic scattering by photon exchange_ ,  Nucl. Phys. B 724, 3 (2005),  10.1016/j.nuclphysb.2005.06.020,  hep-ph/0504242.
  * [28] T. Ueda, B. Ruijl and J. A. M. Vermaseren,  _Calculating four-loop massless propagators with Forcer_ ,  J. Phys. Conf. Ser. 762(1), 012060 (2016),  10.1088/1742-6596/762/1/012060,  1604.08767.
  * [29] T. Ueda, B. Ruijl and J. A. M. Vermaseren,  _Forcer: a FORM program for 4-loop massless propagators_ ,  PoS LL2016, 070 (2016),  10.22323/1.260.0070,  1607.07318.
  * [30] B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt,  _Four-loop QCD propagators and vertices with one vanishing external momentum_ ,  JHEP 06, 040 (2017),  10.1007/JHEP06(2017)040,  1703.08532.
  * [31] P. A. Baikov and K. G. Chetyrkin,  _Four Loop Massless Propagators: An Algebraic Evaluation of All Master Integrals_ ,  Nucl. Phys. B 837, 186 (2010),  10.1016/j.nuclphysb.2010.05.004,  1004.1153.
  * [32] R. N. Lee, A. V. Smirnov and V. A. Smirnov,  _Master Integrals for Four-Loop Massless Propagators up to Transcendentality Weight Twelve_ ,  Nucl. Phys. B 856, 95 (2012),  10.1016/j.nuclphysb.2011.11.005,  1108.0732.