<figure><img src="image_1.png" /></figure>

HAL
open science

# A Fully Quantum Algorithm for Hydrodynamic Lattice Gas Cellular Automata

Niccolo Fonio, Giuseppe Di Molfetta, Pierre Sagaut

<figure><img src="image_2.png" /></figure>

## To cite this version:

Niccolo Fonio, Giuseppe Di Molfetta, Pierre Sagaut. A Fully Quantum Algorithm for Hydrodynamic Lattice Gas Cellular Automata. 2024. hal-04543971

HAL Id: hal-04543971

https://hal.science/hal-04543971

Preprint submitted on 12 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Quantum collision circuit, quantum invariants and quantum phase estimation procedure for fluid dynamic lattice gas automata

Niccolo Fonio, Pierre Sagaut, Giuseppe Di Molfetta

► To cite this version:

Niccolo Fonio, Pierre Sagaut, Giuseppe Di Molfetta. Quantum collision circuit, quantum invariants and quantum phase estimation procedure for fluid dynamic lattice gas automata. 2025. hal-04543971v2

HAL Id: hal-04543971

https://hal.science/hal-04543971v2

Preprint submitted on 27 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A Fully Quantum Algorithm for Hydrodynamic
Lattice Gas Cellular Automata

Niccolo Fonio¹,²*, Giuseppe Di Molfetta¹ and Pierre Sagaut³

¹Aix Marseille Univ, Université de Toulon, CNRS, LIS, Marseille, France.

²Università di Torino, Dipartimento di Fisica, Torino, Italia.

³Aix Marseille Univ, Centrale Marseille, CNRS, M2P2 UMR 7340, 38 rue Joliot-Curie, Marseille, 13013, France.

*Corresponding author(s). E-mail(s): niccolo.fonio@edu.unito.it;

Contributing authors: giuseppe.dimolfetta@lis-lab.fr;

pierre.sagaut@univ-amu.fr;

**Abstract**

Lattice Gas Cellular Automata (LGCA) are a computational model widely known and applied for the simulation of many physical phenomena. Their implementation requires an amount of resources and operations which scale linearly versus the system size and number of time steps. We propose a quantum-pointers-based quantum algorithm able to simulate LGCA while exhibiting an exponential advantage in space complexity and a number of quantum operations independent from the system size. We propose a collision circuit for the FHP lattice-gas automata considering the 2-, 3-, and 4-body collisions. These are implemented with two methodologies that suggest the procedure for finding quantum circuits for LGCA with more collisions. We also propose a phase estimation algorithm to retrieve information about a single cell, whose application can be expanded for implementing other collisions. A general methodology to identify the invariants associated to quantum LGCA is also proposed.

**Keywords:** Cellular Automata, Quantum Computing, Quantum Circuits, Quantum Advantage, Fluid-dynamic
# Quantum collision circuit, quantum invariants and quantum phase estimation procedure for fluid dynamic lattice gas automata

Niccolò Fonio<sup>a,b</sup>, Pierre Sagaut<sup>a</sup>, Giuseppe Di Molfetta<sup>b</sup>

<sup>a</sup>M2P2, Aix-Marseille University, Central Marseille, M2P2 UMR 7340, 38 rue Joliot-Curie, Marseille, 13013, , France

<sup>b</sup>LIS, Aix-Marseille university, LIS UMR 7020, Campus de Luminy, 163 avenue de Luminy, Marseille, 13288, , France

## Abstract

Lattice Gas Cellular Automata (LGCA) is a classical numerical method widely known and applied to simulate several physical phenomena. In this paper, we study the translation of LGCA on quantum computers (QC) using computational basis encoding (CBE), developing methods for different purposes. In particular, we clarify and discuss some fundamental limitations and advantages in using CBE and quantum walk as streaming procedure. Using quantum walks affect the possible encoding of classical states in quantum orthogonal states, feature linked to the unitarity of collision and to the possibility of getting a quantum advantage. Then, we give efficient procedures for optimizing collisional quantum circuits, based on the classical features of the model. This is applied specifically to fluid dynamic LGCA. Alongside, a new collision circuit for a 1-dimensional model is proposed. We address the important point of invariants in LGCA providing a method for finding how many invariants appear in their QC formulation. Quantum invariants outnumber the classical expectations, proving the necessity of further research. Lastly, we prove the validity of a method for retrieving any quantity of interest based on quantum phase estimation (QPE).

**Keywords:** Cellular Automata, Lattice gas, Quantum Computing, Quantum Circuits, Fluid-dynamic

Email address: niccolo.fonio@lis-lab.fr (Niccolò Fonio)

Preprint submitted to Computers and fluids

January 16, 2025
# 1 Introduction

Since Computational Fluid Dynamics (CFD) algorithms for solving Navier-Stokes (NS) equations have been developed through the last decades depending on available technologies, the application of quantum computing is of major interest for fluid dynamic simulation [1], seeking a quantum advantage. For Quantum CFD there are, currently, several promising proposals that follow two main approaches.

The first one is based on computing the numerical solutions of NS equations, and it is addressed as the general field of solving linear or non-linear PDEs on a quantum computer. Several algorithms were proposed concerning this approach. This approach generally starts with the discretization of the PDE of interest and uses a quantum solver to find a solution. Concerning this possibility, a remarkable research [2] obtained an exponential advantage using the quantum amplitude estimation algorithm (QAEA) [3], that has been recently implemented in an efficient way [4], and has been also applied specifically to Burger's equation [5]. Quantum Walks have also been used as quantum numerical tools to simulate hydrodynamical equations [6, 7]. Other proposals regarding quantum numerical methods use a hybrid classical-quantum for accelerating parts of classical solvers or variational approach [8–11].

Alongside numerical methods for solving PDE, we find lattice methods, the aim of which is to simulate a fluid. Lattice Gas Cellular Automata (LGCA) and Lattice Boltzmann Methods (LBM) [12] belong to this class. LGCA considers a gas of particles propagating on a lattice with specific rules at a microscopic level, and its application for CFD started with 2D incompressible flows, with the landmark contributions of Frisch, Hasslacher, and Pomeau (FHP) [13] and Wolfram [14]. However, FHP for 2 dimensions and FCHC (Face Centered Hyper-Cube) [15–17] for 3 dimensions recovered NS equations with additional nonphysical spurious terms. The main issues were anomalous invariants, the lack of Galilean invariance, and the dependency of the pressure on the density. For FCHC these problems were overcome with a multi-speed model [18]. Regardless of more sophisticated models, the issues of LGCA moved the attention to LBM methods, which consider the propagation of probability density functions at the mesoscopic level according to the discrete velocity Boltzmann equation. Both LGCA and LBM were applied in several fields of application in physics, chemistry, biology [19–30]. The main technical issue of LGCA and LBM is scalability: large grids are often needed to simulate physical phenomena, and they quickly reach the limits of classical computers. In addition, LGCA and LBM modeling quantum phenomena may take advantage of quantum computers for many-body simulations. Scalability and many-body simulations are two of the reasons that motivated the formulation of Quantum LGCA [31–34]. However, the formulation proposed for Quantum LGCA does not provide a gate-based quantum algorithm, considering a stop-and-go procedure, with quantum collision and classical streaming, providing no advantage. It was experimentally proved with NMR [35], but has not been further developed. Up to now, to the best knowledge of the authors, no general quantum algorithm for LGCA has been proposed. For LBM, on the other side, there are promising proposals adopting Carleman linearization [36–38], and hybrid classical-quantum methods [39]. In particular, in [36], they consider unitary and non-unitary collision and streaming, comparing the space and computational complexity, without getting any advantage over classical
# 1. Introduction

A Lattice Gas Cellular Automata (LGCA), addressed as $DnQv$ model, is a gas of particles propagating in a discretized space of $n$ dimensions where the components of the system (particles) can exhibit $v$ discrete velocities [1]. Each $DnQv$ model consists of a lattice and a discrete evolution rule. The lattice gas is encoded with a bit string, i.e. *cell*, for each lattice point. In every cell each bit represents the presence of a particle with the corresponding velocity, as in Fig.1. The evolution rule is made of a collision step when

<figure><img src="image_3.png" /><figcaption>Figure 1: Example of a 1D lattice gas with $N$ cells and 3 velocities $[-1, 0, +1]$. The presence of a rest particle is a black dot. The presence of a moving particle is an arrow.</figcaption></figure>

particles scatter in each cell enforcing some conservation laws, and a streaming step when particles move to neighboring cells according to their velocity. The LGCA models we are going to consider, introduced in Sec.2, are a 1D model, namely D1Q3, and a 2D model, namely D2Q6. The latter has been studied by Friesch, Hasslacher and Pomeau [2], and thus it is addressed as FHP. The evolution of a LGCA is capable of simulating various physical nonlinear phenomena [39] and, most interestingly for us, they can be used for computational fluid dynamics (CFD). In particular FHP has been the first LGCA model capable of retrieving Navier-Stokes-like equations. For the opportunities of $DnQv$ models for CFD and given the possible advantages that quantum computing (QC) is showing directly in this field [4, 5], this paper focuses on the quantum computing formulation of LGCA, named herinafter QLGCA.

It is possible to look at QLGCA from different perspectives. First, they can be seen as a subclass of quantum cellular automata (QCA) [26], considering qubits instead of bits. This has been the contribution of some seminal works [10, 11], that provide no advantage in terms of computational efficiency, but show the possibility of simulating PDEs using quantum systems. Another perspective is to look at QLGCA and compare them with the corresponding classical numerical methods. In this case, we are interested in getting an advantage in terms of computational resources (number of qubits
algorithms. The problems caused in the LBM approach by the non-linearity of the col-
lision term and then the need for a linear embedding technique, and the necessity to
stream PDF do not appear in LGCA. In LGCA, since we stream bits instead of PDF,
the evolution in the case of detailed balance is linear: the non-linearity is given by the
complex behavior of the whole system. This is the essential difference in the use of the
encoding of the space regarding previous works. LGCA achieved good results for the
simulation of several phenomena, despite issues with the initial formulation, overcome
by more sophisticated models. Furthermore, LGCA with square lattice (HPP) and
the FHP model (on the hexagonal lattice) we consider were found to be P-complete
[40], thus predicting HPP or FHP-III (that considers specific collisions) is equivalent
to calculating an arbitrary Boolean circuit. These reasons motivated the present work,
whose aim is to present a quantum gate-based algorithm for the simulation on a quan-
tum computer of LGCA, giving a specific example for the FHP model, since it can be
considered as a starting point.

In this article we propose a quantum algorithm for simulating any detailed-balance
LGCA getting a quantum exponential advantage in the number of qubits adopting the
encoding of space proposed in [36] rephrased from a computational point of view with
quantum pointers, and freeing the computational complexity from any dependency
on the size of the system. This opens up the possibility of simulating LGCA with
a number of cells never approachable with a classical computer. The perspectives
for using this advantage will rely also on the capability of future research to model
chemical, biological, and physical phenomena with LGCA. Further developments for
LBM can be considered.

## 2 Results

In the following, we present our main results: 1) a quantum gate-based algorithm for LGCA simulation with unitary collision and streaming that provides an exponential advantage in space complexity, 2) a specific collision circuit for FHP in its first formulation, 3) a phase estimation algorithm to recover the single cell quantities such as mass and momentum and 4) considerations about quantum invariants of a quantum LGCA and a method to calculate their amount.

Quantum LGCA

We first consider that any LGCA with detailed balance can be implemented on a quan-
tum computer. This is ensured by the Hermiticity of its transition matrix, from which
we can formulate the unitary collision operator Ĉ. This is possible assuming the stan-
dard Computational Basis Encoding (CBE) of the cell. CBE is explained in Methods
section. The encoding of the lattice is performed with quantum pointers. For a DnQv
model (with *n* the spatial dimension and *v* the number of discrete velocities, respec-
tively), such pointers coincide with *n* superposed quantum registers for the position
and one superposed quantum register for the velocity set. Then, using generalized Tof-
foli gates, each term of the superposition gets entangled, according to the presence of
a particle. The entanglement allows us to perform the streaming among the pointers
of the position and it also moves the particle that is pointed. As the first application,
and operations), and we aim to find the best representation of the classical
system for a quantum algorithm. In this sense, much attention has been
given to *quantum Lattice Boltzmann Methods* (QLBM).

LBM [28, 29] is a family of classical lattice-based numerical methods that were born from LGCA, and solved some issues for hydrodynamic simulations (e.g. preservation of Galilean invariance), becoming one of the most used CFD methods. Instead of microscopic particles, in LBM we stream mesoscopic probability distribution functions, solving the lattice Boltzmann equation under specific assumptions and simplifications. A significant difference with LGCA is that, in most cases, the collision term of LBM is nonlinear. The first QLBM was developed by Yepez [6, 8, 9, 7] and showed that it was possible to replicate LBM on a quantum computer and to solve PDEs. Since then, different approaches have been proposed and are in rapid expansion. One approach is to linearize the non-linear collision operator of LBM adopting Carleman linearization [12, 13]. Another one gets an exponential advantage in space complexity, requiring measurement and reinitialization at each time step [14, 15]. These first alternatives are probabilistic since they rely on a linear combination of unitaries [16, 17]. Alternative encodings and applications can be found, each with different advantages and drawbacks [19, 18]. The principle of each QLBM proposed is to interpret the quantum amplitudes as the classical probability distribution functions. Despite being one of the most promising foreseen ways for QCFD, the non-unitarity of the process and the necessity of measurements and reinitialization hinder a real advantage at the current state of the art. Thus, it is worth looking at features and methods for QLGCA, as we are going to discuss in this paper.

In particular, LGCA can exhibit 2 advantages over LBM for QC imple-
mentation: (P1) the collision consists of a non-deterministic or deterministic
correspondence of input/output states, and (P2) the collision is the same
in each cell. The property (P2) allows to leverage quantum parallelism, de-
tailed in Sec.3.1, for a computational advantage. The property (P1) can
be advantageous in QC assuming the computational basis encoding (CBE)
that we introduce in Sec.3.1. This encoding of classical states into quan-
tum states allows to carry out deterministic collisions as unitary operations,
and non-deterministic collisions as unitary operations followed by measure-
ments. Unitary operations are the fundamental operations used in QC, and
performing unitary collision and streaming with a small number of qubits is
a crucial aspect of the quantum advantage. Usually, this is not possible for
the non-unitarity of the collisions, as in [14, 15, 5]. QLGCA with CBE, on
let us consider the D2Qv model, where each node lives on a two-dimensional grid. In such situation, we consider three pointers : one for the row $|r\rangle$, one for the column $|c\rangle$, and a last one for the velocity $|v\rangle$. The presence of the particle is coded by one qubit $|n\rangle$ and we also need a supplementary cell register $|q\rangle$

<figure><img src="image_4.png" /><figcaption>Fig. 1 Quantum gate-based circuit for simulation of LGCA on a quantum computer. $I_0$ is the initialization of the lattice made of generalized Toffoli gates. $C_0$ brings on a parallel ordered register of the occupation qubits: a computational basis encoding is adopted. $C$ is the collision operator that depends on the LGCA to simulate. $S$ is the streaming operator, characterized by a sequence of controlled adders that depend on the boundary conditions, number of dimensions, and number of quantum pointers</figcaption></figure>

Pointers' registers undergo a series of Hadamard gates for creating a superposition of $Nv$ indices, where $N$ is the number of cells in the lattice and $v$ is the number of velocities of the LGCA. We entangle the pointers to each particle in the lattice with $\hat{I}_0$, which is performed by applying generalized Toffoli gates. Thus, we get the following state for the lattice

$$|\Psi\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} |r\rangle |c\rangle \sum_{k=0}^{v-1} |k\rangle |n_k\rangle_{r,c} |q=0\rangle \quad (1)$$

