# Two-dimensional quantum lattice gas algorithm for anisotropic Burger-like equations

Niccoló Fonio, Pierre Sagaut, Giuseppe Di Molfetta

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

## To cite this version:

Niccoló Fonio, Pierre Sagaut, Giuseppe Di Molfetta. Two-dimensional quantum lattice gas algorithm for anisotropic Burger-like equations. 2026. hal-05523609

HAL Id: hal-05523609

https://hal.science/hal-05523609v1

Preprint submitted on 23 Feb 2026

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 CC BY 4.0 - Attribution - International License
Two-dimensional quantum lattice gas algorithm for anisotropic Burger-like equations

Niccolo Fonio*

Aix-Marseille Université, CNRS, Centrale Méditerranée, M2P2, Marseille, France and
Aix-Marseille Université, CNRS, LIS, Marseille, France

Pierre Sagaut

Aix-Marseille Université, CNRS, Centrale Méditerranée, M2P2, Marseille, France

Giuseppe Di Molfetta

Aix Marseille Univ, CNRS, LIS, Marseille, France

(Dated: February 20, 2026)

Building on hybrid quantum lattice gas algorithm, we revisit the possibilities of this quantum
lattice model. By deriving a correction to the predicted viscosity, we provide analytical and nu-
merical results that refine original formulation. We introduce a minimal 2D generalization of the
algorithm, which allows to simulate anisotropic Burgers-like equations while retaining only two lat-
tice velocities. This approach opens a promising route toward embedding momentum conservation
and advancing toward Navier-Stokes dynamics in 2D, going beyond Frisch, Hasslacher and Pomeau
(FHP) with a quantum native model.

# I. INTRODUCTION

The application of quantum computing (QC) to com-
putational fluid dynamics (CFD) seeks to uncover po-
tential advantages over classical methods. Several paths
have been explored to this end. Quantum lattice gas
(QLG) models share a common structure: a discrete lat-
tice in which, at each site, one can define particles (bi-
nary lattice gas, BLG, or integer lattice gas, ILG) or a
probability density function (lattice Boltzmann method,
LBM). These undergo a local collision step followed by
streaming along a discrete set of velocities. An example
of the one-dimensional evolution is depicted in Fig.1. All
these methods are capable of capturing nonlinear behav-
ior. In particular, BLG [1] and LBM [2, 3] reproduce the
Navier-Stokes equations.

There are various ways to translate a classical lattice
gas model to a quantum computer. One approach as-
signs one qubit per bit, establishing a direct correspon-
dence between classical and quantum cellular automata.
This computational basis encoding (CBE) and its limi-
tations for a quantum advantage have been studied in
some works [4–6]. Alternatively, the LBM probability
distributions can be encoded into the amplitudes of a
superposed quantum state. This amplitude encoding is
widely used [7–10], not only for LBM but also for other
lattice gas models [11–13]. Its main limitation arises from
the non-unitarity of the collision operator, which neces-
sitates multi-time-step implementations and complicates
the estimation of quantum amplitudes. These issues also
affect linearized QLBM schemes based on Carleman’s
linearization[14–16], whose possibility of simulating tur-
bulence has been analyzed in [17]. Multi-time-step im-
plementations have been demonstrated with some success

<figure><img src="image_3.png" /><figcaption>FIG. 1. Evolution of a one-dimensional lattice gas cellular automaton with three discrete velocities (D1Q3). Each cell is represented by three bits $[b_2b_1b_0]$, indicating the presence of particles moving with velocities $[-1, 0, +1]$. The collision rule that conserves both mass and momentum (assuming the rest particle has mass 2) is $[101] \leftrightarrow [010]$, which occurs at positions $x = 0, 1, 3$, and 5. All other sites remain unaffected by the collision. During the streaming step, particles propagate according to their respective velocities under periodic boundary conditions.</figcaption></figure>

for advection–diffusion and other linearized models[9, 18–20]. Alternative encodings, such as space–time formulations, are also being explored [6, 21–24].

In all these approaches, the primary objective is to
achieve a more efficient or expressive representation of
the classical lattice. From this point of view, the QLG
model proposed by Yepez in [25] represents a distinct and
conceptually complete alternative, being one of the ear-
liest unitary QLG formulations capable of reproducing
hydrodynamic behavior. Rather than encoding a clas-
sical lattice model, it directly defines a programmable
quantum lattice gas, ensuring strict unitarity of the col-
lision operator by construction. This property makes it
particularly appealing for genuine quantum simulation,
as it embeds dissipative macroscopic phenomena,such as

<footnote>* niccol.fonio11@gmail.com</footnote>
viscosity, within a reversible collisional microscopic dynamics.

At the same time, this model raises several open questions. Its classical simulability, as discussed in [26], challenges the extent to which it provides a genuine quantum advantage, while its structural simplicity and numerical stability make it an attractive lattice scheme in its own right. Nevertheless, classical simulability can still be reconsidered in light of later analyses [27], particularly when including a fully quantum streaming operator, for which an analytical treatment remains unavailable. This issue is of special relevance for quantum simulations of fluid systems [28, 29].

Revisiting and extending this model, especially to two spatial dimensions, is thus of both foundational and practical interest. On one hand, it deepens our understanding of how hydrodynamic behavior emerges from a unitary, quantum-compatible lattice dynamics; on the other, it opens the possibility of constructing reduced, resource-efficient models capable of capturing 2D nonlinear dynamics. In this sense, the quantum lattice model we focus on serves not only as a testing ground for exploring the interface between quantum computation and fluid dynamics, but also as a promising pathway toward alternative mesoscopic formulations with tunable hydrodynamic parameters and lower computational costs than classical lattice gas automata, such as the Frisch-Hasslacher-Pomeau (FHP) model.

In summary, the QLG we analyze stands out as a quantum-native numerical method for simulating nonlinear dynamics. Its possible fully quantum formulation remains unexplored, while its degree of classical simulability is still open to discussion. In this sense, the model could either enable simulations beyond classical computational capabilities, or offer classical simulations with features unattainable in existing lattice schemes, such as tunable viscosity in FHP-like models, or lower viscosities and unconditional stability compared to LBM. These considerations motivate the present work, which has two principal outcomes. In Sec.II, we derive a correction to the viscosity, providing a detailed analytical procedure elaborated in the Appendix. In Sec.III, we obtain the general partial differential equation simulated by a minimal two-dimensional extension of the model, leading to anisotropic Burgers-like dynamics. Finally, Sec.IV presents the corresponding numerical results.

## II. Q-D1Q2

A Q-DnQv model is a quantum lattice gas characterized by a lattice in n dimensions where at each site there are v qubits. We assume, for the sake of the algorithmic procedure we will follow, that the qubits are initially disentangled. The k-th qubit at site x and time t is in the

state $|q_k(x, t)\rangle$. The cell is represented by

$$ egin{aligned} |\psi(x,t)\rangle &= \bigotimes_{k=0}^{v-1} |q_k(x,t)\rangle \\ &= |q_0(x,t)\rangle \otimes \cdots \otimes |q_{v-1}(x,t)\rangle, \end{aligned} \quad (1) $$

and the state of the lattice at time t is

$$ |\Psi(t)\rangle = \bigotimes_x |\psi(x,t)\rangle. \quad (2) $$

We first introduce the 1D model of [25], since it will be the starting point for our considerations and results. We write the most important parts of the derivation, which is shown in its integrity in the Appendix.