Where $|\Psi\rangle$ is the state of the lattice, $N_r$ and $N_c$ are the number of rows and columns. Finally, we apply $\hat{C}_0$, which brings all the occupation qubits in an ordered cell register $|q\rangle$.

$$|\Psi\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} |r\rangle |c\rangle \sum_{k=0}^{v-1} |k\rangle |0\rangle |q_0 \dots q_v\rangle_{r,c} \quad (2)$$

This operation can be made with controlled swaps and a quantum adder for superposition states [41]. Now, after $C_0$ we apply the collision to the $q$ register. The cell register is initialized with CBE, so the unitary collision deriving from the transition matrix previously mentioned is able to perform the collision step with $\hat{C}$. Then, we recreate the state of Eq.1 with the updated values of qubits. This operation, called $C_1$ takes up to $O(2^v)$ controlled operations. We conjecture that this cost can be lowered depending on the specific LGCA to simulate. This is a central feature of the algorithm to be made more efficient for the implementation of more precise models such as multi-speed models [18]. As a final step, we stream with $\hat{S}$, composed by a series of controlled adders previously mentioned. Each velocity streams differently. The explicit
the other hand, shows a true correspondence between linearity of the model and unitarity of the collision. This is promising for the development of a *multi-timestep* algorithm, thus without need for measurement and reinitialization as computational parts of the algorithm. However, it presents other challenges.

In fact, the CBE we use has already shown some limitations, as outlined in [18]. It was proved that using CBE for a D1Q3 model does not allow for streaming and collision representing the space with a logarithmic number of qubits. We show in Sec.3 that this limitation goes beyond the CBE used in [18], including any encoding of the classical states in a set of orthonormal quantum states, also called *orthogonal states encoding* (OSE). Considering that these limitations hinder the achievement of a quantum advantage at the current state, our results prove that this is not strictly linked to the CBE, but to the coexistence of CBE and the streaming procedure using quantum walks. Thus, future research may still prove an advantage considering CBE or OSE to which all the methods in this paper can be applied.

Regarding the property (P1), we said that we can translate the classical collision process into a unitary operator. The decomposition of this operator into quantum gates, which are the basic operations in QC, is necessary for executing the algorithm on real devices and can be expensive [30]. We propose and prove the validity of different algorithmic and quantum algorithmic methods for carrying out the collision of D1Q3 in Sec.4 and some collisions of FHP in Sec.5 with an optimal decomposition in quantum gates. In particular for FHP we use an operator whose default decomposition using Qiskit[38] needs an order of $10^4$ operations, while our methods find an overall quantum circuit that uses just 291 operations, giving an improvement of ≈ 100 times.

In addition to computational costs, another crucial aspect when transitioning from the classical model to the quantum model regards the conservation laws. LGCA are based on preserving specific quantities, called *invariants*, such as mass and momentum. In QC, quantities are considered to be Hermitian operators, called *observables*. An observable is conserved, being called *quantum invariant*, if it commutes with the collision operator [20]. In classical and quantum LGCA there can be some *spurious invariants* [1, 33, 32, 21], which are anomalous conserved quantities that can affect the behavior of the system with additional coupled conservation equations. In Sec.3 we develop a method for counting the number of quantum invariants given a collision operator, and we discover that spurious quantum invariants are numerous compared to the classical case. This finding questions
form of each operation is provided in the discussion. This performs the mechanism of any LGCA with $O(\log_2(N)) + O(v)$ qubits and $O(I_0(N, v, T)) + O(2^vT)$ number of quantum gates. The cost of $\hat{I}_0$ depends on the lattice to be initialized, while the second term represents the cost of the evolution, which is independent of the system size. In fact, the collision and streaming operators act in parallel on each cell and on each pointer : that is the advantage of quantum parallelism. This leads to an exponential advantage in space and performs the evolution with a number of operations independent of the system size.

## Quantum collision circuit

LGCA works with a discrete evolution separated into two steps: the collision step and the streaming step. The collision step rearranges the configuration of the cells: particles scatter within the site and change direction according to the conservation laws. The collision step can be seen as a rule that sends the state $s$ to the state $s'$, where $s'$ can be decided according to conservation laws to be applied. With the computational basis encoding, the collision step becomes an operator wit the following effect

$$ \hat{C}|s\rangle = |s'\rangle $$

Love [42] tried to figure out the operations to apply in order to get the desired collision, without implementing an overall circuit for the desired collision. As we show in Fig.2, we propose a series of operations that is able to perform 0-momentum collisions. The reasoning behind this collision algorithm is expanded in 3.3. It is based on a logic formulation of the 3-body collision and an equivalence class formulation for 2- and 4-body collisions. This collision subroutine takes advantage of quantum parallelism to lower the number of operations, which becomes independent of the system size. However, it is not possible to apply non-uniform random collisions, which is a fundamental aspect to avoid the classical model being chiral. The solution to the randomness of collision is left as a future perspective. The application of this collision operation is endorsed for the computational basis encoding. For applying it, is more straightforward to use pre-collision and post-collision operations involving the velocity register $|v\rangle$, the occupation register $|n\rangle$, and the cell register $|q\rangle$. While the pre-collision is composed of a series of controlled swaps and superposition quantum adder, the post-collision can be practiced in two ways. The first one is easier, and it considers reinitialization of the cell register with a complexity of $O(vt)$ where $t$ is the number of time steps. Reinitialization can be avoided with additional qubits, which becomes $O(vt)$. This also allows us to store the state of the lattice at each time step, but increases the space complexity, conserving the advantage over the classical algorithm.

## Quantum Phase Estimation

Regarding CBE, we propose a quantum phase estimation algorithm for retrieving general information about the cell without measuring it. We define, as in [42], quantum observables for the mass and momentum in x- and y-direction $\hat{m}, \hat{P}_x, \hat{P}_y$. These are generally sums of single-qubit gates, which are Hermitian. Thus, we can measure the
the "simple" translation from classical quantities to Hermitian operators, i.e. quantum observables, and opens perspectives never tackled in previous works, highlighting the need of further research on this topic.

The existence of quantum observables and the unexpected number of quantum invariant show the wide possibilities of QLGCA, going to the core of interpreting information in LGCA. Thus, the last aspect we address regards the retrieving of information. Classically, information is accessible at any moment, making it feasible to know the state of the lattice at each time step. In QC retrieving information is possible through measurements, that directly affect the quantum state. Procedures of (re)initialization and measurements were proposed for QLBM algorithms [14, 15, 27]. These procedures, despite being promising, rely on efficient protocols that have yet to be discovered and are necessary for assessing quantum advantage. In our paper, we introduce a novel approach for retrieving quantities of interest, such as mass and momentum, using a *quantum phase estimation* (QPE) algorithm [35, 30]. This algorithm allows information to be accessed probabilistically during the computation without necessarily measure the quantum state of the cell. This does not solve the problem of measurement and reinitialization for an advantageous encoding, but is a general method applicable to any quantity of interest for any algorithm using CBE.

To summarize, the contributions of this paper are fourfold:

* expand the fundamental limitations of adopting CBE and quantum walk streaming protocol;

* we develop an optimal collision circuit for FHP and a novel collision circuit for D1Q3;

* we develop a new method for counting conserved quantities in quantum LGCA, applying it to D1Q3 and FHP and showing unexpected results;

* propose novel QPE protocols for retrieving important physical quantities.

The paper is structured as follows: in Sec.2 we introduce extensively the classical models; in Sec.3 we introduce the CBE and the encodings of the space, discussing the advantages and the methods we develop; in Sec.4 we apply the methods to D1Q3 and in Sec.5 to FHP.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Collisional circuit for 0-momentum collisions of FHP. These are performed depending on a conditional qubit $|b\rangle$ that gets flipped if the input state is a collisional state, and an additional qubit $|a\rangle$ that introduces the non-deterministic character of 2- and 4-body collisions. This circuit uses the respective invariances of collisional states.</figcaption></figure>

mass or other observables of interest $\hat{O}$ defining $\hat{U} = e^{2i\pi\hat{O}}$. This will turn the detection of those quantities of interest into a phase estimation problem.

## Quantum invariants

The last result that opens some perspectives comes with the generalization of invariants. Classically a conserved quantity is a local quantity $O = f(\vec{x}, t)$ such that $f(\vec{x}, t) = f(\vec{x}, t+1)$, e.g. mass ($m = \sum_{i=0}^{v-1} m_i n_i(\vec{x}, t)$) and momentum ($\vec{p} = \sum_{i=0}^{v-1} \vec{c}_i n_i(\vec{x}, t)$). We can find a quantum counterpart of these quantities as proposed by Love [42]. This generalizes to the definition of a *quantum invariant*. A quantum invariant is an observable that commutes with the collision operator. Based on this definition, we find trivially, thanks to CBE, mass, and momentum as the sum of single quantum gates. We can calculate the number of quantum invariants by counting the number of independent solutions to a linear system, as explained in 3.5. This definition brings the number of quantum invariants to be higher than classical invariants. This is proved for D1Q3, where the quantum invariants were found to be 14. Anomalous classical invariants are a problem in the retrieval of NS behavior. The quantum invariants are not expected to be problematic since the evolution mimics the classical evolution. However, the physical reason behind these anomalous quantum invariants and their role in this quantum simulation of a classical system is left as a perspective for future research.

## Perspectives

This article leaves several perspectives open. First, a way to characterize and optimize the initialization procedure is yet to be found. Since $I_0$ consists of a series of generalized Toffoli gates, it should be possible to rely on lowering the number of controls to initialize more cells with fewer control qubits. Symmetric or regular initialization could take advantage of quantum parallelism as well. Otherwise, an initial mechanism with collisions that do not preserve mass and momentum, followed by streaming, can be considered. In the second place, an improvement of $\hat{C}_0$ and $\hat{C}_1$ is necessary in order to get rid of the dependency on the number of velocities. In third place, future research should investigate the possibility of implementing a unitary collision operator for an LGCA with a semi-detailed balance. Another perspective is related to developing quantum circuits for the implementation of different collisions and different DnQv models. From this general implementation of the algorithm, it is possible to develop a
## 2. Classical models

### D1Q3

The first LGCA we consider is D1Q3. D1Q3 is a 1-dimensional lattice gas with 3 velocities: a right-moving particle of mass 1, a left-moving particle of mass 1, and a rest particle of mass 2. Applying an exclusion principle, such that only one particle per velocity per site is allowed, each site is represented with a bit-string $[n_2n_1n_0]$. Conventionally, $n_2, n_1, n_0$ represent the presence of a particle with velocity, respectively, -1, 0, 1. This system can simulate diffusion, and it is the easiest model to perform mass- and momentum-conserving collision, ensuring non-linearities. The collision allowed (and its inverse) consists of splitting a rest particle into two opposite-moving particles (the opposite is merging two opposite-moving particles into a rest particle). This can be represented as $[101] \leftrightarrow [010]$. Fig.2 provides an example of a 1-time step evolution of this gas.

<figure><img src="image_4.png" /><figcaption>Figure 2: D1Q3 evolution example. In this case, each cell is represented with 3 bits $[n_2n_1n_0]$. Before the collision in the cell at $x = 1, 2, 5$ there are collisional configurations, so they evolve according to the chosen collision. All the other cells are not affected by the collision. The streaming takes place according to respective velocities</figcaption></figure>

Classically, we can write pseudocode as in Alg.1 for the collision step. Here we can anticipate that the quantum analogous of a code or pseudocode is a quantum gate-based circuit, interpreting operations as series of quantum gates. Quantities such as mass $m(x) = \sum_i n_i(x)$ and momentum $p(x) = n_0(x) - n_1(x)$ are conserved. Averaging over neighbors we get mass density $\rho(x)$ and momentum density $\vec{u}(x)$. The behavior of these quantities in the continuum limit follows differential equations, which depend on the conservation laws applied.
circuit for any DnQv model. The crucial question, which is a general problem of using a superposition, is how to retrieve the classical information. We conjecture that it should be possible to get global information about the system without perturbing the quantum state of the lattice. This will be possible with a deep modeling approach. The natural development of this algorithm may bring to the implementation of quantum LBM models, improving the existing proposals. Quantum circuits for different boundary conditions are yet to be found. The usage of the phase estimation algorithm can be expanded for the implementation of other collisions or retrieving information from the cells: it is possible to explore the global information about the system. The generalization of quantum invariants leaves open questions about their physical meaning and their possible role in these quantum simulations. As we said before, this algorithm makes it possible to simulate large grids, with an Avogadro's number of cells for example, using a relatively small amount of qubits, depending on the number of velocities to consider. With the current perspective on the technological development of quantum computers, simulations with this algorithm can take place in the future. The availability of large grids that quantum computing admits could flourish beginning with physical, chemical, and biological modeling with LGCA.

# 3 Methods

## 3.1 Introduction to FHP

A lattice gas cellular automata LGCA can generally be addressed as a DnQv model, where $n$ is the space dimension and $v$ is the number of velocities (bits) per cell, and it is characterized by the following elements [14] [12]

* A regular arrangement (lattice) of $N$ cells in $n$ dimensions

* Each cell has a discrete finite set of possible states represented by $v$ bits

* Each cell state gets updated synchronously with collision and streaming operations

In this article we will address the FHP model, also called D2Q6 model, in two dimensions with 6 velocities on a hexagonal lattice, considering it as a first step toward more sophisticated models such as Digital physics [18]. The FHP set of velocities is the following.

$$ \vec{c}_i = \left( \cos\left(\frac{\pi}{3}i\right), \sin\left(\frac{\pi}{3}i\right) \right) \text{ for } i = 0, 1, \dots, 5 $$

Thus, each node of the lattice consists of a cell whose occupation state is defined by 6 bits $n_i(\vec{r}, t) \in \{0, 1\}$ for $i = 0, 1, \dots, 5$. $n_i(\vec{r}, t) = 1$ if a particle with velocity $\vec{c}_i$ is present at node in position $\vec{r}$ at time $t$, and $n_i(\vec{r}, t) = 0$ otherwise. The complete formulation of the model can be found in [12] [14]

Essentially, the LGCA can be computed as a grid of bits whose values evolve according to different rules. We will encode this grid of bits into a quantum system, performing the same evolution assured with Computational Basis Encoding. The dynamics consist of two steps: *collision* and *streaming*. An example of an evolution step is represented in Fig.3. The collision step modifies the state of each cell. For simulating NS equations, the fundamental constraint is to enforce elastic collisions, i.e.
<table><thead><tr><th colspan="2">Algorithm 1 Collision D1Q3 algorithm</th></tr></thead><tbody><tr><td><b>for</b> cell ∈ lattice <b>do</b></td><td></td></tr><tr><td>m ← get mass(cell)</td><td>▷ Retrieving of quantities</td></tr><tr><td>p ← get momentum(cell)</td><td></td></tr><tr><td><b>if</b> cell = [010] or cell = [101] <b>then</b></td><td></td></tr><tr><td>cell ← new state</td><td>▷ Collision process</td></tr><tr><td><b>end if</b></td><td></td></tr><tr><td><b>end for</b></td><td></td></tr></tbody></table>

## FHP

The second LGCA we consider, which is more interesting for computational fluid dynamics (CFD) simulations, was proven by Frisch, Hasslacher, and Pomeau (FHP)[2] and Wolfram [3] to converge in the continuous limit to the Navier-Stokes equations. It is a 2-dimensional LGCA on a triangular lattice, showing a collision that changes the configurations of the cell, followed by a streaming step. An example of the evolution step is drawn in Fig. 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Evolution step for Frisch, Hasslacher, and Pomeau model (FHP). Before collision is the starting state of the lattice. In the center, we see that the particles in cells with the collisional states of Table.1 are rearranged. In the last part we see that particles have been streamed to neighboring cells applying periodic boundary conditions</figcaption></figure>