The evolution of the system is composed of a collision, the measurement of observables of interest, and a reinitialization that carries out the classical streaming, as represented in Fig.2. The algorithmic procedure can also be written as in Alg.1. The qubits are initialized as follows

$$ |q_v(x,t)\rangle = \sqrt{1-f_v(x,t)}|0\rangle + \sqrt{f_v(x,t)}|1\rangle, \quad (3) $$

where $f_i(x,t)$ are the classical probability distribution functions of the i-th velocity. Thus, the state of each cell is

$$ egin{aligned} |\psi(x,t)\rangle = &\sqrt{(1-f_0(x,t))(1-f_1(x,t))}|00\rangle \\ &+ \sqrt{(1-f_0(x,t)f_1(x,t)}|01\rangle \\ &+ \sqrt{f_0(x,t)(1-f_1(x,t))}|10\rangle \\ &+ \sqrt{f_0(x,t)f_1(x,t)}|11\rangle \end{aligned} \quad (4) $$

The populations $f_i$ are modified by the collision, thus they need to be measured to be updated for streaming. They result from the measurement of number operators

$$ \hat{n}_0 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}, \quad \hat{n}_1 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}. \quad (5) $$

having

$$ f_i(x,t) = \langle \psi(x,t) | \hat{n}_i | \psi(x,t) \rangle. \quad (6) $$

It is straightforward to verify the consistency of Eq. 6 with Eq. 4.

<table><caption>**Algorithm 1 QLGCA for Burger's equation**</caption><tr><td>c_lattice ← initialize lattice</td><td>▷ f<sub>i</sub>(x, t)</td></tr><tr><td colspan="2"><b>for</b> t do</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;<b>for</b> x ∈ c_lattice do</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;|ψ(x, t)⟩ ← initialize(f<sub>0</sub>(x, t), f<sub>1</sub>(x, t))▷ from Eq. 4</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;|ψ(x, t)⟩ ← Û|ψ(x, t)⟩▷ Collision</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;f'<sub>i</sub>(x, t) ← &lt;ψ(x, t)|Û† ˆn<sub>i</sub>Û|ψ(x, t)⟩▷ Measurement</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;<b>end for</b></td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;<b>for</b> x ∈ c_lattice do</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;f<sub>i</sub>(x, t) ← f'<sub>i</sub>(x - c<sub>i</sub>, t)▷ Streaming</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;<b>end for</b></td></tr><tr><td colspan="2"><b>end for</b></td></tr></table>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Algorithmic scheme for Q-D1Q2 model</figcaption></figure>

After the initialization, each cell undergoes a collisional operation $\hat{C}$ that must preserve the mass density of the cell, defined as $\rho(x,t) = f_0(x,t) + f_1(x,t)$. The most general unitary that ensures this conservation can be parametrized as follows

$$ \hat{C} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & e^{i\xi}\cos\theta & e^{i\zeta}\sin\theta & 0 \\ 0 & -e^{-i\zeta}\sin\theta & e^{-i\xi}\cos\theta & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}, \quad (7) $$

with $\theta, \zeta, \xi$ being real angles, commonly said Euler's angles. After the collision, $\hat{n}_i$ are measured and give the post-collisional populations $f'_i(x,t)$. Considering also the streaming, which is classical, we know that $f_i(x - c_i \delta x, t + \delta t) = f'_i(x, t)$, being $\delta x$ the lattice spacing, $\delta t$ the time step, $c_i$ the corresponding velocity. We can then write the quantum lattice Boltzmann equation

$$ f_i(x - c_i \delta x, t + \delta t) = f_i(x, t) + \Omega_i(x, t), \quad (8) $$

where $\Omega_i$ is the collision term. $\Omega_i$ in general depends on $f_i(x, t)$. The equilibrium condition is the equation $\Omega_i = 0$, which allows us to calculate $f_i^{eq}$. In quantum terms, we write

$$ \Omega[\psi(x,t)] = \langle \psi(x,t) | \hat{C}^{\dagger} \hat{n}_i \hat{C} - \hat{n}_i | \psi(x,t) \rangle. \quad (9) $$

For this model, given by the mass conservation imposed by the collision, we have $-\Omega_0 = \Omega_1 = \Omega$, where

$$ \Omega = (f_0 - f_1) \sin^2 \theta + \sin(2\theta) \cos(\zeta - \xi) \sqrt{f_0(1 - f_0) f_1(1 - f_1)}. \quad (10) $$

Solving the equilibrium condition brings us to

$$ \begin{aligned} f_0^{eq} &= \frac{\rho}{2} + \frac{1}{2\alpha} \left[ \sqrt{\alpha^2 + 1} - \sqrt{(\alpha^2 + 1) - 2\alpha^2\rho + \alpha^2\rho^2} \right] \\ f_1^{eq} &= \frac{\rho}{2} - \frac{1}{2\alpha} \left[ \sqrt{\alpha^2 + 1} - \sqrt{(\alpha^2 + 1) - 2\alpha^2\rho + \alpha^2\rho^2} \right], \end{aligned} \quad (11) $$

where $\alpha = \cot \theta \cos(\zeta - \xi)$. It is very convenient to obtain the closed-form solution of this QLG model. We will use this feature for the generalization of the model

to 2D. We stress out, in particular, that the expressions in Eqs. 11 depend on the collision model and upon the fact that the populations are not mixed, thus $f_i$ streams to $f_i$. In the case of mixed streaming, e.g. $f_0$ goes to $f_1$, the equilibrium would be affected.

To obtain the PDE we are simulating with this model, we start from Eq. 8. This will also be the starting point for the 2D derivation. We first do the Taylor expansion of the finite difference up to the second order in $\delta x$, and first order in $\delta t$. This is due to the fact that we will find ourselves in the diffusive regime where $\delta x \approx \epsilon$ and $\delta t \approx \epsilon^2$ where $\epsilon$ is the Knudsen number. After the Taylor expansion, we do the Chapman-Enskog expansion. At the first-order in $\epsilon$, we can solve the corresponding equation to find $f_i^{(1)}$ as an expression of $f_i^{eq}$. In 1D we obtain

$$ \epsilon f_i^{(1)} = - \frac{1}{J_1 - J_0} c_i \delta x \partial_x f_i^{eq} \quad (12) $$

with $J_i = \frac{\partial \Omega}{\partial f_i}|_{eq}$. In the appendix, we show explicitly how to calculate $J_1 - J_0$, which results in

$$ J_1 - J_0 = -2 \sin^2(\theta) \left[ 1 - \alpha^2(1 - \rho) - \alpha^2 \frac{\rho^2 - u^2}{2} \right]. \quad (13) $$

Including also the second order, we obtain a PDE for each $f_i$ as follows

$$ \delta t \partial_t f_i^{eq} - c_i \delta x \partial_x (f_i^{eq} + \epsilon f_i^{(1)}) + \frac{\delta x^2}{2} \partial_{xx} f_i^{eq} + O(\epsilon^3) = \Omega_i. \quad (14) $$

This allows us to find a clear PDE for the mass density summing Eq. 14 for $(i = 0)$ and $(i = 1)$. The collision term in the RHS vanishes due to the conservation of mass. We obtain