The set of velocities {$\vec{c}_i$} is defined as follows, and each $c_i$ is linked to $n_i$ and ordered as [$n_5n_4n_3n_2n_1n_0$]

$$ \vec{c}_i = \left( \cos\left(\frac{\pi}{3}i\right), \sin\left(\frac{\pi}{3}i\right) \right) \quad \text{for } i = 0, 1, \dots, 5 \qquad (1) $$
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Evolution step for FHP. Before collision is the starting state of the lattice. In the centre we see that the particles in cells with the collisional states of **Table 1** are rearranged. In the last part we see that particles have been streamed to neighbouring cells applying periodic boundary conditions</figcaption></figure>

to introduce a collision operator such that the mass, i.e. $m = \sum_{i=0}^{6} n_i$, and the total momentum of the cell, i.e. $\vec{p} = \sum_{i=0}^{6} \vec{c}_i n_i$ and the kinetic energy $K = \sum_{i=0}^{6} \vec{c}_i^2 n_i$ are conserved. This is needed for ensuring the conservation laws, that are treated with a Chapman-Enskog expansion, and retrieve at the second order NS-like equations. It is worth noting that on FHP's hexagonal lattice kinetic energy is conserved as soon as mass (i.e. the number of particles) is conserved since all discrete velocities have the same modulus. The conserved quantities in a classical simulation constrain the system to evolve preserving more quantities than needed. When we use a Quantum LGCA, we can define quantum invariants [42]. We show in 3.5 that quantum invariants of D1Q3 are 14 instead of 3, applying a generalized method. It is an open question to know if these quantum invariants have physical meaning. We conjecture they do not affect the simulation since the simulation of the classical system just depends on the collision, not explicitly on the conservation laws. Examples of possible collisions ensuring conservation of mass, momentum, and kinetic energy are given in **Table 1**

<figure><img src="image_7.png" /><figcaption>**Table 1** Collisions of FHP. These are the 0-momentum collisions of the model</figcaption></figure>

In general, as we said, any collision preserving mass (and then automatically kinetic energy) and momentum can be considered, because the form of the final equations does not depend on the collision. The transport coefficients, such as viscosity, generally depend on the collisions that we implement. The collisions given in **Table 1** are those that we are going to implement in a quantum circuit in 3.3.
<figure><img src="image_6.png" /><figcaption>Table 1: These are the 0-momentum collisions of the model FHP model. We can notice that $B_2$ and $B_4$ are probabilistic rotations of 120° or 240°, taking place with equal probability, while $B_3$ is a rotation of 180°. $B_2$ and $B_4$ collisional states are invariant under $B_3$, and vice versa</figcaption></figure>

Among different elastic collisions, we focus on 0-momentum collisions represented in Table1. These are the collisions that involve cells with a configuration that has momentum $p = 0$, and they are represented in Table 1. A rest particle can be added for stability of the numerical simulations and for adding new collisions, causing a faster thermalization to the equilibrium of the solution.

If any cell in the lattice is found to be in a collisional state represented in Table 1, then the scattering process takes place and particles are rearranged. We can see these collisions as rotations. Rotations can also have invariant states. Trivially the full and empty cells are invariant respect to any rotation. We notice, as reported in Table 3, that the states involved in 2- and 4-body collisions, called $B_2$ and $B_4$, are invariant respect to a rotation of 180°, that is going to be the 3-body collision called $B_3$. These invariances are going to be used for the quantum implementation of the algorithm. Specific features used for finding the optimized quantum circuit are reported in Sec.5.1.

Likewise the case of D1Q3, we show a quantum gate-based implementation for the collision step involving 0-momentum collisions. The quantities that are conserved classically are, as before, mass $m(x) = \sum_{i=0}^{5} n_i(x)$ and momentum $\vec{p}(x) = \sum_{i=0}^{5} \vec{c}_i n_i(x)$. These quantities are averaged over neighbouring cells and their evolution is mapped to the continuous limit with a Chapman-Enskog expansion. This retrieves Euler's equations at the first order and Navier-Stokes equations at the second order [1, 2]. We define the same quantities as quantum observables, verifying that they are conserved and we show how we can count the total amount of quantum invariants given a specific collision.
We provide here the quantum version of an LGCA. First of all, we consider the computational basis encoding (CBE) for the bit-string of the cell, thus, any bit in the bit-string of a cell is turned into a qubit:

$$n_i \rightarrow |n_i\rangle \forall i \in \{0, \dots, v\}$$

In the second place, we consider the transition matrix of an LGCA, which is defined as a square real matrix A with elements $a_{s,s'}$ such that $\sum_{s'} a_{s,s'} = 1$ and $a_{s,s'} \ge 0 \forall s, s'$. In particular, $a_{s,s'}$ is the probability that the cell state starts from $s$ and evolves into $s'$. A deterministic transition matrix is a transition matrix where $a_{s,s'} = 1$ or $0, \forall s, s'$: there is not any randomness and each input state evolves into itself or into only one other state. A symmetric transition matrix is a transition matrix such that $a_{s,s'} = a_{s',s} \forall s, s'$, representing a detailed-balanced LGCA: the probability of $s \rightarrow s'$ is the same probability that $s \rightarrow s'$. FHP is a detailed-balanced LGCA. Since the transition matrix of a detailed-balanced LGCA is real and symmetric, it is trivially Hermitian. Since it is Hermitian, we can say that $e^{iA}$ is unitary. Thus, we can simulate any detailed-balanced LGCA on a quantum computer using $\hat{C} = e^{iA}$ as the collision operator that acts on a quantum register given by CBE. If the LGCA is deterministic, the transition matrix corresponds to the matrix representation of the collision operator in the orthonormal basis provided by CBE. Thus, we can state the following property **Proposition 1.** The collision operator of any LGCA with detailed balance can be simulated on a quantum computer

We can generally decompose any unitary operator in a quantum circuit with different methodologies, such as [43]. We will implement the collision for FHP with two methods based on the model itself. The methods to be preferred depend on the physical implementation of the quantum computer being used.

In general, an LGCA without detailed balance has a transition matrix that is not symmetric and thus is not Hermitian. In this case, it can be also possible to find a unitary able to simulate on a quantum computer the corresponding LGCA, but we leave it as an open perspective.

## 3.2 Quantum D2Q6

In this section, we give the specific formulation of the different operations previously mentioned that compose the overall circuit Fig.1. With quantum pointers encoding space and CBE, we have the following state for the lattice

$$|\Psi\rangle = \frac{1}{\sqrt{N2^{v'}}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'}-1} |r\rangle |c\rangle |k\rangle |n_k\rangle_i |q=0\rangle \quad (3)$$

Where $|r\rangle$ is the register of the row index, $|c\rangle$ is the register of the column index, $|k\rangle$ is the register of the velocity, $N_r$ is the number of rows, $N_c$ is the number of columns, $v' = \log_2(v)$ approximated by excess, $|n_k\rangle_{r,c}$ is the occupation qubit of the particle with velocity $\vec{c}_k$ in position $(r, c)$, $|q\rangle$ is the cell-register, that must be in 0 state. Since the initialization procedure depends on the lattice, we give a general example for the following lattice
### 3. Encodings and methods

A classical $DnQv$ model needs $Nv$ bits to represent the state of the lattice with $N$ grid points. These are $Nv$ values of 0s or 1s, as in Fig.1. The fundamental element in quantum computing, analogous to the classical bit, is the qubit. A qubit is a vector in a Hilbert space $|\psi\rangle \in \mathcal{H}^2$ that can be written as $|\psi\rangle = \alpha|0\rangle + \beta|1\rangle$, where $|0\rangle$ and $|1\rangle$ compose the computational basis, and $\alpha, \beta \in \mathbb{C}$ and $|\alpha|^2 + |\beta|^2 = 1$. A qubit represents a quantum states where the physical system is at the same time in $|0\rangle$ and $|1\rangle$. If we carry out a measurement, the system collapses in $|0\rangle$ with probability $|\alpha|^2$, or in $|1\rangle$ with probability $|\beta|^2$. Considering column vector notation $|0\rangle = (1, 0)^T$ and $|1\rangle = (0, 1)^T$, being $|\psi\rangle = (\alpha, \beta)^T$. Multiqubit states are represented with the cross product of the single qubits, usually implicit, and can use binary notation (e.g. a three-qubit state can be $|0\rangle \otimes |1\rangle \otimes |0\rangle = |010\rangle = |2\rangle$). To change the state of a qubit we need to apply unitary operations $\hat{U}$, called gates. These can be single-qubit gates, thus representable as $2 \times 2$ unitary matrices, or v-qubit gates, representable as $2^v \times 2^v$ unitary matrices. In this section, we present possible encoding of LGCA in quantum states and methods for three purposes: optimal implementation of collision, calculation of quantum invariants, and retrieving quantities of interest.

#### 3.1. Encodings

In this paper we use the Computational basis encoding (CBE). We encode the cell populated by particles with $v$ different velocities in a $DnQv$ model with $v$ qubits. We consider the quantum states $|0\rangle$ and $|1\rangle$ for the absence and presence of a particle with respective velocity.

$$ \begin{gathered} \text{Classical encoding} \rightarrow \text{Quantum register} \\ [n_v n_{v-1} \dots n_0] \rightarrow |n_v n_{v-1} \dots n_0\rangle \end{gathered} \qquad (2) $$

An encoding with $v$ qubits per cell is used in different works for QLBM [6, 8, 9] and QLGA [27], and our methods can be applied to different encodings of the space. In particular, we can consider a linear encoding and a sublinear encoding of the space.

The linear encoding of the space utilizes $Nv$ qubits. With this encoding, the state $|\Psi(t)\rangle$ of the entire lattice at time $t$ is as follows

$$ |\Psi(t)\rangle = \bigotimes_{x=0}^{N} |\psi(x, t)\rangle = \bigotimes_{x=0}^{N} |n_v(x, t) \dots n_0(x, t)\rangle \qquad (3) $$
$$ \begin{bmatrix} 010010 & 100001 \\ 000000 & 101010 \end{bmatrix} $$

With quantum pointers encoding of space and CBE the initial state of this lattice can be written as follows

$$ |\Psi\rangle = \frac{1}{\sqrt{32}} [ |0\rangle |0\rangle (|0\rangle |0\rangle + |1\rangle |1\rangle + |2\rangle |0\rangle + |3\rangle |0\rangle + |4\rangle |1\rangle + |5\rangle |0\rangle ) |q=0\rangle + \\ |0\rangle |1\rangle (|0\rangle |1\rangle + |1\rangle |0\rangle + |2\rangle |0\rangle + |3\rangle |0\rangle + |4\rangle |0\rangle + |5\rangle |1\rangle ) |q=0\rangle + \\ |1\rangle |0\rangle (|0\rangle |0\rangle + |1\rangle |0\rangle + |2\rangle |0\rangle + |3\rangle |0\rangle + |4\rangle |0\rangle + |5\rangle |0\rangle ) |q=0\rangle + \\ |1\rangle |1\rangle (|0\rangle |1\rangle + |1\rangle |0\rangle + |2\rangle |1\rangle + |3\rangle |0\rangle + |4\rangle |1\rangle + |5\rangle |0\rangle ) |q=0\rangle ] $$

Where $|i\rangle$ is intended to be the bit-representation of $i$ for the velocity pointer.

The $I_0$ operation for this lattice will be the following one represented in Fig.4. Different combinations of controls over the same superposed registers make us access

<figure><img src="image_8.png" /><figcaption>Fig. 4 Example of state preparation with quantum pointers</figcaption></figure>

specific occupation qubits, that become entangled with the specific pointer. This is an example of massive initialization. As we said, a way to optimize this procedure is one of the open perspectives. After initializing the lattice, we act with $C_0$, the pre-collision operation, which is an operation with the following action

$$ C_0 \sum_{k=0}^{2^{v'-1}} |k\rangle |n_k\rangle |000000\rangle = \sum_{k=0}^{2^{v'-1}} |k\rangle |0\rangle |n_0 n_1 n_2 n_3 n_4 n_5\rangle $$

We first use a set of controlled-swaps to get the following state.

$$ \sum_{k=0}^{v'-1} |k\rangle |n=0\rangle |n_k 2^k\rangle $$

Then, with the quantum adder of [41], we get the sum of all these superposed terms, describing the state of the cell in the computational basis encoding. Since the
where $|\psi(x, t)\rangle$ is the state of the cell in $x$ at time $t$. This encoding does not allow for any advantage of interest for quantum simulations of classical CFD schemes. Nevertheless, this system corresponds to a Watrous quantum cellular automata [26, 34]. Thus, the methods proposed here are of interest for quantum simulations of QCA.

The sublinear encoding, on the other hand, leverages superposition, entanglement, and quantum parallelism to seek an advantage. If we have $n = \log_2 N$ qubits, we are capable of representing $N$ states at the same time. This corresponds to creating a superposition in the space register. Then, with an initialization procedure, we entangle each state of the space register to the corresponding state of the cell register, using CBE. The state of the lattice is as follows

$$|\Psi\rangle = \frac{1}{\sqrt{N}} \sum_{x=0}^{N} |x\rangle |\psi(x)\rangle = \frac{1}{\sqrt{N}} \sum_{x=0}^{N} |x\rangle |n_v(x) \dots n_0(x)\rangle \quad (4)$$

We can now carry out the collision on the cell register once, and it will affect each cell. This is called quantum parallelism [30], and can lower the cost of the operations needed and of the computational resources dramatically. However, some problems arise. With the advantageous encoding in Eq.4, we cannot perform streaming and collision as we need. This is rooted in the following assumptions:

* occupation states of the classical cell correspond to orthogonal states in the quantum register: this ensures the unitarity of the collision;

* the streaming operation is performed with a quantum walk procedure as in [14, 15];

* Different velocities must be represented in distinguishable states in the cell register, in order to apply a controlled shift for moving information in the correct direction.

Under these assumptions, coming from transposing the classical algorithm in quantum terms, it is not possible to perform unitary collision and streaming for multi-time step implementation. The full proof can be found in Appendix A, and expands the limitations already outlined in [18], where an analogous proof was given for a specific encoding. We extend that proof to any possible encoding of $|\psi(x)\rangle$ in Eq.4, going beyond CBE.
quantum adder needs two copies of the same superposed state to be summed, as said before we need a parallel copy of the registers involved that are initialized identically. This is also caused by the no-cloning theorem. After the collision, the $C_1$ operation, called the post-collision operation, has to split the different particles giving them back to the respective velocity pointers. This can be done with specific operations depending on the number of velocities. For example, if we had two velocities we would need a circuit with the following effect

$$C_1 |q_0 q_1\rangle \rightarrow \frac{|q_0 0\rangle |0\rangle + |0 q_1\rangle |1\rangle}{\sqrt{2}}$$

followed by the controlled swaps previously mentioned. This can be done with a quantum circuit using two additional qubits. The idea is to manipulate each possible configuration of the cell bringing it to the desired state. This procedure increases the number of operations, which becomes equivalent to the number of possible states of the cell. From here we have the exponential dependency on velocity $O(2^v)$ for each $C_1$ operation. The lowering of this cost is left as a future perspective. The streaming operator for FHP will be the one displayed in Fig.5. Moving a particle is interpreted as changing the corresponding pointer according to the velocity. In particular, referring to the grid, we will have the following analogies

$$ \begin{aligned} \vec{c}_0 &\Rightarrow \text{rows} + 0, \text{cols} + 1 \\ \vec{c}_1 &\Rightarrow \text{rows} - 1, \text{cols} + 0 \\ \vec{c}_2 &\Rightarrow \text{rows} - 1, \text{cols} - 1 \\ \vec{c}_3 &\Rightarrow \text{rows} + 0, \text{cols} - 1 \\ \vec{c}_4 &\Rightarrow \text{rows} + 1, \text{cols} + 0 \\ \vec{c}_5 &\Rightarrow \text{rows} + 1, \text{cols} + 1 \end{aligned} \tag{5} $$

<figure><img src="image_9.png" /><figcaption>Fig. 5 Streaming step. The convention on the velocity index is the same as before. The adders are unitary operations such that the modulo ensures the periodic boundary condition</figcaption></figure>

We see that we modify the pointers in space according to the velocity. The whole evolution of the lattice is represented as follows.
As discussed in the introduction, this is not proof that a quantum advantage cannot be obtained using CBE. It rather highlights some limitations for paving the way to new algorithms. The methods proposed in the following apply to any algorithm that includes the first assumption. The 2nd and 3rd assumptions can be adapted and changed for further research.

## 3.2. Collisional quantum circuits

Considering CBE, we give here the general methods adopted for optimal decomposition of the collision operator into a series of quantum gates, being this a challenge for optimisation of QC algorithms. The key idea is to use one or more ancillary qubits to be flipped if the cell is in a collisional configuration. Then, we apply a controlled collision with target the cell register, as shown in Fig. 4.

<figure><img src="image_7.png" /><figcaption>Figure 4: $|\psi\rangle$ is the cell register, $|a\rangle$ is the ancillary register, $\hat{V}$ is the verification procedure, flipping the ancillary register only if there is a collisional state, $\hat{C}$ is the controlled collision applied to the cell register</figcaption></figure>

For D1Q3 we use two ancillary qubits for a QPE procedure with some specific operators, which identify a specific set of states. This is an example on how we can use QPE for the verification procedure.

For FHP we leverage equivalence classes and logical exclusions, considerably lowering the decomposition cost of the collision operator. We implement 2-,3- and 4-body collisions through rotations of the cell, and the verification procedure is made with logical or equivalence class criteria. This gives the first proposed optimal overall circuit for FHP, allowing for a quantum simulation of FHP using CBE. The idea for carrying out non-deterministic collisions is to create a superposition of the two different outcomes and then make the system collapse using a measurement. This is reminiscent of a random extraction. With this method, we have different options, depending on the measurement we make. The first option is to create a superposition directly in the cell register. This means that if classically $s_0$ scatters into $s_1$ or $s_2$ with 50% probability, the collision does the following

$$ \hat{C} |s_i\rangle = |s'_i\rangle = \frac{|s_j\rangle + |s_k\rangle}{\sqrt{2}} \qquad (5) $$
$$|\Psi_0\rangle = |0\dots0\rangle |0\dots0\rangle |000\rangle |0\rangle |000000\rangle \qquad (7)$$

$$|\Psi_1\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t)\rangle |c(t)\rangle |k\rangle |0\rangle |000000\rangle \qquad (8)$$

$$|\Psi_2\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t)\rangle |c(t)\rangle |k\rangle |n_k(t=0)\rangle_{r,c} |000000\rangle \qquad (9)$$

$$|\Psi_3\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t)\rangle |c(t)\rangle |k\rangle |0\rangle |n_0 n_1 n_2 n_3 n_4 n_5 (t=0)\rangle_{r,c} \qquad (10)$$