$$ \begin{aligned} \delta t \partial_t \rho - \delta x \partial_x [f_1^{eq} - f_0^{eq} + \epsilon(f_1^{(1)} - f_0^{(1)})] + \\ + \frac{\delta x^2}{2} \partial_{xx} \rho + O(\epsilon^3) = 0. \end{aligned} \quad (15) $$

Defining $u = f_1^{eq} - f_0^{eq}$ and remembering Eq. 12, we have

$$ \begin{aligned} \delta t \partial_t \rho - \delta x \partial_x u + \delta x^2 (\partial_x \rho) \partial_x \left( \frac{1}{J_1 - J_0} \right) + \\ + \frac{\delta x^2}{2} \left( 1 + \frac{2}{J_1 - J_0} \right) \partial_{xx} \rho + O(\epsilon^3) = 0. \end{aligned} \quad (16) $$
We can now rewrite the PDE only in terms of $\rho$, having

$$ 
\begin{aligned}
& \delta t \partial_t \rho - \delta x \frac{\alpha(\rho-1) \partial_x \rho}{\sqrt{1+\alpha^2(\rho-1)^2}} + \\
& \delta x^2 (\partial_x \rho) \frac{1}{2 \sin^2(\theta)} \frac{\alpha^2 [-\partial_x \rho + \rho \partial_x \rho - u \partial_x u]}{(1+\alpha^2(1-\rho)+\alpha^2 \frac{\rho^2-u^2}{2})^2} + \\
& + \frac{\delta x^2}{2} \left( 1 - \frac{1}{2 \sin^2(\theta)} \frac{1}{1+\alpha^2(1-\rho)+\alpha^2 \frac{\rho^2-u^2}{2}} \right) \partial_{xx} \rho \\
& + O(\epsilon^3) = 0.
\end{aligned} 
\quad (17) $$

The next step consists of applying the low-Mach approximation, i.e. considering $\rho - 1 = O(\epsilon)$. Usually the low-Mach approximation results in an incompressible regime having $\rho \approx \rho_0$. In this case it may seem artificial to consider $\rho_0 = 1$. The question may become then if this model is effective in simulating systems where $\rho_0 \neq 1$. The answer is positive, but we need to rescale the mass density with $\rho_0$, thus if we want to simulate a system with density $\rho \approx \rho_0$, we could consider simulating with our lattice gas the system where $\rho' = \rho/\rho_0$, thus normalizing the mass.

The explicit calculation for the low-Mach regime can be found in the appendix, while here we just remind the important fact that we use $O((\rho-1)) = O(\epsilon)$. This allows us to obtain the following equation

$$ 
\begin{aligned}
& \delta t \partial_t \rho + \delta x \alpha (1 - \rho) \partial_x \rho + \\
& + \frac{\delta x^2}{2} \left( 1 - \frac{1}{\sin^2(\theta) \sqrt{\alpha^2 + 1}} \right) \partial_{xx} \rho + O(\epsilon^4) = 0.
\end{aligned} 
\quad (18) $$

In [25], the same equation is obtained but with a different viscosity, in particular

$$ 
\begin{aligned}
& \delta t \partial_t \rho + \delta x \alpha (1 - \rho) \partial_x \rho \\
& - \frac{\delta x^2}{2} \cot^2(\theta) \partial_{xx} \rho + O(\epsilon^3, \alpha^2 \epsilon) = 0
\end{aligned} 
\quad (19) $$

This difference is due to the approximation for $J_1 - J_0$, requiring the correctness of the above equation only for $O(\epsilon^3, \epsilon\alpha^2)$.

Since every term is at the same order, we can consider the formal limit of Eq. 18 and write the equation followed by any one-dimensional mass-preserving quantum lattice gas as defined

$$ \partial_t \rho + c_s(1-\rho)\partial_x \rho = \nu \partial_{xx} \rho \quad (20) $$

$$ c_s = c\alpha = \frac{\delta x}{\delta t} \cot(\theta) \cos(\zeta - \xi) \quad (21) $$

$$ \nu = - \frac{\delta x^2}{\delta t} \frac{1}{2} \left( 1 - \frac{1}{\sin^2(\theta) \sqrt{\alpha^2 + 1}} \right). \quad (22) $$

We can connect this equation to Burger's equation with a change of variable. Using

$$ w = c_s(1 - \rho), \quad (23) $$

it is straightforward to see that we obtain

$$ \partial_t w + w \partial_x w = \nu \partial_{xx} w, \quad (24) $$

which is exactly the Burger's equation. For the change of variable in Eq. 23, taking the limit for $\rho - 1 \approx 0$ corresponds to taking the limit $w/c_s \approx 0$, and this is well known as the low-Mach regime. We point out that $w$ corresponds to a macroscopic velocity field, but it is not identical to the microscopic velocity field.

This result was possible thanks to a series of factors. First of all the collision allowed for calculating the exact microscopic expression of the equilibrium populations. This is not generally guaranteed, and usually in LGCA it is necessary to apply an approximation to get their expression [1] as function of mass and momentum. In the second place, we notice that the expressions we got for $J_1 - J_0$, which was fundamental, depends on the collision we apply. This means that the 2D generalization that we propose, as any model having the same collision and streaming structure can then inherit the same equilibrium distributions and the same $J_1 - J_0$.

### III. Q-D2Q2

In this section, we introduce a generalization of the presented QLG model to 2 dimensions. In particular, 2D QLGA have already been studied with FHP model in [30], where fluid dynamic simulations were carried out. We are interested, instead, in studying the behavior of QLGA on a square lattice where, at each site, there are 2 qubits, in order to obtain the 2D Burger's equation. Using only two velocities will inevitably bring us to anisotropic equations. However, this will allow for a complete analytical result without considering further approximations, extending naturally the Burger-like behavior to 2D.

For simplicity we assume that the qubits are arranged in a square grid. We will see that our analytical result can be valid on any Bravais lattice in 2D. Each pair of qubits is initialized as in the 1D model and undergo the same collision. We allow the populations to stream also in the y-direction. We keep the same number of qubits per site and collision because, as we said, a closed-form solution for the equilibrium distributions is convenient.

We start from the following expression, given by the algorithmic procedure, considering for simplicity $\delta x = \delta y = \delta s$,

$$ f_i(x - c_{i,x} \delta s, y - c_{i,y} \delta s, t + \delta t) - f_i(x, y, t) = \Omega_i. \quad (25) $$

We do the Taylor expansion, remembering the diffusive scaling $\delta s \approx \epsilon$ and $\delta t \approx \epsilon^2$, and the Chapman-Enskog expansion, obtaining
$$f_i(x - c_{i,x}\delta s, y - c_{i,y}\delta s, t + \delta t) - f_i(x, y, t) = \delta t \partial_t f_i^{eq} - \epsilon \delta s (c_{i,x}\partial_x f_i^{(1)} + c_{i,y}\partial_y f_i^{(1)}) + \\ - \delta s (c_{i,x}\partial_x f_i^{eq} + c_{i,y}\partial_y f_i^{eq}) + \frac{\delta s^2}{2} (c_{i,x}^2 \partial_{xx} f_i^{eq} + c_{i,y}^2 \partial_{yy} f_i^{eq}) + \\ + 2c_{i,x}c_{i,y}\partial_{xy}f_i^{eq}). \quad (26)$$

Eq. 26 is the LHS of Eq. 25. The RHS, since it is not affected by the streaming because it depends on $f_i(x, t)$,

is the same as in 1D. We solve the equation at the first order, which now depends on $c_i$'s, obtaining

$$c_{i,x}\partial_x f_i^{(1)} + c_{i,y}\partial_y f_i^{(1)} = -\frac{1}{J_1 - J_0}(c_{i,x}^2\partial_{xx}f_i^{eq} + 2c_{i,x}c_{i,y}\partial_{xy}f_i^{eq} + c_{i,y}^2\partial_{yy}f_i^{eq}) + \\ -(c_{i,x}\delta x f_i^{eq} + c_{i,y}\delta y f_i^{eq})\delta s(c_{i,x}\partial_x + c_{i,y}\partial_y)\frac{1}{J_1 - J_0}. \quad (27)$$

Now we substitute Eq. 27 in Eq. 26 and we sum the two equations for $(i=0)$ and $(i=1)$. We sum the two equations because we are interested in studying the behavior of the conserved quantity of interest, the mass density in this case. For this quantity, we notice as before that the RHS is 0, because we chose the collision as in 1D to have $-\Omega_0 = \Omega_1$.

We can then rewrite this expression using $f_1^{eq} = \frac{\rho+u}{2}$ and $f_0^{eq} = \frac{\rho-u}{2}$, with $u = f_1 - f_0$. We notice that $u$ in this case is just a convenient variable, it is not connected necessarily to the momentum in any direction.

We then apply the low-Mach approximation as before. Since we haven't modified the equilibrium distributions, the definition of $u$ does not change, as well as the definition of $J_1 - J_0$. The low-Mach approximation always consists of $\rho \approx 1$, and the Eq. 28 is valid in that limit

$$\partial_t \rho + \frac{c}{2}[(c_{0,x} + c_{1,x})\partial_x \rho + (c_{0,y} + c_{1,y})\partial_y \rho] + \\ + \frac{c_s}{2}(1 - \rho)[(-c_{0,x} + c_{1,x})\partial_x \rho + (-c_{0,y} + c_{1,y})\partial_y \rho] + \\ - \frac{\nu}{2}[(c_{0,x}^2 + c_{1,x}^2)\partial_{xx}\rho + (c_{0,y}^2 + c_{1,y}^2)\partial_{yy}\rho + 2(c_{0,x}c_{0,y} + c_{1,x}c_{1,y})\partial_{xy}\rho] + O(\epsilon^4) = 0. \quad (28)$$

If we define the constant vectors and anisotropic diffusion tensor as follows

$$\mathbf{a} = \frac{c}{2}(c_{0,x} + c_{1,x}, c_{0,y} + c_{1,y}), \quad (29)$$

$$\mathbf{b} = \frac{c_s}{2}(c_{0,x} - c_{1,x}, c_{0,y} - c_{1,y}), \quad (30)$$

$$D = \frac{\nu}{2} \begin{pmatrix} c_{0,x}^2 + c_{1,x}^2 & c_{0,x}c_{0,y} + c_{1,x}c_{1,y} \\ c_{0,x}c_{0,y} + c_{1,x}c_{1,y} & c_{0,y}^2 + c_{1,y}^2 \end{pmatrix}. \quad (31)$$

Then Eq. 28 can be written compactly as

$$\partial_t \rho + \mathbf{a} \cdot \nabla \rho + \mathbf{b} \cdot [(1 - \rho)\nabla \rho] - \nabla \cdot (D\nabla \rho) = 0. \quad (32)$$

or

$$\partial_t \rho + \mathbf{d} \cdot \nabla \rho - \nabla \cdot (D\nabla \rho) = 0, \quad (33)$$

with

$$\mathbf{d} = \mathbf{a} + (1 - \rho)\mathbf{b}. \quad (34)$$

This is the final expression we obtain, thus the generalization to 2D for the Q-D1Q2 model, in the limit $\rho - 1 \approx 0$. In the first place, we can see that this includes the 1D case. In fact for $\vec{c}_0 = (-1, 0)$ and $\vec{c}_1 = (1, 0)$, we obtain Eq. 28. Secondly, we can notice an additional advective term. This appears only when the velocities are not symmetric, as we would expect. We also notice that there is a cross-term proportional to $\partial_{xy}\rho$ that can disappear for specific choices of $\vec{c}_i$. The non-linear term is preserved. This disappears in the trivial case $\vec{c}_0 = \vec{c}_1$.

This expression is valid for an arbitrary discrete velocity set that can identify a Bravais lattice. Thus, we can use any 2D Bravais lattice. We consider, as an example, the different choices represented in Fig.3. We first consider a square lattice. We can have $c_{0,i} = -c_{1,i} = c_i$, obtaining

$$\partial_t \rho + c_s(1 - \rho)[2c_x \partial_x \rho + 2c_y \partial_y \rho] - \nu[2c_x^2 \partial_{xx} \rho + 2c_y^2 \partial_{yy} \rho + 4c_x c_y \partial_{xy} \rho] = 0. \quad (35)$$

However, this is just a 1-dimensional dynamics. There is no real mixing of information in different directions. Thus, the appearance of terms in $x$ and $y$ is just a matter of frame of reference choice.

We can get rid of the cross-derivatives term imposing $c_{0,x}c_{0,y} + c_{1,x}c_{1,y} = 0$. For example we could define $c_{0,x} = 1, c_{0,y} = 0, c_{1,x} = 0, c_{1,y} = -1$, thus having the population $f_0$ moving horizontally and the population $f_1$ moving vertically. We then obtain

$$\partial_t \rho + c[\partial_x \rho - \partial_y \rho] + c_s(\rho - 1)\nabla \rho = \nu \nabla^2 \rho. \quad (36)$$

This is the Burger's equation with an additional advective term. Thus, we can remove the anisotropy, but we add an advective term.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Different possible velocity sets. Continuous (dashed) arrows correspond to $f_1(f_0)$ streaming direction. On the left we have $c_{0,i} = -c_{1,i} = c_i$. At the center we have $c_{0,x} = 1, c_{0,y} = 0, c_{1,x} = 0, c_{1,y} = -1$. On the right we have a triangular lattice and $\vec{c}_0 = (-1/2, \sqrt{3}/2)$ and $\vec{c}_4 = (1/2, \sqrt{3}/2)$</figcaption></figure>

On a triangular lattice if we chose, for example, $\vec{c}_0 = (-1/2, \sqrt{3}/2)$ and $\vec{c}_1 = (1/2, \sqrt{3}/2)$, then we would obtain

$$ \partial_t \rho + \frac{c}{2} \sqrt{3} \partial_y \rho + c'_{s}(\rho - 1) 2 \partial_x \rho + \nu' [2 \partial_{xx} \rho + \frac{3}{2} \partial_{yy} \rho] = 0. \quad (37) $$

Also in this case we have an additional advective term, but only in the y-direction. We got the non-linear term only in the x-direction, and the diffusive term is anisotropic, but it does not involve any cross derivatives.

All these velocity sets show that it is possible to obtain a Burger-like nonlinear PDE for this QLG in 2 dimensions. These equations can show some additional anisotropies or advective terms depending on the discrete set of lattice velocities, and also on the lattice itself. Since

Eq. 28 is the most general result we can obtain with 2 velocities, the features we saw, like the anisotropies or the advective terms, may be overcome only by adding more qubits. This would completely change the model, since it would need to define an appropriate mass-conserving collision. A 3-velocity model may also be useful for introducing momentum conservation. This is left as a future perspective.