$$|\Psi_4\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t)\rangle |c(t)\rangle |k\rangle |0\rangle |n_0 n_1 n_2 n_3 n_4 n_5 (t=1)\rangle_{r,c} \qquad (11)$$

$$|\Psi_5\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t)\rangle |c(t)\rangle |k\rangle |n_k(t+1)\rangle_{r,c} |0\rangle \qquad (12)$$

$$|\Psi_6\rangle = \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t) \oplus a_k\rangle |c(t) \oplus b_k\rangle |k\rangle |n_k(t+1)\rangle_{r,c} |0\rangle \qquad (13)$$

$$= \frac{1}{\sqrt{Nv}} \sum_{r=0}^{N_r-1} \sum_{c=0}^{N_c-1} \sum_{k=0}^{2^{v'-1}} |r(t+1)\rangle |c(t+1)\rangle |k\rangle |n_k(t+1)\rangle_{r,c} |0\rangle \qquad (14)$$

The state starts in $|\Psi_0\rangle$ with a series of $|0\rangle$ in Eq. 7. We create the superposition of the pointers Eq. 8 with a series of Hadamard gates, getting $|\Psi_1\rangle$. We initialize the system entangling the presence qubit with the pointers through the application of generalized Toffoli gates Eq. 9, getting $|\Psi_2\rangle$. An example of a circuit practicing this preparation is given in Fig. 4. We bring the state of each cell on a parallel quantum register to perform the collision Eq. 10, getting $|\Psi_3\rangle$. We perform the collision Eq. 11 bringing each cell to the updated state at $t+1$, getting $|\Psi_4\rangle$. We bring back the updated qubits to the values pointed by the row, column, and velocity pointers Eq. 12, getting $|\Psi_5\rangle$. We perform the streaming with controlled adders updating the value of the row and the column with different values depending on the velocities Eq. 13, getting $|\Psi_6\rangle$. In this step, the $a_k$ and $b_k$ are the ones previously mentioned for different velocities. For FHP this circuit is given in Fig. 5. We can reorder the sum matching the row and column pointers and get the updated positions Eq. 14. Then, we are ready to perform another time step, that will start back from Eq. 9.

### 3.3 Collisional circuit for FHP

For performing the desired collisions in our algorithm, we start defining these in terms of quantum operations, interpreting them as rotations of the cell. Then, we propose an overall circuit for the implementation of 2-3-4-body collisions with 0-momentum,
With $i, j, k \in \{0, 1, 2\}$ and different from each other. It is easy to prove that this operation is not unitary, since $\langle s'_i | s'_j \rangle \neq \langle s_i | s_j \rangle \forall i, j$. The alternative we are going to consider for calculating quantum invariants is to slightly change the classical collision, introducing the (low) probability of no collision. This means that the collision does the following

$$ \hat{C} |s_i\rangle = \frac{-|s_i\rangle + 2|s_j\rangle + 2|s_k\rangle}{3} \qquad (6) $$

This operation is unitary. This fact will be used in Sec.5.2. If you use this unitary collision, you can measure one of the qubits and the state collapses into one of the possibilities, analogously to a random extraction. This carries out the desired collision and is used for studying the quantum invariants of the model.

The alternative we use in Sec.5.1, more feasible as an algorithmic procedure, is to add an ancillary qubit $|a\rangle$ initialized to $|0\rangle$, thus performing the following operation

$$ \hat{C} |s_i\rangle |a\rangle = \frac{|s_{i+1}\rangle |0\rangle + |s_{i+2}\rangle |1\rangle}{\sqrt{2}} \qquad (7) $$

Where $i + 1 = \text{mod}_3(i + 1)$ and $i + 2 = \text{mod}_3(i + 2)$. This corresponds to the 2- and 4-body collisions of FHP. Using this method, we can measure the ancillary qubit obtaining the random extraction. This is the operation we are going to decompose optimally for the quantum collisional circuit of FHP. Measuring an ancilla does not solve the problem of reinitialization in case of sublinear encoding of the space, but offers a convenient decomposition of the unitary collision.

### 3.3. Quantum invariants

We look at the conservation of local quantities as defined by Love [20], such as mass and momentum. Classically, for the discretization of the space, some nonphysical invariants arise [21]. The purpose of many studies has been to get rid of them for getting better results. This is necessary since any additional invariant can bring to a conservation equation that couples with mass and momentum conservation, vanishing the simulation purpose of the algorithm. In our quantum computing framework we can calculate the exact number of invariants, differently from the classical case.

As said in the introduction, the quantum counterpart of classical quantities, such as mass and momentum, are *quantum observables*, which are
represented in Table1. We apply two different methodologies to find such a circuit. Each collision is implemented in two parts. The first part is the discrimination of the collisional states, and it uses two conditional qubits: $|b\rangle$ that is supposed to be $|1\rangle$ if we have a collisional input state, $|0\rangle$ otherwise, and $|a\rangle$ is used for creating the superposition for $B_2$ and $B_4$.

As we can see from Table1, collision $B_3$ can be interpreted as a rotation of 60° or 180°. Collisions $B_2$ and $B_4$ correspond to different rotations. We consider $B_2$ and $B_4$ to be implemented by one rotation of 120° and a successive rotation of 120° taking place with a probability of 0.5, resembling the model of FHP-I [12]. We point out that $B_3$ collisional states are invariant under a rotation of 120°, and $B_2$ and $B_4$ collisional states are invariant under a rotation of 180°. This will allow us to merge the circuit in a unique overall circuit for the three collisions. In terms of quantum operations, we can see rotations of the cell as a series of swaps operations, as shown in Fig.6

<figure><img src="image_10.png" /><figcaption>Fig. 6 Rotations with quantum circuits</figcaption></figure>

To implement the discrimination of collisional states of $B_3$ we consider a logic methodology, starting from the expression that can be found in [12], adapted to our convention.

$$b = (q_0 \land q_1) \& (q_1 \land q_2) \& (q_2 \land q_3) \& (q_3 \land q_4) \& (q_4 \land q_5)$$

Where $\land$ is a XOR operation and $\&$ is an AND operation. We can turn this logic expression in a quantum circuit interpreting $\land$ as CNOT and $\&$ as a generalized Toffoli gate with the cell's qubits as controls and $|b\rangle$ as a target. Then we restore the original cell and we apply the C-Swap gates for a rotation of 180° with $|b\rangle$ as a control qubit. The circuit implementing $B_3$ is showed in Fig.7

<figure><img src="image_11.png" /><figcaption>Fig. 7 $B_3$ collision</figcaption></figure>
Hermitian operators. An operator $\hat{O}$ is Hermitian if $\hat{O}^\dagger = \hat{O}$, where $\hat{O}^\dagger$ is the transpose conjugate of $\hat{O}$. Hermitian operators can be expressed as a linear combination of the Pauli operators, that span their orthonormal basis. For 1-qubit-observables these are $\mathcal{A}_1 = \{I, X, Y, Z\}$ where

$$I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \quad X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \quad Y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \quad Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \quad (8)$$

If we consider 2-qubit-observables, their orthonormal basis is given by the cross product of the basis of the 1-qubit-observables $\mathcal{A}_2 = \{II, IX, IY, IZ, XX, XY, XZ, YY, YZ, ZZ\}$ where the implicit operation is the cross product, i.e. $XX = X \otimes X$. This holds also for bigger dimensions, thus generally, a basis for $v$-qubit-observables $\mathcal{A}_v$ has $4^v$ elements, so that any quantum $v$-qubit-observable $\hat{O}$, i.e. any $2^v \times 2^v$ Hermitian matrix, is

$$\hat{O} = \sum_{i=0}^{4^v-1} \alpha_i \hat{O}_i \qquad (9)$$

where $\hat{O}_i \in \mathcal{A}_v$, and $\alpha_i$ are the corresponding coefficients.

How do observables, i.e. counterparts of classical quantities, evolve? If the state $|\psi\rangle$ undergoes an evolution $\hat{C}$, the evolved state is $|\psi'\rangle = \hat{C}|\psi\rangle$. Analogously, the evolved observables can be written as $\hat{O}' = \hat{C}^\dagger\hat{O}\hat{C}$. Considering a cell with $v$ qubits, an observable $\hat{O}$ is a quantum invariant if it commutes with the collision operator $\hat{C}$, so if it satisfies $[\hat{C}, \hat{O}] = \hat{C}\hat{O} - \hat{O}\hat{C} = 0$. We can write this property as follows, considering that for the unitarity of collision we have $\hat{C}\hat{C}^\dagger = \hat{C}^\dagger\hat{C} = I$,

$$\hat{C}^\dagger \hat{O} \hat{C} = \hat{O} \qquad (10)$$

If we look at the rhs as the evolved operator, we can clearly see the invariant property: the post-collision observable (lhs) is equal to the pre-collision observable (rhs). We can say that any operator $\hat{O}$ respecting Eq.10 is a quantum invariant. If we want to know the number of linearly independent quantum invariants for a QLGCA model, we look at the evolution of each basis element $\hat{O}_i \in \mathcal{A}_v$. The number of quantum invariants corresponds to the number of linearly independent solutions of the linear system where each equation is the conservation equation Eq.10 for $\hat{O}_i \in \mathcal{A}_v$. This translates into the following property
This circuit alters only the collisional states of $B_3$, thus we will have the following evolution

$$|101010\rangle|00\rangle \rightarrow |010101\rangle|10\rangle \\ |010101\rangle|00\rangle \rightarrow |101010\rangle|10\rangle$$

To implement $B_2$ and $B_4$, we apply a reasoning based on equivalence classes. We define an asymmetric opposite pair as a pair of opposite bit-velocities that differ from each other (e.g. $|100000\rangle$ has 1 opposite pairs for the asymmetry between $|q_0\rangle$ and $|q_3\rangle$, $|110000\rangle$ has 2 opposite pairs for $|q_0\rangle|q_3\rangle$ and $|q_1\rangle|q_4\rangle$, $|100100\rangle$ has 0 opposite pairs). We consider the equivalence class of states with 0 asymmetric pairs.

As we can see, the collisional states of $B_2$ and $B_4$ belong to this class, and the other states are invariant under these collisions. Thus, we can target $|b\rangle$ in order to be $|1\rangle$ if there are no asymmetric pairs. Then with a series of C-Swaps with $|b\rangle$ as a control, we apply the first 120° rotation. Then, to recall the non-deterministic character of the rule, we use $|a\rangle$ to create a superposition between different outcomes. We apply a controlled-H gate with $|b\rangle$ as control and $|a\rangle$, initialized to $|0\rangle$, as target. Then we apply the same controlled rotation of 120° with $|a\rangle$ as a control.

<figure><img src="image_12.png" /><figcaption>Fig. 8 $B_{2,4}$ (Head-on) collision</figcaption></figure>

The circuit in Fig. 8 alters the collisional states of $B_2$ and $B_4$, thus we will have the following evolution

$$|100100\rangle \rightarrow \frac{|001001\rangle|10\rangle + |010010\rangle|11\rangle}{\sqrt{2}} \\ |010010\rangle \rightarrow \frac{|100100\rangle|10\rangle + |001001\rangle|11\rangle}{\sqrt{2}} \\ |001001\rangle \rightarrow \frac{|010010\rangle|10\rangle + |100100\rangle|11\rangle}{\sqrt{2}} \\ |110110\rangle \rightarrow \frac{|101101\rangle|10\rangle + |011011\rangle|11\rangle}{\sqrt{2}} \\ |101101\rangle \rightarrow \frac{|011011\rangle|10\rangle + |110110\rangle|11\rangle}{\sqrt{2}}$$
**Property.** Consider a collision of a quantum DnQv model as a unitary operator $\hat{C}$. The number of linearly independent quantum invariants of a DnQv model is $l = 4^v - r$, where r is the rank of the evolution matrix M, with elements

$$M_{i,j} = \beta_{i,j} - \delta_{i,j} \qquad (11)$$

where $\delta_{i,j}$ is the Kronecker delta and $\beta_{i,j}$ is the coefficient of the $\hat{O}_j$ operator for the Pauli decomposition of the evolved $\hat{O}_i$ operator, explicitly

$$\hat{C}^{\dagger}\hat{O}_i\hat{C} = \sum_{j=0}^{4^v-1} \beta_{i,j}\hat{O}_j$$

for $\hat{O}_i \in \mathcal{A}_v$.

This methodology is applied to D1Q3 and FHP and confirms the expected conservation of mass and momentum operators, proving the existence of several more and unexpected quantum invariants.

### 3.4. Quantum Phase Estimation

The QPE algorithm is one of the standard subroutines in QC [30]. Consider a unitary operator $\hat{U}_q$ with eigenvector $|s\rangle$ and corresponding eigenvalue $e^{i2\pi q(s)}$. We start with n qubits in an additional register initialized to $|0\rangle$. These ancillae are controls of controlled-$\hat{U}_q$ operators, being $|s\rangle$ the target register, as shown in Fig.5. The state of the ancillary register acquires a rel-

<figure><img src="image_8.png" /><figcaption>Figure 5: Phase estimation algorithm for mass detection. $|s\rangle$ is the cell register, while $|q_2q_1q_0\rangle$ results in the approximated binary value of the quantity detected. $FT^\dagger$ corresponds to the inverse Fourier Transform</figcaption></figure>

ative phase according to $\hat{U}_q$ and $|s\rangle$. In the end, using the inverse quantum
$$|011011\rangle \rightarrow \frac{|110110\rangle|10\rangle + |101101\rangle|11\rangle}{\sqrt{2}}$$

As a natural consequence, these non-deterministic rules create entangled qubits within the cell. Superposition and entanglement are not physical features of a fluid dynamic simulation. Thus, a measurement of $|a\rangle$ can make the states collapse, exactly like a random extraction that takes place classically. However, there is a bias of such an implementation of non-deterministic rules: all the rotations will take place in the same way, so at each time step all the collisional states will have a rotation of the same angle, that can be different at the next time step. This represents a problem with the chirality of the models to be simulated. The consequences of this effect need to be explored in further studies. To apply this procedure we then need to reinitialize $|b\rangle$ and $|a\rangle$. This can be made by reapplying the discrimination of the collisional state part. These two circuits can be merged together for the invariances previously mentioned, as shown in Fig. 2 and their correct behavior has been verified on Quiskit.

We improved the circuits presented in [42], obtaining an overall circuit for the implementation of 0-momentum collisions. Moreover, we gave two possible methodologies to implement a collision, which can be used for other collisions. The formulation of other collisions is not trivial and requires a specific definition of equivalence classes or logic rules. The effect of avoiding the measurement of $|a\rangle$, propagating entangled qubits, is yet to be studied since the previous calculation of the evolution of the system holds if the cell can be represented as a product state. Thus, models of LGCA where there is a collision that creates entangled particles within the cell must be the subject of future studies.

## 3.4 Phase estimation algorithm

Staying within the cell, we can calculate different quantities of interest, such as momentum in the x and y direction and mass. In a QLGCA this retrieving of classical information has majorly been performed with a measurement and re-initialization of the cell [31] [39]. We show that with a phase estimation algorithm, we can avoid this stop-and-go procedure.

In [42] some specific operators whose eigenvalues can be interpreted as the respective quantities are defined. We can use these definitions since they depend on the CBE applied.

$$\hat{P}_x = Z_0 - Z_3 + \frac{Z_1 + Z_5}{2} - \frac{Z_2 + Z_4}{2} \qquad (15)$$

$$\hat{P}_y = \frac{\sqrt{3}}{2}(Z_1 + Z_2 - Z_4 - Z_5) \qquad (16)$$

$$\hat{m} = Z_0 + Z_1 + Z_2 + Z_3 + Z_4 + Z_5 \qquad (17)$$

In general, these are not unitary operators, so we cannot apply them directly to the quantum register of a cell. However, we can apply complex exponentiation to get unitary operators $U_{\hat{O}} = e^{-2\pi i \hat{O}}$. In this way detecting the corresponding eigenvalue, thus the corresponding quantity of interest, becomes a phase estimation problem. With
Fourier transform (IQFT), this relative phase translates into a binary approximation of the eigenvalue corresponding to s. Schematically, it provides the following operation

$$|0 \dots 0\rangle |s\rangle \rightarrow |\tilde{q}(s)\rangle |s\rangle \qquad (12)$$

Where $\tilde{q}(s)$ is an approximation of $q(s)$. We can look at this procedure as a way of getting the information $q(s)$ about the state $|s\rangle$ depending on the operator $\hat{U}_q$. The idea underlying the method we propose is to define a diagonal operator $\hat{U}_q$ with elements $u_q(s, s') = \delta_{s,s'}e^{iq(s)}$ for $s, s' = 0, \dots, 2^v$ for a classical quantity $q$. An example of a mass operator used for QPE is given in Eq.17. Being diagonal, the eigenstates of $\hat{U}_q$ correspond to the states of the cell adopting CBE. Thus, we can define the eigenvalues $e^{iq(s)}$ according to $s$, i.e. the state of the cell. In Eq.17 we see that states with the same mass have the same eigenvalue, so that when carrying out QPE we obtain the same result on the ancillary register. Another method for defining operators that are suitable for a QPE procedure involves the imaginary exponentiation of a quantum observable. If we consider an observable $\hat{Q}$, we can run the QPE using $\hat{U}_Q = e^{i\hat{Q}}$. Given the correspondence between classical quantities and quantum observables already outlined, this is always possible. For this paper, we preferred to use the first method as it gives more precise spectra for different values of the same quantity, as we see in the following sections.

Moving relevant physical information to an additional register allows us to measure the additional register instead of the cell register. For a linear encoding of the space, this avoids reinitialization in the case it is needed. For a sublinear encoding this QPE procedure does not solve the reinitialization obstacle and does not allow for a multi-time step implementation. However, our method allows for accessing information *during the computation*: conditional operations for which a quantity is needed can be carried out, suggesting new optimizations. For QLBM a QPE procedure was proposed [36], but used in a different way. In this other work, since the information on the probability distributions is stored in relative phases, QPE was conceived for direct information retrieving of the distributions of the cell.

In general, CBE allows us to retrieve quantities of interest with an additional register also with arithmetic operations, thus our method is not the only one for such a purpose. However, the method we propose can be applied to any quantity, going beyond the ones accessible with arithmetic operations. This is of interest for considering a change of encoding or an investigation for
this procedure, we take advantage of the fact that all the possible states of the cell are eigenstates of the mentioned operators since we are using a CBE. In the appendix, a specific example of momentum in x-direction can be found. In general, it is clear that for the CBE, the $Z_i$ operator and their combination detect the state of the cell, which can be collisional. Thus, in an alternative scheme for the discrimination of collisional states, the conditional qubit $|b\rangle$ in Fig.2 can be influenced by a phase estimation procedure for detecting a collisional configuration. For example, the operator $\hat{B}_3 = Z_0 - Z_1 + Z_2 - Z_3 + Z_4 - Z_5$ will have eigenvalue 0 if and only if there is a collisional state of $B_3$. Thus, the PE procedure together with the interpretation of $Z_i$ operators as operators that detect the presence of the i-th particle, due to the CBE adopted, can be considered another general procedure to practice a particular collision. The application of the phase estimation procedure here proposed does not provide a retrieval of classical information with quantum pointers' encoding. In fact, all the information about the cell is entangled to a superposed term of the quantum pointer register: this makes the retrieving of information particularly hard and should be investigated in future research.

## 3.5 Quantum invariants

In this section, we explain some fundamental features of this model arising from the implementation of LGCA on a quantum computer. This has a principal consequence in the conservation of local quantities [42] such as mass and momentum. Classically, for the discretization of the space, some unphysical invariants arise [44]. The purpose of many studies has been to get rid of them for getting better results. Moolvig and Texeira made extensive work [18] [45] for FCHC, which models a fluid in 3 dimensions.

Quantumly, we are dealing with observables. A physical interpretation of conserved quantities for a unitary evolution has been given [46] [47]. However, in this article, we address the computation of LGCA on a quantum computer, we do not propose a study of a quantum lattice gas. Thus, we will rely on the local conserved quantities defined by the commutation with the collision operator, as explained in the following. The conservation of these observables is the main topic of causal cellular automata, well defined in other works [48]. We define an observable $\hat{O}$ as a quantum invariant if it commutes with the collision operator $\hat{C}$, so if it satisfies $[\hat{C}, \hat{O}] = 0$. We rewrite this property expanding the commutator as follows

$$ \hat{C}\hat{O}\hat{C}^{\dagger} = \hat{O}' = \hat{O} \qquad (18) $$

We can say that any operator $\hat{O}$ respecting Eq.18 is a quantum invariant. If we want to know how many quantum invariants there are for a QLGCA we apply the following theorem.

**Theorem 2.** Consider a collision of a quantum DnQv model as a unitary operator $\hat{C}$. Given a set of Pauli gates on the product space of the v qubits involved and their linear combination, the number of quantum conserved quantities of a DnQv model is defined by the rank of the evolution matrix $M$, with elements

$$ M_{i,j} = \alpha_{i,j} - \delta_{i,j} $$
non-classical quantities, that raise interest following the quantum invariants we outlined in the previous subsection. Thus, our method proposes new tools for future research, and is validated on D1Q3 and FHP.

## 4. Application to the D1Q3 model

### 4.1. Collision circuit

The collision exchanges a rest particle and 2 moving particles, doing $[010] \leftrightarrow [101]$. As we said, CBE consists of seeing the occupational states as a set of orthogonal quantum states. Thus, the collision can be represented as a unitary matrix

$$ \hat{C} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \qquad (13) $$

A gate-based implementation of this is given in [20]. Even if Love's implementation is already optimal, we show that it is also possible to use the method in Sec.3.2 providing the verification procedure with QPE. We start from the operators $ZIZ$ and $IZZ$. If we look at their eigenvalues, practicing a QPE allows us to identify specific configurations. Dividing the space of possibilities again and again as shown in Fig. 6, we arrive in the end to have a subset of the collisional state, $\{010, 101\}$ in our case. We can call $ZIZ$ and $ZZI$ the discrimination operators, and apply a QPE scheme as shown in Fig.7. The operators were chosen to look at the different possible eigenvalues of operators on 3 qubits involving Z-gates. Operators involving X-gates and Y-gates were excluded from the CBE adopted. The collision circuit for D1Q3 in Fig.7 represents a novel method for executing collisions.

### 4.2. Quantum invariants

The first thing to do for calculating the quantum invariants is writing the collision operator. The collision operator of D1Q3 presented in [20] can be written as
where

$$ \hat{C}^{\dagger} \hat{O}_i \hat{C} = \sum_{j=0}^{4^v - 1} \alpha_{i,j} \hat{O}_j $$

This is proved by the fact that the rank of the evolution matrix correspond to the number of independent solution of the linear system where each equation is the conservation equation for the combination of Pauli gates on v qubits. To give an example, the D1Q3 collision operator presented in [42] can be written in the following way

$$ \hat{C}_{tot} = \frac{1}{4}(3III + IZZ + XXX + XYY - YXY + YYX - ZIZ + ZZI) \quad (19) $$

The previous procedure gives a rank for the evolution matrix equal to 14. This is the number of quantum invariants. Three of them are mass, momentum, and the identity. These are caused by the CBE applied, since they just "count" two quantities that are conserved by construction. The meaning of all the others has to be the subject of future studies. Thus, the formulation of a collision operator and the intuitive definition of a quantum invariant may hide several more invariants than mass and momentum for any quantum LGCA. However, the evolution of the system simulates the classical model. Thus, is probably true that the presence of multiple quantum conserved quantities does not affect the simulation model. This is true as long as this algorithm is used to simulate a classical LGCA. The quantum conserved quantities may affect the modeling of LGCA for other systems. Thus, the study of these conservation laws needs to be explored in more detail.

# 4 Conclusion

We proved that any detailed-balanced LGCA can be simulated with a quantum computer. A generalization to any LGCA, including semi-detailed balance, should be investigated in the future. Adopting the computational basis encoding for the cell state and a quantum pointers' encoding for the space allowed us to find an exponential advantage in space complexity. This can bring large-scale simulations using a relatively small number of qubits, accessible on quantum computers in the future, and with the possibility to characterize the interaction of particles within the cell in a classical or quantum way. The perspectives about simulating large-scale LGCA should be in the future to discover which phenomena can be modeled with an LGCA procedure. The advantage in computational complexity depends on the initialization: with massive initialization, needed to access each single particle, we do not get a strong advantage. Applying an alternative initialization, to be investigated in future works, brings to a cost of the initialization independent of the system size. A specific categorization of initialization procedures, also for ordered, semi-ordered, or random initial states, should be studied in the future. For the post-collisional operation, a quantum circuit without an exponential dependency on the number of velocities is to be found. For the streaming step of our algorithm using quantum pointers, we took inspiration from the classical reasoning of the algorithm, thus, we think that this can inspire also correct
<figure><img src="image_9.png" /><figcaption>Figure 6: Tree scheme to see the partition of states depending on the detected eigenvalue of associated operator for QPE.</figcaption></figure>

$$ \hat{C}_{tot} = \frac{1}{4}(3III + IZZ + XXX + XYY + \\ - YXY + YYX - ZIZ + ZZI) \qquad (14) $$

The second step is to calculate the evolution of each Pauli operator on 3 qubit space. The results are in Table B.5. From this it is possible to verify the conservation of mass *m* and momentum *p* as defined by Love :

$$ m = IIZ + 2IZI + ZII \qquad (15) $$

$$ p = IIZ - ZII \qquad (16) $$

Looking at all the evolved operators we see that there are more conserved quantities than expected. We can compute their number by calculating the evolution matrix and its rank, as explained in the Sec.3.3. For D1Q3 the rank is equal to 14, meaning that there are 50 linearly independent conserved quantities.

The conservation of all these quantities is unexpected. The quantum invariants *IZZ*, *ZIZ*, *ZZI* have explainable classical counterparts as $n_0 \oplus n_1$, $n_0 \oplus n_2$, $n_1 \oplus n_2$. This has a clear explanation in the CBE adopted. The states $|0\rangle$ and $|1\rangle$ are the eigenstates of *Z* operator. Thus, the *Z* operator is reminiscent of the presence of a particle. Others are due to a symmetry of the collision operator with respect to a change of basis, as suggested in
streaming in the presence of different boundary conditions or obstacles. We presented a collisional circuit for the implementation of FHP collisions with 0 momentum. The proposed circuit is built up with two different methodologies for the discrimination of collisional states: the first one is based on logic, and the second one is based on equivalence classes of collisional states. These two methodologies can be applied to look for implementing new and different collisions. General considerations about the usage of a phase estimation algorithm have been proposed. This can be used for different purposes, including to implement an arbitrary collision, and to retrieve information about the cell on an additional register. This procedure does not admit to retrieving information from the quantum system with a quantum pointers' encoding of the space, for which a general procedure is yet to be found. We gave a general methodology for counting the quantum conserved quantities of a QLGCA, applying it to the case of the quantum version of D1Q3, finding more invariants than expected. While the conservation of mass and momentum can be deducted from the CBE encoding, the meaning of the other quantum-conserved quantities is yet to be found.

**Acknowledgments.** This work is supported by the PEPR integrated project EPiQ ANR-22-PETQ-0007, by the ANR JCJC DisQC ANR-22-CE47-0002-01 founded from the French National Research Agency and with the support of the french government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université - A*MIDEX AMX-21-RID-011.

**Authors' contributions.** NF designed the circuit-based model supervised by GDM and PS, and was a major contributor in writing the manuscript. All authors contributed to write, read and approved the final manuscript.

**Data availability.** Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

## Appendix A Phase estimation example

We provide a specific example for momentum in x-direction. The operator defined by Love, with the present convention is Eq.15 The eigenvalues of this operator are 0, ±1, ±2, ±3. We slightly modify $P_x$ in order to get natural eigenvalues. Thus we define

$$ \hat{P}_x = Z_0 - Z_3 + \frac{Z_1 + Z_4}{2} - \frac{Z_2 + Z_4}{2} + 3I \quad (A1) $$

Thus, we get the eigenvalues to be 0, 1, 2, 3, 4, 5, 6. In this way we have integer eigenvalues, and their bit-representation is precise with 3 qubits in the register for the phase estimation algorithm. We define $\hat{U}_x = e^{-2\pi i \frac{P_x}{8}}$. We add a quantum register of 3 qubits where we store the eigenvalue, and we implement the circuit Fig.A1

We get $|p_{x_2}p_{x_1}p_{x_0}\rangle$ on the additional register, that is the bit representation of the momentum. Thus, we can measure the additional register in order to get the desired quantity. If $|u\rangle$ is a superposition of two eigenstates of $\hat{U}$, as happens after $B_2$ or $B_4$, then we get on the additional register the superposition of the corresponding eigenvalues.
<figure><img src="image_10.png" /><figcaption>Figure 7: Collision circuit for D1Q3. The first two qubits are conditional: the practicing of the collision, which is a series of Toffoli, depends on their state. This is a phase estimation algorithm with the first two as additional registers. The result, for the simplicity of the model, is deterministic.</figcaption></figure>

<table><caption>Table 2: We can see that the quantum formulation introduces a symmetry on conserved quantities. This happens also in FHP</caption><thead><tr><th>Classical counterpart</th><th>Symmetric counterpart</th></tr></thead><tbody><tr><td>III</td><td>XXX</td></tr><tr><td>IZZ, ZIZ, ZZI</td><td>XXY, YXY, YYX</td></tr><tr><td>m = IIZ + 2IZI + ZII</td><td>XXY + 2YXX + YXX</td></tr><tr><td>p = IIZ - ZII</td><td>XXY - YXX</td></tr><tr><td></td><td>IXI + XIX</td></tr></tbody></table>

Table2. However, the remaining ones are unexpected, and a clear and unique correspondence between quantum and classical quantities is left as a future perspective. Our finding is intended to highlight this feature of QC algorithm for LGA that has never been considered in such detail before.

### 4.3. Quantum phase estimation for quantities

The methodology consists of defining the quantum operator to have the desired eigenvalues for a phase estimation procedure, as explained in the methods section. For example, the matrix representation of the mass operator is the following
<figure><img src="image_13.png" /><figcaption>Fig. A1 Phase estimation algorithm for detecting momentum in x-direction</figcaption></figure>

$$PE(\hat{U}) |000\rangle |u_1\rangle = |\phi_1\rangle |u_1\rangle \\ PE(\hat{U}) |000\rangle |u_2\rangle = |\phi_2\rangle |u_2\rangle \\ PE(\hat{U}) |000\rangle \frac{|u_1\rangle + |u_2\rangle}{\sqrt{2}} = \frac{|\phi_1\rangle |u_1\rangle + |\phi_2\rangle |u_2\rangle}{\sqrt{2}}$$

Thus, a measurement on the additional register would cancel any effect of super-position or entanglement that we can have in the cell. The phase estimation algorithm, merged with the interpretation of the operators given by Love [42], can be applied also using operators representing specific configuration on the cells.

An appendix contains supplementary information that is not an essential part of the text itself but which may be helpful in providing a more comprehensive understanding of the research problem or it is information that is too cumbersome to be included in the body of the paper.

## References

[1] Bharadwaj, S.S., Sreenivasan, K.R.: Quantum computation of fluid dynamics. arXiv preprint arXiv:2007.09147 (2020)

[2] Gaitan, F.: Finding flows of a navier-stokes fluid through quantum computing. npj Quantum Information 6(1), 61 (2020)

[3] Brassard, G., Hoyer, P., Mosca, M., Tapp, A.: Quantum amplitude amplification and estimation. Contemporary Mathematics 305, 53-74 (2002)

[4] Oz, F., San, O., Kara, K.: An efficient quantum partial differential equation solver with chebyshev points. Scientific Reports 13(1), 7767 (2023)

[5] Oz, F., Vuppala, R.K., Kara, K., Gaitan, F.: Solving burgers' equation with quantum computing. Quantum Information Processing 21, 1-13 (2022)

[6] Hatifi, M., Di Molfetta, G., Debbasch, F., Brachet, M.: Quantum walk hydrodynamics. Scientific reports 9(1), 2989 (2019)
$$ \hat{U}_m = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & e^{-i} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & e^{-2i} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & e^{-3i} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & e^{-i} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & e^{-2i} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & e^{-3i} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & e^{-4i} \end{pmatrix} \quad (17) $$

We clearly see that the eigenvalues of this operator correspond to $e^{-im(s)}$ where $m(s)$ is the mass of the state $s$ (e.g. $m(\{001\}) = 1$). This turns the problem of measuring the mass of the cell into a phase estimation problem. Using $\hat{U}_m$ as an operator for QPE the circuit in Fig.5, we get the results in Fig. 8. We see that the rows of states $|001\rangle, |100\rangle$ and $|101\rangle, |010\rangle$, and

<figure><img src="image_11.png" /><figcaption>Figure 8: On the left: for each input state (y-axis) we plot the measurement outcome probability on the mass register (x-axis) carrying out the phase estimation algorithm in Fig. 5. We can see the success of the QPE algorithm in identifying the same rows for states with the same quantities. On the right: measurement outcome frequencies depending explicitly on the mass. We can distinguish different peaks that identify different masses with different outputs in the QPE register. The probabilities reported come from the cell states {000, 001, 101, 011, 111}</figcaption></figure>

$|110\rangle, |011\rangle$ are the same. This means that the algorithm can detect the same eigenvalues, i.e. the same mass. This procedure can be seen as a possible quantum implementation of `getPhi()` and `getMass()` functions in Algorithm1, that is probabilistic but manages to avoid direct measurement of the cell, and can be applied to any quantity of interest. The best result that can be obtained is having one peak for each eigenvalue, i.e. for each
[7] Zylberman, J., Di Molfetta, G., Brachet, M., Loureiro, N.F., Debbasch, F.: Quantum simulations of hydrodynamics via the madelung transformation. Physical Review A 106(3), 032408 (2022)

[8] Steijl, R., Barakos, G.N.: Parallel evaluation of quantum algorithms for computational fluid dynamics. Computers & Fluids 173, 22–28 (2018)

[9] Demirdjian, R., Gunlycke, D., Reynolds, C.A., Doyle, J.D., Tafur, S.: Variational quantum solutions to the advection–diffusion equation for applications in fluid dynamics. Quantum Information Processing 21(9), 322 (2022)

[10] Chen, Z.-Y., Xue, C., Chen, S.-M., Lu, B.-H., Wu, Y.-C., Ding, J.-C., Huang, S.-H., Guo, G.-P.: Quantum approach to accelerate finite volume method on steady computational fluid dynamics problems. Quantum Information Processing 21(4), 137 (2022)

[11] Leong, F.Y., Ewe, W.-B., Koh, D.E.: Variational quantum evolution equation solver. Scientific Reports 12(1), 10817 (2022)

[12] Wolf-Gladrow, D.A.: Lattice-gas Cellular Automata and Lattice Boltzmann Models: an Introduction. Springer, ??? (2004)

[13] Frisch, U., Hasslacher, B., Pomeau, Y.: Lattice-gas automata for the navier-stokes equation. Physical review letters 56(14), 1505 (1986)

[14] Wolfram, S.: Cellular automaton fluids 1: Basic theory. Journal of statistical physics 45, 471–526 (1986)

[15] Higuera, F.J., Jiménez, J.: Boltzmann approach to lattice gas simulations. Europhysics letters 9(7), 663 (1989)

[16] Rivet, J.-P., Henon, M., Frisch, U., d'Humières, D.: Simulating fully three-dimensional external flow by lattice gas methods. Europhysics Letters 7(3), 231 (1988)

[17] d'Humières, D., Lallemand, P., Frisch, U.: Lattice gas models for 3d hydrodynamics. Europhysics Letters 2(4), 291 (1986)

[18] Molvig, K., Donis, P., Miller, R., Myczkowski, J., Vichniac, G.: Multi-species lattice-gas automata for realistic fluid dynamics. In: Cellular Automata and Modeling of Complex Physical Systems: Proceedings of the Winter School, Les Houches, France, February 21–28, 1989, pp. 206–231 (1989). Springer

[19] Mainster, M.A.: Cellular automata: retinal cells, circulation and patterns. Eye 6(4), 420–427 (1992)

[20] F Brandner, A., Timr, S., Melchionna, S., Derreumaux, P., Baaden, M., Sterpone, F.: Modelling lipid systems in fluid with lattice boltzmann molecular dynamics
classical value of the corresponding quantity. This optimization, left as a future perspective, can allow to calculate directly quantities that can be used for other quantum subroutines.

## 5. Application to the FHP LGCA model

### 5.1. Collision circuit

We have $v$ qubits representing the cell, these compose the *cell register*. Each collision is a rotation in the cell, as represented in Table 1. It is possible to define these rotations in terms of quantum operations, as a series of swaps gates in Figure 9.

<figure><img src="image_12.png" /><figcaption>Figure 9: Rotations with quantum circuits</figcaption></figure>

Applying the method explained in Sec.3.2, we develop an optimal overall circuit for the implementation of 0-momentum collisions in Table1, as shown in Figure 10. Each collision is implemented in two parts. We have first the verification of the collisional states that uses one conditional qubit $|b\rangle$. The second part is the collision, which is carried out as a controlled operation on the conditional qubit using an ancilla $|a\rangle$ for simulating random outcomes of $B_2$ and $B_4$.

To implement the verification of collisional states of $B_3$ we consider a logic methodology, starting from the expression that can be found in [1], adapted to our convention.

$$b = (n_0 \land n_1) \& (n_1 \land n_2) \& (n_2 \land n_3) \& (n_3 \land n_4) \& (n_4 \land n_5) \quad (18)$$

Where $\land$ is a XOR operation and $\&$ is an AND operation. We can turn this logic expression in a quantum circuit interpreting $\land$ as CNOT and $\&$ as a generalized Toffoli gate with the cell's qubits as controls and $|b\rangle$ as a target. Then we restore the original cell and we apply the C-Swap gates for a rotation of 180° with $|b\rangle$ as a control qubit.
simulations and hydrodynamics. Scientific reports 9(1), 1-14 (2019)

[21] McCullough, J., Coveney, P.: An efficient, localised approach for the simulation of elastic blood vessels using the lattice boltzmann method. Scientific Reports 11(1), 24260 (2021)

[22] Nava-Sedeño, J.M., Hatzikirou, H., Klages, R., Deutsch, A.: Cellular automaton models for time-correlated random walks: derivation and analysis. Scientific reports 7(1), 16952 (2017)

[23] Roy, I., Ray, P.K., Balasubramanian, G.: Examining oxidation in β-nial and β-nial+ hf alloys by stochastic cellular automata simulations. npj Materials Degradation 5(1), 55 (2021)

[24] Stockman, H., Stockmant, C., Carrigan, C.: Modelling viscous segregation in immiscible fluids using lattice-gas automata. Nature 348, 523-525 (1990)

[25] Appert, C., Rothman, D.H., Zaleski, S.: A liquid-gas model on a lattice. Physica D: Nonlinear Phenomena 47(1-2), 85-96 (1991)

[26] Chopard, B., Droz, M.: Cellular automata model for heat conduction in a fluid. Physics Letters A 126(8-9), 476-480 (1988)

[27] Frenkel, D., Ernst, M.: Simulation of diffusion in a two-dimensional lattice-gas cellular automaton: a test of mode-coupling theory. Physical review letters 63(20), 2165 (1989)

[28] Simons, N.R., Bridges, G., Podaima, B., Sebak, A.: Cellular automata as an environment for simulating electromagnetic phenomena. IEEE microwave and guided wave letters 4(7), 247-249 (1994)

[29] Chen, S., Chen, H., Martnez, D., Matthaeus, W.: Lattice boltzmann model for simulation of magnetohydrodynamics. Physical Review Letters 67(27), 3776 (1991)

[30] Boghosian, B.M., Levermore, C.D.: A cellular automaton for burgers' equation. Complex Systems 1(1), 17-30 (1987)

[31] Yepez, J.: Quantum lattice-gas model for computational fluid dynamics. Physical Review E 63(4), 046702 (2001)

[32] Yepez, J.: Quantum lattice-gas model for the diffusion equation. International Journal of Modern Physics C 12(09), 1285-1303 (2001)

[33] Yepez, J.: Quantum lattice-gas model for the burgers equation. Journal of Statistical Physics 107, 203-224 (2002)

[34] Yepez, J., Boghosian, B.: An efficient and accurate quantum lattice-gas model for
<figure><img src="image_13.png" /><figcaption>Figure 10: Collisional circuit for 0-momentum collisions of FHP. These are performed depending on a conditional qubit $|b\rangle$ that gets flipped if the input state is collisional, and an additional qubit $|a\rangle$ that introduces the non-deterministic character of 2- and 4-body collisions and is measured at each time-step. This circuit uses the invariance of collisional states for the merging of all the collisions in one circuit.</figcaption></figure>

To implement the verification of collisional states of $B_2$ and $B_4$, we apply a reasoning based on equivalence classes. We define an *asymmetric opposite pair* as a pair of opposite bit-velocities that differ from each other (e.g. $|100000\rangle$ has 1 opposite pairs for the asymmetry between $|n_5\rangle$ and $|n_2\rangle$, $|110000\rangle$ has 2 asymmetric opposite pairs for $|n_5\rangle,|n_4\rangle$ and $|n_2\rangle,|n_1\rangle$, $|100100\rangle$ has 0 asymmetric opposite pairs). We consider the equivalence class of states with 0 asymmetric opposite pairs represented in Table 3.

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Invariant states</td><td></td></tr><tr><td>$B_2$ collisional states</td><td></td></tr><tr><td>$B_4$ collisional states</td><td></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Table 3: States with 0 asymmetric pairs. These states are invariant under a rotation of 180°, which corresponds to a $B_3$ collision.</figcaption></figure>

The collisional states of $B_2$ and $B_4$ belong to this class, and the other states are invariant under these collisions and under $B_3$. Thus, we can target $|b\rangle$ to be $|1\rangle$ if there are no asymmetric opposite pairs. Then, with a series of C-Swaps with $|b\rangle$ as a control, we apply the first 120° rotation. We apply
the many-body schrödinger wave equation. Computer Physics Communications 146(3), 280-294 (2002)

[35] Pravia, M.A., Chen, Z., Yepez, J., Cory, D.G.: Experimental demonstration of quantum lattice gas computation. Quantum Information Processing 2, 97-116 (2003)

[36] Itani, W., Sreenivasan, K.R., Succi, S.: Quantum algorithm for lattice boltzmann (qalb) simulation of incompressible fluids with a nonlinear collision term. arXiv preprint arXiv:2304.05915 (2023)

[37] Todorova, B.N., Steijl, R.: Quantum algorithm for the collisionless boltzmann equation. Journal of Computational Physics 409, 109347 (2020)

[38] Itani, W., Succi, S.: Analysis of carleman linearization of lattice boltzmann. Fluids 7(1), 24 (2022)

[39] Budinski, L.: Quantum algorithm for the advection-diffusion equation simulated with the lattice boltzmann method. Quantum Information Processing 20(2), 57 (2021)

[40] Moore, C., Nordahl, M.G.: Predicting lattice gases is p-complete. Technical report, Citeseer (1997)

[41] Lu, X., Jiang, N., Hu, H., Ji, Z.: Quantum adder for superposition states. International Journal of Theoretical Physics 57, 2575-2584 (2018)

[42] Love, P.: On quantum extensions of hydrodynamic lattice gas automata. Condensed Matter 4(2), 48 (2019)

[43] Krol, A.M., Sarkar, A., Ashraf, I., Al-Ars, Z., Bertels, K.: Efficient decomposition of unitary matrices in quantum circuit compilers. Applied Sciences 12(2), 759 (2022)

[44] Bernardin, D.: Global invariants and equilibrium states in lattice gases. Journal of statistical physics 68, 457-495 (1992)

[45] Teixeira, C.M.: Continuum limit of lattice gas fluid dynamics. (1994)

[46] Meyer, D.A.: Quantum lattice gases and their invariants. International Journal of Modern Physics C 8(04), 717-735 (1997)

[47] Meyer, D.A.: Quantum mechanics of lattice gas automata: One-particle plane waves and potentials. Physical Review E 55(5), 5261 (1997)

[48] Schumacher, B., Werner, R.F.: Reversible quantum cellular automata. arXiv preprint quant-ph/0405174 (2004)
a controlled-H gate with $|b\rangle$ as control and $|a\rangle$, initialized to $|0\rangle$, as target. Then we apply the same controlled rotation of 120° with $|a\rangle$ as a control. In this way, the state $|0\rangle$ of $|a\rangle$ is entangled to a rotation of 120°, while the state $|1\rangle$ is entangled to a rotation of 240°, but only if $|b\rangle = |1\rangle$. Measuring $|a\rangle$ will cause a collapse of the cell state into one of the two, resembling a random extraction. The two collision can be made on the same circuit and with the same ancilla $|b\rangle$ for the invariances introduced in Sec.2. The circuit in Figure 10 was verified in Qiskit applying a measurement on the cell register. The results are reported in Fig.11.

<figure><img src="image_15.png" /><figcaption>Figure 11: Probabilities of outcome from a measurement of the cell register at the end of the collisional circuit in Fig.10. We can see that the collision operator applied is a diagonal operator except for $B_3$ collisional states (21,42), $B_2$ collisional states (9,18,36), $B_4$ collisional states (27,45,54). This gaurantees the correctness of the decomposition.</figcaption></figure>

The progress we report stands in the computational cost of this quantum procedure. Each algorithm must be decomposed in a set of universal gates to be executed on a quantum computer. The simulation tool we used, Qiskit[38], provides a decomposition method for this purpose. The algorithm we developed was decomposed by Qiskit in 291 gates. If we try to decompose directly the collision operator that should give the same result as in Fig.11, we get an order of $10^4$ gates. Our method shows an optimization over the default decomposition of Qiskit of ≈ 100 times. We also propose to use for the first time the features of the classical configurations to find optimal quantum circuits, that can be further applied to other collisions with the use of more ancillary $|b\rangle$ qubits.
## 5.2. Quantum invariants

We apply the methodology of the evolution matrix $M$ for finding the number of quantum invariants. For FHP we need first to introduce random collisions. For calculating the number of quantum invariants, we need to write the unitary operator $\hat{C}$. We cannot take it directly from Figure 10 because it involves ancillary qubits and applies a measurement, which makes the operator non-unitary. However, we can consider the analogous unitary operation that, instead of relying on an ancilla, creates a superposition of desired states. This is the collision stated in Sec.3.2, specifically in Eq.6. This collision is not *precisely* the one of FHP. However, it conserves mass and momentum, and it can still be used for seeing spurious quantum invariants. A method for calculating quantum invariants of non-unitary operator is left as a future perspective.

As we can see from Table 4 we verified that the number of quantum invariants becomes smaller if we increase the number of collisions. This was a result expected classically. Moreover, we verified the conservation of

<table><caption>Table 4: Rank of evolution matrix corresponding to different collisions introduced in the unitary. The number of quantum invariants corresponds to $4^6 - r$ where $r$ is the rank of the evolution matrix</caption><thead><tr><th>Collisions</th><th>Rank of M</th><th>Quantum invariants</th></tr></thead><tbody><tr><td>$B_3$</td><td>126</td><td>3970</td></tr><tr><td>$B_{2,4}$</td><td>488</td><td>3608</td></tr><tr><td>$B_{2,3,4}$</td><td>590</td><td>3506</td></tr></tbody></table>

mass and momentum as the linear combination of Pauli operators defined by Love [20]

$$M = Z_0 + Z_1 + Z_2 + Z_3 + Z_4 + Z_5 \qquad (19)$$

$$P_x = Z_0 - Z_3 + \frac{1}{2}(Z_1 + Z_5 - Z_2 - Z_4) \qquad (20)$$

$$P_y = \frac{\sqrt{3}}{2}(Z_1 + Z_2 - Z_4 - Z_5) \qquad (21)$$

Where $Z_i$ is the Z operator on the i-th qubit (e.g. $Z_2 = IIIIZII$, $Z_{0,3} = IIZIIIZ$). Beyond these operators, analogously to D1Q3, other operators with a classical counterpart are conserved, considering $B_{2,3,4}$, as $I^{\otimes 6}$,
$Z_{0,3}, Z_{1,4}, Z_{0,1,3,4}$ etc., as well as the symmetric respect to the $X^{\otimes 6}$ multipli- cation. All these quantities are unexpected, but we should consider that in the end we are emulating exactly the classical system with its conservation laws. Thus, we believe that the measurement process makes these additional invariants negligible, and the resulting algorithm can be used for FHP sim- ulation on a quantum computer. A more detailed analysis of the quantum invariants is out of the scope of this paper.

<figure><img src="image_16.png" /><figcaption>Figure 12: Above: for each input state (x-axis) we plot the measurement outcome proba- bility on the mass register (y-axis) carrying out the phase estimation algorithm in Fig. 5. We can see the success of the QPE algorithm in identifying the same columns for states with the same quantities. Below: measurement outcome frequencies depending explicitly on the mass. We can distinguish different peaks that identify different masses with dif- ferent outputs in the QPE register. The probabilities reported come from the cell states $\{2^i\}$ for $i \in \{0, 1, 2, 3, 4, 5\}$</figcaption></figure>

5.3. *Quantum phase estimation for quantities*

We applied the QPE procedure to FHP for mass and momentum in x and y directions. For each procedure we used the same operators as previ- ously defined, with elements $u_q(s, s') = \delta_{s,s'}e^{iq(s)}$. The results we obtained are shown in Figure 12. Analogous simulations were run for momentum in x and y direction. The success of the algorithm can be qualitatively seen
by the presence of one peak for each mass. We proved that a phase estima-
tion algorithm offers an alternative to retrieving quantities of interest. This
procedure can be equally applied to any DnQv model.

6. Conclusion

In this paper, we propose different ideas and methodologies for executing
LGCA on quantum computers. In the first place, we highlight the limita-
tions of an advantageous quantum representation of classical DnQv models
using CBE, expanding previous results. We show that the limitation is not
given specifically by the CBE adopted, but relies also on the definition of
the streaming operator. This opens perspectives on future directions to con-
sider for moving towards quantum advantage. We develop and prove the
validity of different methodologies for finding quantum collision circuits. In
particular, we develop a new collision circuit for D1Q3, proving the validity
of a method that uses QPE, and that can be expanded to any DnQv model.
Furthermore, we develop a quantum collision circuit for 0-momentum col-
lisions of FHP. The implementation of non-deterministic collisions can be
done for algorithms with CBE presenting linear or sublinear encoding of the
space, providing reinitialization in this second case. This is due to the mea-
surement of an ancillary qubit as a quantum analogy to a classical random
extraction. Our algorithm for FHP collisions proves to be optimal in the
number of universal gates needed, with a lower cost, compared to the de-
fault Qiskit decomposition, of ≈ 100 times. The method we use, consisting
of relying on classical features of the verification process, can be further ap-
plied to finding optimal circuits for executing other collisions. We develop
and prove the effectiveness of a method for counting quantum invariants,
showing unexpected results that contribute to the study of quantum DnQv
models. This method can be applied to any unitary collision operator that
can be defined for any DnQv model. A method for quantum invariants of
non-unitary operators is left as a future perspective. Finding an explanation
and the possible consequences of the presence of numerous quantum spuri-
ous invariants will be subject to further studies. In the last place, QPE was
used for the first time as a subroutine for detecting quantities of interest
using CBE. This procedure can be carried out in quantum algorithms that
use CBE, with any encoding of the space, and for any purpose that involves
moving physical information to an additional register. All the methods here
presented can be applied to quantum algorithms that involve CBE, or that
generally use an orthogonal set of quantum states for representing different classical occupation states.

## Acknowledgments

This work is supported by the PEPR integrated project EPiQ ANR-22-PETQ-0007, by the ANR JCJC DisQC ANR-22-CE47-0002-01 founded from the French National Research Agency and with the support of the french government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université - A*MIDEX AMX-21-RID-011.

## Authors' contribution

**Niccolo Fonio**: conceptualization, data curation, methodology, formal analysis, investigation. **Giuseppe di Molfetta** and **Pierre Sagaut**: conceptualization, funding acquisition, writing - review & editing. All authors contributed to write, read and approved the final manuscript.

## Declaration of competing interest

The author declare no competing interests

## References

[1] Wolf-Gladrow, D. Lattice-gas cellular automata and lattice Boltzmann models: an introduction. (Springer, 2004)

[2] Frisch, U., Hasslacher, B. & Pomeau, Y. Lattice-gas automata for the Navier-Stokes equation. *Physical Review Letters*. **56**, 1505 (1986)

[3] Wolfram, S. Cellular automaton fluids 1: Basic theory. *Journal Of Statistical Physics*. **45**, pp. 471-526 (1986)

[4] Bharadwaj, S. & Sreenivasan, K. Quantum computation of fluid dynamics. *ArXiv Preprint arXiv:2007.09147*. (2020)

[5] Succi, S., Itani, W., Sreenivasan, K. & Steijl, R. Quantum computing for fluids: Where do we stand?. *Europhysics Letters*. **144**, 10001 (2023)
[6] Yepez, J. Quantum lattice-gas model for computational fluid dynamics. *Physical Review E*. **63**, 046702 (2001)

[7] Yepez, J. Quantum lattice-gas model for the Burgers equation. *Journal Of Statistical Physics*. **107** pp. 203-224 (2002)

[8] Yepez, J. Quantum lattice-gas model for the diffusion equation. *International Journal Of Modern Physics C*. **12**, 1285-1303 (2001)

[9] Yepez, J. & Boghosian, B. An efficient and accurate quantum lattice-gas model for the many-body Schrödinger wave equation. *Computer Physics Communications*. **146**, 280-294 (2002)

[10] Meyer, D. Quantum lattice gases and their invariants. *International Journal Of Modern Physics C*. **8**, 717-735 (1997)

[11] Meyer, D. Quantum mechanics of lattice gas automata: One-particle plane waves and potentials. *Physical Review E*. **55**, 5261 (1997)

[12] Itani, W., Sreenivasan, K. & Succi, S. Quantum algorithm for lattice Boltzmann (QALB) simulation of incompressible fluids with a nonlinear collision term. *Physics Of Fluids*. **36** (2024)

[13] Itani, W. & Succi, S. Analysis of Carleman Linearization of Lattice Boltzmann. *Fluids*. **7**, 24 (2022)

[14] Budinski, L. Quantum algorithm for the advection-diffusion equation simulated with the lattice Boltzmann method. *Quantum Information Processing*. **20**, 57 (2021)

[15] Budinski, L. Quantum algorithm for the Navier Stokes equations by using the streamfunction vorticity formulation and the lattice Boltzmann method. *ArXiv Preprint arXiv:2103.03804*. (2021)

[16] Childs, A. & Wiebe, N. Hamiltonian simulation using linear combinations of unitary operations. *ArXiv Preprint arXiv:1202.5822*. (2012)

[17] Low, G. & Chuang, I. Hamiltonian simulation by qubitization. *Quantum*. **3** pp. 163 (2019)
[18] Schalkers, M. & Möller, M. On the importance of data encoding in quantum Boltzmann methods. *Quantum Information Processing*. **23**, 20 (2024)

[19] Todorova, B. & Steijl, R. Quantum algorithm for the collisionless Boltzmann equation. *Journal Of Computational Physics*. **409** pp. 109347 (2020)

[20] Love, P. On quantum extensions of hydrodynamic lattice gas automata. *Condensed Matter*. **4**, 48 (2019)

[21] Bernardin, D. Global invariants and equilibrium states in lattice gases. *Journal Of Statistical Physics*. **68** pp. 457-495 (1992)

[22] Wing-Bocanegra, A. & Venegas-Andraca, S. Circuit implementation of discrete-time quantum walks via the shunt decomposition method. *Quantum Information Processing*. **22**, 146 (2023)

[23] Shakeel, A. Efficient and scalable quantum walk algorithms via the quantum Fourier transform. *Quantum Information Processing*. **19**, 323 (2020)

[24] Razzoli, L., Cenedese, G., Bondani, M. & Benenti, G. Efficient implementation of discrete-time quantum walks in NISQ devices. *ArXiv Preprint arXiv:2402.01854*. (2024)

[25] Schumacher, B. & Werner, R. Reversible quantum cellular automata. *ArXiv Preprint quant-ph/0405174*. (2004)

[26] Farrelly, T. A review of quantum cellular automata. *Quantum*. **4** pp. 368 (2020)

[27] Zamora, A., Budinski, L., Niemimäki, O. & Lahtinen, V. Efficient quantum lattice gas automata. *Computers & Fluids*. **286** pp. 106476 (2025)

[28] Kruger, T. Lattice Boltzmann Method-Principles and Practice. (Springer International Publishing, 2016)

[29] Mohamad, A. Lattice boltzmann method. (Springer, 2011)

[30] Nielsen, M. & Chuang, I. Quantum computation and quantum information. (Cambridge University Press, 2010)
[31] Rivet, J. & Boon, J. Lattice Gas Hydrodynamics. (Cambridge University Press, 2001)

[32] Zanetti, G. Counting hydrodynamic modes in lattice gas automata models. *Physica D: Nonlinear Phenomena* **47**, 30-35 (1991)

[33] Zanetti, G. Hydrodynamics of lattice-gas automata. *Physical Review A* **40**, 1539 (1989)

[34] Watrous, J. On one-dimensional quantum cellular automata. *Proceedings Of IEEE 36th Annual Foundations Of Computer Science*. pp. 528-537 (1995)

[35] Kitaev, A. Quantum measurements and the Abelian stabilizer problem. *ArXiv Preprint quant-ph/9511026*. (1995)

[36] Kocherla, S., Song, Z., Chrit, F., Gard, B., Dumitrescu, E., Alexeev, A. & Bryngelson, S. Fully quantum algorithm for mesoscale fluid simulations with application to partial differential equations. *AVS Quantum Science*. **6** (2024)

[37] Shende, V., Bullock, S. & Markov, I. Synthesis of quantum logic circuits. *Proceedings Of The 2005 Asia And South Pacific Design Automation Conference*. pp. 272-275 (2005)

[38] Javadi-Abhari, A., Treinish, M., Krsulich, K., Wood, C., Lishman, J., Gacon, J., Martiel, S., Nation, P., Bishop, L., Cross, A. & Others Quantum computing with Qiskit. *ArXiv Preprint arXiv:2405.08810*. (2024)

[39] Doolen, G. Lattice gas methods for partial differential equations. (CRC Press, 2019)

## Appendix A. Sublinear encoding of the space

Observing the possibilities of a sublinear encoding of the space is of major interest for reducing drastically the number of qubits required, giving a practical advantage for these algorithms since they rely on large grids. Previous works [14, 15] have provided quantum algorithms for LBM simulations. Their remark is to practice the streaming step with a quantum walk procedure. We show that the same procedure applied to LGCA results in an exact
evolution of the system, but constrains the collision step, as partly treated in [18]. We show, expanding this previous result, that a unitary streaming is possible does not allow a collision step for the intrinsic indistinguishability of the cell imposed.

The most general sublinear encoding of the lattice is

$$|\Psi(t)\rangle = \frac{1}{\sqrt{N}} \sum_x |x\rangle |\psi(x, t)\rangle \qquad (A.1)$$

where N is the number of cells. The information about the cell is contained in $|\psi(x)\rangle$. We can see that each cell needs to have some fundamental features. There must be information about each velocity and information about the occupation state of the corresponding velocity. The separability of velocities allows us to use a quantum walk procedure. The separability of the occupation state is necessary for the execution of arbitrary collisions. One example is the following encoding of the cell. The first register $|v\rangle$ represents the velocity, while the second qubit represents the occupation state being V the number of velocities,

$$|\psi(x, t)\rangle = \frac{1}{\sqrt{V}} \sum_v |v\rangle |n_v(x, t)\rangle \qquad (A.2)$$

With this encoding, we have a correspondence between any classical state of the cell and the quantum counterpart, as we can see in this example of D1Q2 model, equal to D1Q3 without rest particle.

$$[00] \rightarrow |\psi_{00}\rangle = \frac{|00\rangle + |10\rangle}{\sqrt{2}} = |+0\rangle \qquad (A.3)$$

$$[01] \rightarrow |\psi_{01}\rangle = \frac{|00\rangle + |11\rangle}{\sqrt{2}} = |\beta_{00}\rangle \qquad (A.4)$$

$$[10] \rightarrow |\psi_{10}\rangle = \frac{|01\rangle + |10\rangle}{\sqrt{2}} = |\beta_{01}\rangle \qquad (A.5)$$

$$[11] \rightarrow |\psi_{11}\rangle = \frac{|01\rangle + |11\rangle}{\sqrt{2}} = |+1\rangle \qquad (A.6)$$

The direct advantage of an encoding like this is the streaming step, which can be translated in a series of controlled shift operators, implementing a quantum walk with the following operator
$$ \hat{S} = \sum_{v} \hat{\Delta}_{v} \otimes |v\rangle \langle v| \otimes I = \sum_{v} \hat{\Delta}_{v} \otimes \hat{I}_{v} \qquad (A.7) $$

Where $\hat{I}_v = |v\rangle \langle v| \otimes \hat{I}$, and the identity acts on the occupation register because the streaming does not change the information about the occupation register, but only its location. Also,

$$ \hat{\Delta}_v = \sum_x |x+v\rangle \langle x| $$

The streaming procedure does not apply only to the encoding encoding in Eq. A.3-A.6, which is not considered in the following discussion. If we apply $S$ to the lattice, in the cell-register we are going to have the subsystem of the cell that is divided depending on the velocity, according to the streaming operator.

$$ \begin{aligned} |\Psi(t+1)\rangle &= \hat{S} |\Psi(t)\rangle \\ &= \left(\sum_v \hat{\Delta}_v \otimes \hat{I}_v\right) \frac{1}{\sqrt{N}} \sum_x |x\rangle \langle\psi(x,t)\rangle \\ &= \frac{1}{\sqrt{N}} \sum_x \sum_v \hat{\Delta}_v |x\rangle \otimes \hat{I}_v |\psi(x,t)\rangle \\ &= \frac{1}{\sqrt{N}} \sum_x \sum_v |x+v\rangle \otimes \hat{I}_v |\psi(x,t)\rangle \end{aligned} $$

In general, we can write

$$ \begin{aligned} |\Psi(t+1)\rangle &= \sum_x |x\rangle \otimes \sum_v \hat{I}_v |\psi(x-v, t)\rangle \\ &= \sum_x |x\rangle \otimes |\psi(x, t+1)\rangle \end{aligned} $$

Where

$$ |\psi(x, t+1)\rangle = \sum_v \hat{I}_v |\psi(x-v, t)\rangle \qquad (A.8) $$

The Eq.A.8 depends directly on the choice of performing a quantum walk for moving the information through the lattice. It says that the information after the streaming (lhs) merges information before the streaming (rhs). We
can show that this condition, which holds for any DnQv model with the general sublinear encoding presented, forbids the post-streaming state to belong to an orthogonal set, thus forbidding a unitary collision. To prove this, we consider the D1Q2 case, but the same procedure can be extended to other DnQv models

We define the most general encoding as follows, in matrix notation for simplicity, we consider a cell register of 2 qubits

$$|\psi_{i,j}\rangle = \begin{bmatrix} a_{i,j} \\ b_{i,j} \\ c_{i,j} \\ d_{i,j} \end{bmatrix} \qquad (A.9)$$

Using $|v_0\rangle = |0\rangle$, $|v_1\rangle = |1\rangle$ and we try to solve Eq.A.8. We can formulate Eq.A.8 for D1Q2 as follows $\forall k, k' \in 0, 1, 2, 3$

$$|\psi_{ij}\rangle = \hat{I}_0 |\psi_{ik}\rangle + \hat{I}_1 |\psi_{k'j}\rangle \qquad (A.10)$$

This is a set of 16 equations. Some of them are trivial and if we use Eq.A.9, we find that there are 8 conditions to satisfy

$$ \begin{array}{l} c_{00} = c_{10} \qquad d_{00} = d_{10} \qquad a_{00} = a_{10} \qquad b_{00} = b_{10} \\ c_{11} = c_{01} \qquad d_{11} = d_{01} \qquad a_{11} = a_{10} \qquad b_{11} = b_{10} \end{array} \qquad (A.11) $$

We can rewrite our states as follows

$$ \begin{align*} |\psi_{00}\rangle &= \begin{pmatrix} a \\ b \\ c \\ d \end{pmatrix} & |\psi_{01}\rangle &= \begin{pmatrix} a \\ b \\ e \\ f \end{pmatrix} \\ |\psi_{10}\rangle &= \begin{pmatrix} g \\ h \\ c \\ d \end{pmatrix} & |\psi_{11}\rangle &= \begin{pmatrix} g \\ h \\ e \\ f \end{pmatrix} \end{align*} $$

Now we can apply the orthogonality condition and we get these 6 equations

$$ \langle\psi_{i,j}|\psi_{i',j'}\rangle = \delta_{i,i'}\delta_{j,j'} \qquad (A.12) $$
We apply the following definitions to rewrite A.12 and normalization conditions

$$A = |a|^2 + |b|^2$$

$$B = |c|^2 + |d|^2$$

$$C = |e|^2 + |f|^2$$

$$D = |g|^2 + |h|^2$$

$$E = a^*g + b^*h$$

$$F = c^*e + d^*f$$

And we get the following set of equations

$$A + F = 0$$

$$E + B = 0$$

$$E + F = 0$$

$$E + F^* = 0$$

$$E + C = 0$$

$$D + F = 0$$

$$A + B = 1$$

$$A + C = 1$$

$$B + D = 1$$

$$C + D = 1$$

The first two lines are the orthogonality conditions, and the last line are the normalization conditions. It is easy now to see that these equations are not solvable, since according to the first three $A+B=0$, while the normalization condition imposes $A+B=1$. The peculiar aspect is that even increasing the space, i.e. increasing the qubits for an encoding, brings to the same set of equations. This is intrinsically caused by the streaming operator that separates the information, making it impossible to merge it in a distinguishable way. These considerations are intended as a hint that a sublinear encoding preserving orthogonality in the cell register and performing the evolution of LGCA through quantum walk is not possible.

In conclusion, this encoding does not permit the practice of the collision step, for the impossibility of ensuring the orthogonality of the cell states. Our results supposed for simplicity to have the case of the chosen velocities $|0\rangle$ and $|1\rangle$, and of the streaming operator Eq.A.7. The same proof holds in case of different velocities, that must anyway always be orthogonal. This orthogonality of the velocities, needed for the streaming procedure, brings to similar equations to EqsA.11 for any DnQv model. A general proof considering an alternative and more general streaming is yet to be found. Alternative streaming would change the character of the algorithm, which is conceived as a naive translation of the classical algorithm using a superposition.

The circuit that performs the streaming procedure we just explained, considering the example encoding in Eq. A.3-A.6, is given in Fig.A.13. Here
the register $|x\rangle$ is the space register. The cell occupational state is represented with two qubits $|v\rangle$ and $|n\rangle$. $|v\rangle$ has the information about the velocity ($|0\rangle \rightarrow v = +1$ and $|1\rangle \rightarrow v = -1$) and will be used by the streaming operator. $|n\rangle$ has information about the occupational state of the corresponding velocity. The classical-quantum correspondence is in Eq. A.3-A.6. First, the space register undergoes a series of Hadamard gates, creating the superposition. Then the operation $A_0$ is the initialization procedure, which is a series of generalized Toffoli gates with space register qubits and $|v\rangle$ as controls, and $|n\rangle$ as target. More optimal operations could be found [37]. Then, the operators $\Delta_r$ and $\Delta_l$ changes the position according to the velocity state of $|v\rangle$, which is the control. This collisionless multi-time step algorithm was run for 24 time steps and the results are shown in Fig.A.14.

<figure><img src="image_17.png" /><figcaption>Figure A.13: Quantum circuit with unitary steaming for D1Q2. The *x* register is the space register, *v* is the velocity register and *n* is the occupation register. Here we count $2^4 = 16$ sites. The $A_0$ operation is a series of generalized Toffoli that initialize the occupation qubits in the *n* register. The controlled $\Delta_v$ operations are repeated for each time step</figcaption></figure>

## Appendix B. Complete evolution of D1Q3 operators

In TableB.5 we give all the evolved operators as calculated with the collision provided by Love. Thus, it is possible to verify directly the conservation of the quantities reported in the article. Analogous evolved operator for FHP can be provided upon request
<figure><img src="image_18.png" /><figcaption>Figure A.14: D1Q2 lattice mass density at different time-steps with 64 cells run with the quantum circuit in Fig.A.13. This simulation was run with mps simulator by Qiskit. At increasing time steps we can see the system diffuses and particles move according to respective velocities. After each value of *t* the system was measured with 1000 shots.</figcaption></figure>
<table><caption>Table B.5: Table of the evolution of all the D1Q3 Pauli operators</caption><thead><tr><th>Input state <span style="text-decoration: overline;">O</span></th><th>Output state <span style="text-decoration: overline;">C</span><sup>†</sup><span style="text-decoration: overline;">O</span><span style="text-decoration: overline;">C</span></th><th>Input state <span style="text-decoration: overline;">O</span></th><th>Output state <span style="text-decoration: overline;">C</span><sup>†</sup><span style="text-decoration: overline;">O</span><span style="text-decoration: overline;">C</span></th></tr></thead><tbody><tr><td>III</td><td>III</td><td>XZZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXX - IYY + XII + XZZ)</td></tr><tr><td>IIX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIX + XXI + YYI + ZZX)</td><td>YII</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXY + IYX + YII + YZZ)</td></tr><tr><td>IYI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYI + XYI - YXI + ZZY)</td><td>YIX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYI - XIY + YIX + ZYZ)</td></tr><tr><td>IIZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIZ + IZI - ZII + ZZZ)</td><td>YIY</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXI + XIX + YIY - ZXZ)</td></tr><tr><td>IXI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXI + XIX - YIY - ZXZ)</td><td>YIZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(YIZ + YZI + ZXY - ZYX)</td></tr><tr><td>IXX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXX - IYY + XII - ZZZ)</td><td>YXI</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IYI + XYI + YXI + ZZY)</td></tr><tr><td>IXY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXY + IYX - YII + YZZ)</td><td>YXX</td><td><sup>1</sup>&frasl;<sub>2</sub>(-XXY + XYX + YXX - YYYY)</td></tr><tr><td>IXZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXZ + XZX - YZY - ZXI)</td><td>YXY</td><td>YXY</td></tr><tr><td>IYI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYI + XIY + YIX - ZYZ)</td><td>YXZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IZY + XYZ + YXZ + ZIY)</td></tr><tr><td>IYX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYX + IYX + YII - YZZ)</td><td>YYI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIX - XXI + YYI - ZZX)</td></tr><tr><td>IYY</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXX + IYY + XII - XZZ)</td><td>YYX</td><td>YYX</td></tr><tr><td>IYZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYZ + XZY + YZX - ZYI)</td><td>YYY</td><td><sup>1</sup>&frasl;<sub>2</sub>(-XXY + XYX - YXX + YYYY)</td></tr><tr><td>IZI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIZ + IZI + ZII - ZZZ)</td><td>YYZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZX - XXZ + YYZ - ZIX)</td></tr><tr><td>IZX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZX + XXZ + YYZ + ZIX)</td><td>YZI</td><td><sup>1</sup>&frasl;<sub>2</sub>(YIZ + YZI - ZXY + ZYX)</td></tr><tr><td>IZY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZY + XYZ - YXZ + ZIY)</td><td>YZX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYZ - XZY + YZX + ZYI)</td></tr><tr><td>IZZ</td><td>IZZ</td><td>YZY</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXZ + XZX + YZY - ZXI)</td></tr><tr><td>XII</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXX + IYY + XII + XZZ)</td><td>YZZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYX - IYX + YII + YZZ)</td></tr><tr><td>XIX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXI + XIX + YIY + ZXZ)</td><td>ZII</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IIZ + IZI + ZII + ZZZ)</td></tr><tr><td>XIY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYI + XIY - YIX + ZYZ)</td><td>ZIX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZX - XXZ - YYZ + ZIX)</td></tr><tr><td>XIZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(XIZ + ZXI - ZXX - ZYY)</td><td>ZIY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZY - XYZ + YXZ + ZIY)</td></tr><tr><td>XXI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIX + XXI - YYI - ZZX)</td><td>ZIZ</td><td>ZIZ</td></tr><tr><td>XXX</td><td>XXX</td><td>ZXI</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXZ + XZX - YZY + ZXI)</td></tr><tr><td>XXY</td><td><sup>1</sup>&frasl;<sub>2</sub>(XXY + XYX - XXX - YYYY)</td><td>ZXX</td><td><sup>1</sup>&frasl;<sub>2</sub>(-XIZ + XZI + ZXX - ZYY)</td></tr><tr><td>XXZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZX + XXZ - YYZ - ZIX)</td><td>ZXY</td><td><sup>1</sup>&frasl;<sub>2</sub>(YIZ - YZI + ZXY + ZYX)</td></tr><tr><td>XYI</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYI + XYI + YXI - ZZY)</td><td>ZXZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXI + XIX - YIY + ZZX)</td></tr><tr><td>XYX</td><td><sup>1</sup>&frasl;<sub>2</sub>(XXY + XYX + YXX + YYYY)</td><td>ZYI</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IYZ + XZY + YZX + ZYI)</td></tr><tr><td>XYY</td><td>XYY</td><td>ZYX</td><td><sup>1</sup>&frasl;<sub>2</sub>(-YIZ + YZI + ZXY + ZYX)</td></tr><tr><td>XYZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(IZY + XYZ + YXZ - ZIY)</td><td>ZYY</td><td><sup>1</sup>&frasl;<sub>2</sub>(-XIZ + XZI - ZXX + ZYY)</td></tr><tr><td>XZI</td><td><sup>1</sup>&frasl;<sub>2</sub>(XIZ + XZI + ZXX + ZYY)</td><td>ZYZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IYI + XIY + YIX + ZYZ)</td></tr><tr><td>XZX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IXZ + XZX + YZY + ZXI)</td><td>ZZI</td><td>ZZI</td></tr><tr><td>XZY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IYZ + XZY - YZX + ZYI)</td><td>ZZX</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIX - XXI - YYI + ZZX)</td></tr><tr><td>XZZ</td><td><sup>1</sup>&frasl;<sub>2</sub>(-IXX - IYY + XII + XZZ)</td><td>ZZY</td><td><sup>1</sup>&frasl;<sub>2</sub>(IIY - XYI + YXI + ZZY)</td></tr></tbody></table>