# IV. NUMERICAL RESULTS

## A. Viscosities in 1D

To compare the viscosities, we run simulations for different $\theta$s, extract the viscosity of this system numerically, and plot the result. The experimental viscosity is calculated as

$$ \nu_{exp}(t) = \sum_{x} \frac{1}{N_x} \frac{\rho(x, t+1) - \rho(x, t) + \alpha(\rho-1)(\rho(x+1, t) - \rho(x, t))}{\rho(x-1, t) - 2\rho(x, t) + \rho(x+1, t)}. \quad (38) $$

For validating our hypothesis numerically, it is necessary to reduce the numerical instability caused by the shock. For this, we chose a regime with $\rho_a = 0.005$. Even with this initial condition, the presence of the shock could affect the average process. Thus we follow this procedure: (1) simulate the system, (2) look at the viscosity in space for each time step, (3) filter out the points away from 1 standard deviation with respect to the space average for each time step, (4) average the filtered data over space, (5) average the filtered data over time. From Fig.4, we see that the experimental results follow the viscosity we predicted.

The simulations where carried out for $\theta \in [0.05, \frac{\pi}{2}]$, having 30 values. In fact, we see that $\alpha \propto 1/\sin(\theta)$ diverges for $\theta = 0$. Filtering was more effective for increasing angles after $\theta \approx 1$.

However, even filtered data could not explain the disagreement with the tale. For this, we need to get deeper into the inviscid regime. In fact, for $\theta \approx 0$ we have a collision with a small effect on the populations, causing them to change very slowly. For $\theta = 0$ populations are not affected by the collision, avoiding thermalization, thus requiring infinite time to align with our prediction. For $\theta = \pi/2$ we have populations exchange, which also cannot thermalize. However, data are more stable because there is no divergence of $\alpha$. Thus, the last points in the tale necessitate more time to get to the equilibrium. If we run the same simulation for $T = 2000$ time steps and using $\theta \in [1.2, \pi/2]$, we obtain the result in Fig.4.

We align more points to our prediction, reaching agreement for $\theta \approx 1.5$. This shows that eventual disagreement in the inviscid regime can be due to longer relaxation
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. Simulation viscosities calculated according to Eq. 38, initializing the system with $N_x = 64$ and $\rho_a = 0.005$, doing the average for $T = 200$ timesteps for the plot above, and $T = 2000$ for the plot below. The blue line is the viscosities in Eq. 22, while the red line is the viscosity of [25] $\cot^2(\theta)/2$</figcaption></figure>

times. The same effect can be seen considering the same simulations as Yepez, thus having $\rho_a = 0.4$.

### Arbitrarily small viscosities

In single-relaxation-time lattice Boltzmann (LBM) models, the kinematic viscosity is given by

$$ \nu_{\text{LBM}} = c_s^2 \left( \tau - \frac{1}{2} \right), \qquad (39) $$

which imposes a lower bound on the attainable viscosity due to numerical stability constraints: single-relaxation-time LBM becomes unstable as $\tau \to 0.5$ [31, 32]. Hence, the viscosity cannot be reduced arbitrarily. However, it can be controlled. As an alternative, Entropic lattice Boltzmann models provide a more stable numerical method because they satisfy an H-theorem

[33]. However, in these models, you cannot usually control the viscosity.

In contrast, the presented QLG model follows an H-theorem, as proven in [25], thus making it unconditionally stable. Additionally, the viscosity can be controlled as explained by Eq. 22

$$ \nu_{\text{QLG}} = -\frac{\delta x^2}{\delta t} \frac{1}{2} \left( 1 - \frac{1}{\sin^2(\theta)\sqrt{\alpha^2 + 1}} \right). $$

It follows that the viscosity can be made arbitrarily small. However, as discussed in the previous section, the collision strength (controlled by $\theta$) influences the relaxation dynamics. For $\theta \approx 0$ or $\theta \approx \pi/2$, populations thermalize more slowly, delaying shock formation. This causes longer simulation times for smaller or higher viscosities

To illustrate this effect, we computed the maximal steepness of the shock front by evaluating the spatial gradient of $\tilde{w} = w(x, t)/\alpha$ [Eq. 23] and tracking its peak value over a fixed number of time steps. In particular, we calculated numerically

$$ \Delta = \max_{x,t} |\tilde{w}(x,t) - \tilde{w}(x+1,t)|. \qquad (40) $$

Since the steepness is inversely related to the effective viscosity, sharper fronts correspond to smaller viscosities. The results are shown in Fig. 5.

Several trends emerge. Increasing $\theta$ initially enhances the steepness, corresponding to a decrease in viscosity. Near $\theta \approx \pi/2$, however, the steepness decreases because the shock has not yet fully developed—the relaxation is too slow. Comparing the curves at different times (blue square crosses and orange points) confirms that the front continues to sharpen with longer evolution times.

Finally, Eq. 22 indicates that the viscosity scales with both spatial and temporal discretization. Consistently, simulations with finer grids exhibit higher effective viscosities (green points). Overall, these results demonstrate that the QLG model can achieve arbitrarily small

<figure><img src="image_8.png" /><figcaption>FIG. 5. Maximal spatial gradient of $w(x, t)/\alpha$ for different parameters.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 6. Simulation with $N_x = 64$, $L_x = 2$, $\rho_a = 0.4$, $\rho_b = 1$, $\theta = \pi/3$, $\zeta = \xi = 0$ at different time steps. Comparison with the analytical solution with viscosity as in Eq. 22(blue line), the analytical solution with viscosity as in Eq. 19 (red line), and data from QLGCA simulation (blue crosses). We simulated the system with initial condition as in Eq. 41, but for visual clarity we depicted only points for $x \in [1, 2]$.</figcaption></figure>

viscosities without compromising numerical stability, at the cost of longer equilibration times or reduced spatial resolution.

## B. Analytical comparison

We can now compare the viscosities in Eq. 20 and Eq. 19 against the analytical solution. Our equation can be mapped to a Burger's equation which is analytically solvable. We use the same method as in [25]. We start from Eq. 20

$$ \partial_t \rho + c_s(1-\rho)\partial_x \rho = \nu \partial_{xx} \rho. $$

We use periodic boundary condition over the domain $x \in [0, L_x)$, having the initial condition

$$ \rho(x, 0) = \rho_b + \rho_a \cos(\beta x), \quad (41) $$

with $\beta = 2\pi/L_x$. We do apply the change of variable in Eq. 23, $w = c\alpha(1-\rho)$, thus having the initial condition

$$ w(x, 0) = c\alpha[1 - \rho_b - \rho_a \cos(\beta x)] = w_0(x). \quad (42) $$

We remember that $c = \delta x/\delta t$ and $\delta x = L_x/N_x$ and $\delta t = \delta x^2$. We can solve the Burger's equation with initial condition in Eq. 42 using the Cole-Hopf transformation

$$ w(x, t) = \bar{w} - 2\nu\partial_x \ln \psi, \quad (43) $$

with

$$ \bar{w} = \frac{1}{L_x} \int_0^{L_x} w_0(s) ds = c\alpha(1 - \rho_b). \quad (44) $$

This transformation turns the Burger equation into a heat equation for $\psi(x, t)$. The choosen initial condition
<figure><img src="image_10.png" /><figcaption>FIG. 7. The red dots represent the mean square error between the analytical solution with viscosity of Eq. 19 and the simulation data. The blue crosses represent the mean square error between the analytical solution with viscosity of Eq. 20 and the simulation data.</figcaption></figure>

makes it profitable to write the solution with Bessel functions. For general initial conditions, Fourier coefficients for the heat equation should be calculated. Thus, the solution with our initial condition is explicitly

$$ \psi(x,t) = I_0(A) + 2 \sum_{l=1}^{\infty} I_l(A) e^{-\nu l^2 \beta^2 t} F_l(x), \quad (45) $$

with

$$ F_l(x) = \cos\left(\frac{l\pi}{2}\right) \cos(l\beta x) + \sin\left(\frac{l\pi}{2}\right) \sin(l\beta x), \quad (46) $$

where $A = \frac{c\alpha\rho_a}{2\nu\beta}$. From here we can calculate $w(x,t)$ using Eq. 43 and then $\rho(x,t)$ with Eq. 23. To obtain the same condition as Yepez, we run the simulation truncating the series to 80 terms. We simulate the systems with the parameters specified in Fig.6.

We can see that the analytical solution calculated with the corrected viscosity reproduces better the quantum lattice gas. This explains also the small gap that was already reported in [25]. This comparison can be carried out by calculating the mean squared errors, having the results in Fig.7. Here we see the error decreasing because of the damping of the wave occurring after the formation of the shock.

### C. Simulations in 2D

We initialize the system with the following initial condition

$$ \rho(x_i, y_j, 0) = \rho_b + \rho_a \left[ \cos\left(\frac{2\pi i}{N_x}\right) + \cos\left(\frac{2\pi j}{N_y}\right) \right]. \quad (47) $$

In Fig.8 we can see the time evolution of the quantum lattice gas for different time steps.

In particular, we can appreciate the shock formation in each of the 4 cases. The first two are actually equivalent to 1D dynamics, while the third and fourth velocity sets actually mix the information along two orthogonal axes. We can still appreciate the shock formation with the third velocity set, as predicted in Eq. 36 while having an advection. In the fourth case, it is evident as well the shock formation and the advection, with an additional asymmetry brought by the asymmetry of the velocity set.

In this article, we do not give an explicit solution for the 2D equation we found. However, its correctness can be assessed with the comparison of the same equation solved by another numerical method. Thus, we consider Eq. 28 and solve it with an explicit Finite Difference-based numerical method. As we can see in Fig.9 there is a good qualitative and numerical agreement between the quantum lattice gas and the Finite Difference method solution. The QLG was, as expected, more stable, and the comparison with $L_2$ norm is shown only until the FDM did not explode. This usually occurred, as expected, at the formation of the shock, which can be attributed to as the cause of the instability of the FDM adopted. The agreement, however, did not motivate a further comparison with other classical numerical schemes.

## V. CONCLUSION AND PERSPECTIVES

In this work, we have analyzed the one-dimensional quantum lattice gas model proposed in [25]. In particular, we derived an analytical correction to the viscosity, which extends the validity of the Burgers-like equation to order $O(\epsilon^4)$ in the low-Mach regime. We discuss in detail the implications of this regime and provide numerical evidence supporting the correction, including comparisons with analytical solutions of the Burgers equation.

Furthermore, we highlighted the capability of this model to reproduce arbitrarily small viscosity for the Burger's equation, paying the price of longer simulation times or less grid points.

In the second part, we proposed a minimal 2D extension of the algorithm. This model retains only two discrete velocities, preserving the analytical tractability and conceptual clarity of the 1D case. We demonstrated that, under these conditions, the scheme exhibits two-dimensional anisotropic Burgers-like behavior and derived the explicit form of the corresponding partial differential equation as a function of the chosen velocity set.

The interest in this minimal 2D model arises from certain quantum walks. Specifically, quantum walks on a triangular lattice — with amplitudes defined on edges — which can reproduce 2D transport equations [34, 35]. This suggests a potential path toward a quantum lattice gas model distinct from FHP, relying on fewer qubits while still capturing Navier-Stokes hydrodynamics. The
reduced number of collisional qubits would also be better for simulations on real hardware.

In conclusion, this model is promising for CFD simulations on quantum computers, since it combines analytical simplicity with tunable hydrodynamic parameters. Its unitary structure makes it a candidate framework for implementing quantum-native fluid simulations, potentially bridging the gap between quantum algorithms and classical lattice Boltzmann schemes. Future developments include: introducing a third velocity to incorporate momentum conservation for approaching Navier-Stokes dynamics; studying possible isotropic alternatives with quantum walk-inspired streaming schemes; exploring the full quantum version of the algorithm, adopting a quantum streaming; assess the model's performance on near-

term quantum hardware. Such extensions would clarify whether this quantum lattice gas scheme can provide not only an efficient classical alternative but also a platform for realizing quantum advantages in computational fluid dynamics or quantum simulations.

ACKNOWLEDGMENTS

We are sincerely thankful to Jeffrey Yepez for the insights about the model, and for helpful discussions.

This work was partly supported by the PEPR EPiQ ANR-22-PETQ-0007, and ANR JCJC DisQC ANR-22-CE47-0002-01

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

[2] T. Kruger, *Lattice Boltzmann Method-Principles and Practice* (Springer International Publish, 2016).

[3] S. Succi, *The lattice Boltzmann equation: for fluid dynamics and beyond* (Oxford university press, 2001).

[4] N. Fonio, P. Sagaut, and G. Di Molfetta, arXiv preprint arXiv:2310.07362 (2023).

[5] A. D. B. Zamora, L. Budinski, O. Niemimäki, and V. Lahtinen, *Computers & Fluids* **286**, 106476 (2025).

[6] M. A. Schalkers and M. Möller, *Quantum Information Processing* **23**, 20 (2024).

[7] L. Budinski, *Quantum Information Processing* **20**, 57 (2021).

[8] B. Ljubomir, *International Journal of Quantum Information* **20**, 2150039 (2022).

[9] D. Wawrzyniak, J. Winter, S. Schmidt, T. Indinger, U. Schramm, C. Janßen, and N. A. Adams, arXiv preprint arXiv:2405.13391 (2024).

[10] E. D. Kumar and S. H. Frankel, *AVS Quantum Science* **7** (2025).

[11] N. Fonio, L. Budinski, V. Lahtinen, and P. Sagaut, arXiv preprint arXiv:2504.13549 (2025).

[12] A. D. B. Zamora, L. Budinski, V. Lahtinen, and P. Sagaut, arXiv preprint arXiv:2503.23750 (2025).

[13] B. Wang, Z. Meng, Y. Zhao, and Y. Yang, *npj Quantum Information* (2025).

[14] S. Succi, W. Itani, K. Sreenivasan, and R. Steijl, *Europhysics Letters* **144**, 10001 (2023).

[15] C. Sanavio and S. Succi, *AVS Quantum Science* **6** (2024).

[16] A. A. Zecchi, C. Sanavio, S. Perotto, and S. Succi, *Quantum Science and Technology* **10**, 035039 (2025).

[17] D. Lewis, S. Eidenbenz, B. Nadiga, and Y. Subaşı, *Quantum* **8**, 1509 (2024).

[18] Y. Xiao, L. Yang, C. Shu, Y. Du, H. Dong, and J. Wu, arXiv preprint arXiv:2505.10883 (2025).

[19] L. Xu, M. Li, L. Zhang, H. Sun, and J. Yao, *Physical Review E* **111**, 045305 (2025).

[20] M. Bediche, M. van Waveren, D. Ricot, and P. Sagaut, in *QUEST IS 25* (2025).

[21] M. A. Schalkers and M. Möller, *Journal of Computational Physics* **502**, 112816 (2024).

[22] M. A. Schalkers and M. Möller, *Computers & Fluids* **285**, 106453 (2024).

[23] C. A. Georgescu, M. A. Schalkers, and M. Möller, *Computer Physics Communications*, 109699 (2025).

[24] C. A. Georgescu, M. A. Schalkers, and M. Möller, *arXiv preprint arXiv:2506.12662* (2025).

[25] J. Yepez, *Physical Review A*—Atomic, Molecular, and Optical Physics **74**, 042322 (2006).

[26] P. J. Love and B. M. Boghosian, *Physica A: Statistical Mechanics and its Applications* **362**, 210 (2006).

[27] P. Love, *Condensed Matter* **4**, 48 (2019).

[28] J. Yepez, G. Vahala, and L. Vahala, *The European Physical Journal-Special Topics* **171**, 9 (2009).

[29] J. Yepez, G. Vahala, L. Vahala, and M. Soe, *Physical Review Letters* **103**, 084501 (2009).

[30] M. M. Micci and J. Yepez, *Physical Review E* **92**, 033302 (2015).

[31] G. Wissocq, P. Sagaut, and J.-F. Boussuge, *Journal of Computational Physics* **380**, 311 (2019).

[32] G. Wissocq and P. Sagaut, *Journal of Computational Physics* **450**, 110858 (2022).

[33] S. A. Hosseini, M. Atif, S. Ansumali, and I. V. Karlin, *Computers & Fluids* **259**, 105884 (2023).

[34] P. Arrighi, G. Di Molfetta, I. Márquez-Martín, and A. Pérez, *Physical Review A* **97**, 062111 (2018).

[35] G. Di Molfetta, *Quantum Walks, Limits, and Transport Equations* (IOP Publishing, 2024).
Appendix A: Full 1D calculations

We start from finite difference expression

$$f_i(x - c_i \delta x, t + \delta t) = f_i(x, t) + \Omega_i(x, t). \quad (A1)$$

We know that $\Omega_i = c_i \Omega$ with

$$\Omega = (f_0 - f_1) \sin^2 \theta + \sin(2\theta) \cos(\zeta - \xi) \sqrt{f_0(1 - f_0)f_1(1 - f_1)}. \quad (A2)$$

If we solve $\Omega = 0$, which is done in [25], it is possible to derive the explicit equilibrium populations

$$f_i^{eq} = \frac{\rho}{2} + \frac{c_i}{2\alpha} \left[ \sqrt{\alpha^2 + 1 - 2\alpha^2\rho + \alpha^2\rho^2} - \sqrt{\alpha^2 + 1} \right]. \quad (A3)$$

These bring correctly to $\rho = f_1^{eq} + f_0^{eq}$. We will also explicitly use the momentum $u = f_1^{eq} - f_0^{eq}$, explicitly

$$u = -\frac{1}{\alpha} \left( \sqrt{1 + \alpha^2} - \sqrt{1 + \alpha^2(\rho - 1)^2} \right). \quad (A4)$$

with $J_i = \frac{\partial \Omega}{\partial f_i}|_{eq}$. We know that $\Omega[f_i^{eq}] = 0$. We can focus then on the first order, where we have

$$-c_i \delta x \partial_x f_i^{eq} = \epsilon c_i [J_0 f_0^{(1)} + J_1 f_1^{(1)}]. \quad (A8)$$

Matrix notation as used in [25] is very convenient. However, we can alternatively simplify $c_i$ and $\epsilon$ and use some algebraic manipulations to obtain

$$J_0 f_0^{(1)} + J_1 f_1^{(1)} = \frac{1}{2} (f_1^{(1)} - f_0^{(1)})(J_1 - J_0). \quad (A9)$$

From this it is straightforward to have

$$\epsilon(f_1^{(1)} - f_0^{(1)}) = -\frac{1}{J_1 - J_0} \partial_x \rho. \quad (A10)$$

This is the crucial expression that needs to be calculated explicitly to get the correction to the viscosity. From trigonometric identities we have

$$\sin(2\theta) \cos(\zeta - \xi) = 2\alpha \sin^2(\theta), \quad (A11)$$

knowing $\alpha = \cot(\theta) \cos(\zeta - \xi)$. We have in particular, using $\tilde{f}_i = 1 - f_i$

$$J_0 = \sin^2(\theta) \left[ 1 + 2\alpha \left( \frac{\sqrt{f_0 \tilde{f}_0 f_1}}{2\sqrt{\tilde{f}_0}} - \frac{\sqrt{f_0 f_1 \tilde{f}_1}}{2\sqrt{\tilde{f}_0}} \right) \right], \quad (A12)$$

$$J_1 = \sin^2(\theta) \left[ -1 + 2\alpha \left( \frac{\sqrt{f_0 \tilde{f}_0 f_1}}{2\sqrt{\tilde{f}_1}} - \frac{\sqrt{f_0 \tilde{f}_0 f_1}}{2\sqrt{\tilde{f}_1}} \right) \right]. \quad (A13)$$

We take now Eq. A1 and we do a Taylor expansion. We keep in mind that later we will do a Chapman-Enskog expansion, thus we will have $\delta x \approx \epsilon$ and $\delta t \approx \epsilon^2$, where $\epsilon$ is the Knudsen number. We obtain with the Taylor expansion

$$-c_i \delta x \partial_x f_i + \delta t \partial_t f_i + \frac{\delta x^2}{2} c_i^2 \partial_{xx} f_i = c_i \Omega. \quad (A5)$$

Now we do the Chapman-Enskog expansion, that is considering

$$f_i \approx f_i^{eq} + \epsilon f_i^{(1)} + \epsilon^2 f_i^{(2)} + O(\epsilon^3). \quad (A6)$$

Thus we obtain

$$-c_i \delta x \partial_x f_i^{eq} - c_i \delta x \epsilon \partial_x f_i^{(1)} + \delta t \partial_t f_i^{eq} + \frac{\delta x^2}{2} c_i^2 \partial_{xx} f_i^{eq} + O(\epsilon^3) = c_i \Omega [f_i^{eq}] + \epsilon c_i [J_0 f_0^{(1)} + J_1 f_1^{(1)}], \quad (A7)$$

We define for simplicity $a = \sqrt{f_0 \tilde{f}_0 f_1 \tilde{f}_1}$ and calculate

$$\begin{aligned} J_1 - J_0 &= -2 \sin^2(\theta) + \alpha \sin^2(\theta) a \left[ \frac{1}{f_1} - \frac{1}{\tilde{f}_1} - \frac{1}{f_0} + \frac{1}{\tilde{f}_0} \right] \\ &= -2 \sin^2(\theta) + \alpha \sin^2(\theta) a A. \end{aligned} \quad (A14)$$

It is possible to calculate

$$A = \frac{(f_0 - f_1)}{a^2} (1 - f_0 - f_1 + 2f_0 f_1). \quad (A15)$$

Then, from the equilibrium condition it is possible to see that

$$\frac{f_0 - f_1}{a} = -2\alpha. \quad (A16)$$

Thus using the definitions of $\rho$ and $u$, we have

$$J_1 - J_0 = -2 \sin^2(\theta) - 2\alpha^2 \sin^2(\theta) \left[ 1 - \rho + \frac{\rho^2 - u^2}{2} \right]. \quad (A17)$$

This is consistent with Yepez result, but here we managed to simplify its explicit expression, which will be needed to calculate the correction to the viscosity. It is possible to calculate this relation

$$\begin{aligned} \frac{\rho^2 - u^2}{2} = & -\frac{1}{\alpha^2} + (\rho - 1) + \\ & + \frac{1}{\alpha^2} \left[ \sqrt{\alpha^2 + 1} \sqrt{\alpha^2 + 1 - 2\alpha^2\rho + \alpha^2\rho^2} \right]. \end{aligned} \quad (A18)$$
Thus, the final expression that we aimed for is

$$J_1 - J_0 = -2 \sin^2(\theta) \left[ \begin{aligned} & \alpha^2(\rho - 1) \\ & + \sqrt{\alpha^2 + 1} \sqrt{\alpha^2 + 1 - 2\alpha^2\rho + \alpha^2\rho^2} \end{aligned} \right]. \quad (A19)$$

For further notice we calculate now these derivatives

$$\frac{\partial(J_1 - J_0)}{\partial x} = -2 \sin^2(\theta) \left[ \begin{aligned} & \alpha^2 \partial_x \rho + \\ & + \frac{\sqrt{\alpha^2 + 1}}{2} \frac{-2\alpha^2 \partial_x \rho + 2\alpha^2 \rho \partial_x \rho}{\sqrt{\alpha^2 + 1 - 2\alpha^2 \rho + \alpha^2 \rho^2}} \end{aligned} \right], \quad (A20)$$

$$\frac{\partial(J_1 - J_0)}{\partial \rho} = \sin^2(\theta) \frac{(2 - 2\rho)(\alpha^2 + 1)^{3/2}}{\sqrt{\alpha^2(\rho - 1)^2 + 1}}. \quad (A21)$$

Now we go back to Eq. A7. We need to calculate the sum of those expressions ($i = 0$) + ($i = 1$). This is because we are conserving mass. Technically it would be possible to write this expression for each conserved quantity. The right-hand side (RHS) vanishes in this sum operation, since we had $\Omega_i = c_i \Omega$ and also for any further orders in $\epsilon$. For the left-hand side(LHS) we obtain

$$\begin{aligned} \delta t \partial_t \rho - \delta x \partial_x (f_1^{eq} - f_0^{eq}) - \delta x \epsilon \partial_x (f_1^{(1)} - f_0^{(1)}) + \\ \frac{\delta x^2}{2} \partial_{xx} \rho + O(\epsilon^3) = 0. \end{aligned} \quad (A22)$$

Using Eq. A10, we rewrite

$$\delta t \partial_t \rho - \delta x \partial_x u - \delta x^2 \partial_x \rho \left[ \frac{\partial_x (J_1 - J_0)}{(J_1 - J_0)^2} \right] + \delta x^2 \frac{1}{J_1 - J_0} \partial_{xx} \rho + \frac{\delta x^2}{2} \partial_{xx} \rho + O(\epsilon^3) = 0. \quad (A23)$$

From this expression we need to calculate the low-Mach approximation, that consists in $\rho \approx 1$. Thus we need to expand the different expressions around $\rho = 1$.

Defining $g(\rho) = 1/(J_1 - J_0)$ we have

$$g(\rho) \approx g(\rho = 1) + \left. \frac{\partial g}{\partial \rho} \right|_{\rho=1} (\rho - 1) + O((\rho - 1)^2) \quad (A24)$$

$$\approx -\frac{1}{2 \sin^2(\theta)} \frac{1}{\sqrt{\alpha^2 + 1}} + \left. \frac{\partial \rho(J_1 - J_0)}{(J_1 - J_0)^2} \right|_{\rho=1} (\rho - 1). \quad (A25)$$

We know from Eq. A21 that

$$\left. \frac{\partial \rho (J_1 - J_0)}{(J_1 - J_0)^2} \right|_{\rho=1} = 0.$$

Thus we obtain, supposing $O((\rho - 1)) = O(\epsilon)$,

$$\frac{1}{J_1 - J_0} \approx -\frac{1}{2 \sin^2(\theta)} \frac{1}{\sqrt{\alpha^2 + 1}} + O(\epsilon^2). \quad (A26)$$

We can notice a term that is

$$\epsilon \delta x \partial_x \rho \frac{\partial_x (J_1 - J_0)}{(J_1 - J_0)^2}. \quad (A27)$$

From Eq. A20 we have that $\partial_x(J_1 - J_0) \propto \partial_x \rho$. Considering that in the low-Mach regime $\partial_x \rho \approx O(\epsilon)$, it is clear to see that this term can be neglected. Last but not least we consider the term $\partial_x u$.

$$\partial_x u = \frac{\alpha(\rho - 1)\partial_x \rho}{\sqrt{1 + \alpha^2(\rho - 1)^2}}. \quad (A28)$$

Now again we define $h(\rho) = \partial_x u$ and calculate the expansion around $\rho = 1$. We keep in mind, however, that we will need to go one order further, because this term is at order $\delta x$. Thus we have

$$h(\rho) \approx h(\rho = 1) + \left. \frac{\partial h}{\partial \rho} \right|_{\rho=1} (\rho - 1) + \left. \frac{1}{2} \frac{\partial^2 h}{\partial \rho^2} \right|_{\rho=1} (\rho - 1)^2 \quad (A29)$$

We can calculate these terms, having

$$h(\rho = 1) = 0, \quad (A30)$$

$$\partial_{\rho} \partial_{x} \rho \bigg|_{\rho=1} = 0, \quad (A31)$$

$$\frac{\partial h}{\partial \rho} = \frac{\alpha}{(\alpha^2(\rho - 1)^2 + 1)^{3/2}} \bigg|_{\rho=1} = \alpha, \quad (A32)$$

$$\frac{\partial^2 h}{\partial \rho^2} = \frac{6\alpha^3(\rho - 1)}{(\alpha^2(\rho - 1)^2 + 1)^{5/2}} \bigg|_{\rho=1} = 0. \quad (A33)$$

Thus, we finally get

$$\partial_x u \approx \alpha(\rho - 1) + O(\epsilon^3). \quad (A34)$$

We can finally write Eq. A23 in the low-Mach regime as

$$\partial_t \rho + c\alpha(1 - \rho)\partial_x \rho + \frac{\delta x^2}{2} \frac{1}{\sin^2(\theta)\sqrt{\alpha^2 + 1}} \partial_{xx} \rho + O(\epsilon^4) = 0. \quad (A35)$$
<figure><img src="image_11.png" /><figcaption>FIG. 8. Simulations of a 64×64 square grid for different discrete velocity sets. The blue (red) arrow represents the streaming direction of $f_0$ ($f_1$). The lattice is initialized as in Eq. 47</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 9. Comparison between QLG simulations and FDM solving Eq. 28 for different velocity sets. Simulations were made with a 64×64 square grid for different discrete velocity sets initialized as in Eq. 47. The blue (red) arrow represents the streaming direction of $f_0$ ($f_1$).</figcaption></figure>