## Decomposition of Nonlinear Collision Operator in Quantum Lattice Boltzmann Algorithm

E. Dinesh Kumar ∗ and Steven H. Frankel †

Faculty of Mechanical Engineering

Technion - Israel Institute of Technology, Haifa, 3200003, Israel. ∗ †

edk261@gmail.com , frankel@technion.ac.il

Abstract - We propose a quantum algorithm to tackle the quadratic nonlinearity in the Lattice Boltzmann (LB) collision operator. The key idea is to build the quantum gates based on the particle distribution functions (PDF) within the coherence time for qubits. Thus, both the operator and a state vector are linear functions of PDFs, and upon quantum state evolution, the resulting PDFs will have quadraticity. To this end, we decompose the collision operator for a DmQn lattice model into a product of 2( n +1) operators, where n is the number of lattice velocity directions. After decomposition, the ( n +1) operators with constant entries remain unchanged throughout the simulation, whereas the remaining ( n +1) will be built based on the statevector of the previous time step. Also, we show that such a decomposition is not unique. Compared to the second-order Carleman-linearized LB, the present approach reduces the circuit width by half and circuit depth by exponential order. The proposed algorithm has been verified through the one-dimensional flow discontinuity and two-dimensional Kolmogrov-like flow test cases.

Introduction. -The potential for simulating fluid flows on quantum computers has generated significant interest due to the unique properties of qubits, such as superposition and entanglement. By encoding flow variables onto the probability amplitudes of qubit states, we can achieve significantly reduced memory requirements and a potentially exponential increase in computing speed. However, quantum operations are inherently linear and unitary, while the equations of fluid mechanics exhibit nonlinearity. As a result, existing computational fluid dynamics (CFD) algorithms need to be redesigned and adapted to fully utilize the power of quantum processor units (QPUs).

Traditional CFD methods solve the discretized NavierStokes equation (NSE), and commonly used algorithms are implicit in time-stepping, which requires linear system solvers (LSS) [1-5]. The Lattice Boltzmann Method (LBM) has been considered an alternate approach to NSEbased CFD solvers and has been studied extensively for various flow problems. Due to locality in space and explicit time-stepping, LBM does not require LSS and has proven to be efficient on graphic processor units (GPUs) [6]. Compared to NSE, where the nonlinear term ( u · ∇ u ) is also non-local, the LB equation has nonlinearity ( u · u ) in the equilibrium function and is local. This nonlinearity does not pose any difficulty in classical computers, as u will be saved as a floating point number followed by a square operation [7,8]. However, if we mimic the classical procedure on QPUs, the quantum advantage may be lost, especially when amplitude encoding is chosen.

There are several existing works on Quantum Lattice Boltzmann (QLB) for fluid flow simulation, but most of them ignore nonlinearity and are only used for very low Reynolds number cases [9-13]. In the absence of nonlinearity, the particle distribution functions (PDF) f can be encoded directly as probability amplitudes of qubits. However, evaluating the nonlinear term involves the vector f ⊗ f . [14] highlighted that it is impossible to have a linear operator that performs the squares of probability amplitudes, as the resulting state vector may not be a unit vector. Furthermore, the no-cloning theorem states that it is impossible to copy an arbitrary quantum state to create an identical quantum state [14]. To address this, the Carleman linearization (CL) technique can be utilized which converts the finite-dimensional nonlinear system into a infinite-dimensional linear system which is then truncated for implementation [15]. Previously, CL has been successfully applied to LBM up to second-order truncation for low to moderate Reynolds number cases. When CL is applied to LBM, the resultant state vector becomes infinite, i.e., ( f f , ⊗ f f , ⊗ f ⊗ f , . . . ). In addition, another approximation has to be made for the term 1 ρ ≈ 2 -ρ , where ρ is the density, which is valid for weakly incompressible regimes. Even after successive simplifications, the Carleman LB (CLB) algorithm still requires repeated encoding and read-out of state vectors for each time step [16-19].

In the present work, instead of the CL algorithm, we propose the decomposition of the LB operator as a product of multiple operators derived from the PDFs from the previous time step. We retain the first-order Taylor ap-

proximation for 1 ρ and accept that the encoding and readout process is required for each time step [17-19]. Our approach does not involve any truncation of infinite sets. Instead, the state vector is a finite set, specifically, four sets of f appended sequentially.

In the following, we briefly describe the LB equations, followed by deriving the quadratic formulation of the LB collision operation. Later, we will present the quantum algorithm and the computational complexity of two-qubit gates.

Lattice Boltzmann Method. LBM describes the evolution of flow field through PDFs in a uniform Cartesian grid with the chosen lattice model DmQn , where m and n denotes the number of spatial dimensions and velocity directions, respectively. The single relaxation time LB model for the flow field is given by

$$f _ { i } ( \mathbf x + \mathbf e _ { i } \Delta t, t + \Delta t ) = f _ { i } ( \mathbf x, t ) - \frac { \Delta t } { \tau } \left [ f _ { i } ( \mathbf x, t ) - f _ { i } ^ { e q } ( \mathbf x, t ) \right ] \quad \text{and} \ \delta _ { j k } \\ \text{expressi} \\ \text{$\Phi$} \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi$$

where f i represents the PDF along the i th direction, e i is the lattice velocity, and the kinematic viscosity ν is related to the relaxation parameter τ = 3 ν +0 5. The equilibrium . distribution function (EDF) is given by

$$f _ { i } ^ { e q } = w _ { i } \rho \left [ 1 + \frac { e _ { i } \cdot \mathbf u } { c _ { s } ^ { 2 } } + \frac { ( e _ { i } \cdot \mathbf u ) ^ { 2 } } { 2 c _ { s } ^ { 4 } } - \frac { \mathbf u \cdot \mathbf u } { 2 c _ { s } ^ { 2 } } \right ] \quad ( 2 ) \quad \ f _ { i } ^ { * } \colon$$

where w i is the lattice weight, c s the sound speed, ρ the fluid density, and u is the flow velocity. The zeroth and first moment of PDFs yields the macroscopic flow variables density( ρ ) and momentum ( ρ u ),

$$\rho ( { \mathbf x }, t ) & = \sum _ { i } f _ { i } ( { \mathbf x }, t ) & ( 3 ) \quad \text{when} \\ o ( { \mathbf x }. t ) \mathbf u ( { \mathbf x }. t ) & = \sum \mathbf e _ { i } f _ { i } ( { \mathbf x }. t ) & ( 4 )$$

$$\rho ( x, t ) \mathbf u ( x, t ) = \sum _ { i } ^ { \frac { i } { i } } \mathbf e _ { i } f _ { i } ( x, t ) \quad \quad ( 4 ) \\.. \ \ - \. \ \. \ \. \ \. \ \.. \ \.$$

Typically, LB algorithm is split into collision and streaming steps. In the collision step, the interaction of particles and their relaxation to equilibrium state is modelled. Later the particles move to their neighbouring site according to the lattice velocity, a process termed streaming. Thus, Eq. (1) can be written as

$$f _ { i } ^ { * } ( { \mathbf x }, t ) = f _ { i } ( { \mathbf x }, t ) - \frac { \Delta t } { \tau } \left [ f _ { i } ( { \mathbf x }, t ) - f _ { i } ^ { e q } ( { \mathbf x }, t ) \right ] \quad ( 5 ) \quad \underset { \dots } { \mathbf x } \quad \text{where} }$$

$$f _ { i } ( { \mathbf x }, t + \Delta t ) & = f _ { i } ^ { * } ( { \mathbf x } - { \mathbf e } _ { i } \Delta t, t ) & ( 6 ) \underset { \text{treated} } { \text{ where } }$$

where f ∗ i is the post-collision state. In the presence of a wall boundary, we use half-way bounce back scheme to reconstruct PDFs coming from the wall node ( f -i ),

$$f _ { - i } ( { \mathbf x } _ { b }, t + \Delta t ) = f _ { i } ^ { * } ( { \mathbf x }, t ) \quad \quad ( 7 ) \ \text{ when}$$

Formulation. -

Quadratic form for collision step. First, we rewrite Eq. (2) and using Eqs. (3) and (4), we get the quadratic form for EDF as,

$$\text{sets}. \quad & \text{form for EDF as}, \\ \text{four} \\ \text{ions}, \quad & \text{f} _ { i } ^ { e q } = \frac { 1 } { \rho } w _ { i } \left [ \rho ^ { 2 } + \rho \frac { e _ { i } \cdot \rho \mathbf u } { c _ { s } ^ { 2 } } + \frac { ( e _ { i } \cdot \rho \mathbf u ) ^ { 2 } } { 2 c _ { s } ^ { 4 } } - \frac { \rho \mathbf u \cdot \rho \mathbf u } { 2 c _ { s } ^ { 2 } } \right ] \quad ( 8 ) \\ \text{e LB} \\ \text{return} \\ \text{public} \quad & \text{$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$\quad$
\quad \cdot \quad \cdot \quad \sim \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot$$

where the coefficient α ijk is given by

$$\mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } \\ \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } \\ \mathfrak { d } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s }& \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s } & \mathfrak { s }$$

and δ jk denotes the Kronecker delta function. Next, the expression for collision in Eq. (5) can be rewritten as,

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \delta$$

Using Eqs. (3) and (10), the quadratic form of the collision step can be written as

$$& \text{step can be wrtten as} \\ \text{the} \\ \text{and} \\ \text{ables}$$

where β ijk = γ ijk · ( 1 -∆ t τ ) + χ jk · ∆ t τ α ijk and T denotes the transpose. The indicator functions γ and χ for i th lattice direction is given by,

$$\gamma _ { _ { i j k } } = \delta _ { i j } + \delta _ { i k } - \delta _ { i j } \cdot \delta _ { i k } \quad \quad ( 1 3 )$$

$$\gamma _ { i j k } & = \delta _ { i j } + \delta _ { i k } - \delta _ { i j } \cdot \delta _ { i k } \quad \quad ( 1 3 ) \\ \chi _ { j k } & = \begin{cases} 1 & \text{if $j\leq k$} \\ 0 & \text{otherwise} \end{cases}$$

Matrix formulation for weakly compressible LBM. In case of weakly compressible flows, ρ -1 ≈ 2 -ρ , we then express

$$\int _ { \ell = 1 } ( 2 - \rho ) = ( - 1, - 1, \cdots, - 1, 2 ) \cdot ( f _ { 1 }, f _ { 2 }, \cdots, f _ { n _ { \ell } }, 1 ) ^ { \mathrm T } \ ( 1 5 )$$

where the value 1 appended at the end of PDFs can be treated as a auxiliary constant. Therefore, the collision step can be written as

$$\intertext { e t o } \, \begin{pmatrix} \quad \,. \\ \quad \, \left ( \begin{smallmatrix} f _ { i } ^ { * } \\ 1 \end{smallmatrix} \right ) = \left ( \begin{smallmatrix} f \\ 1 \end{smallmatrix} \right ) ^ { \mathrm T } \left ( \begin{smallmatrix} \beta _ { i } & 0 \\ 0 & 1 \end{smallmatrix} \right ) \left ( \begin{smallmatrix} \text{Diag} ( f ) & 0 \\ 0 & 1 \end{smallmatrix} \right ) \, \mathrm W \left ( \begin{smallmatrix} f \\ 1 \end{smallmatrix} \right ) \, \end{pmatrix} \, \begin{smallmatrix} \quad \, ( 1 6 ) \\ \quad \, \left ( \begin{smallmatrix} \quad \,. \\ 1 \end{smallmatrix} \right ) \, \end{pmatrix}$$

where f = ( f , f 1 2 , · · · , f n e ), and β i represent the upper triangular matrix in Eq. (12). The term Diag ( f ) represents the diagonal matrix with entries from f . The matrix

W is given by

$$\dot { \ } g \text{ven by } & & & & & & & & & &$$

## Quantum Algorithm. -

Encoding. Let n g denote the total number of grid points and the number of PDFs will be n f = n n e g . According to the lattice directions, we arrange and order the PDFs in a vector form, df = ( f 1 , f 2 , . . . , f n e , 1), where the suffix refers to the direction, and each f i is of length n g . Since PDFs are encoded as a quantum state-vector, we define φ = ( df , df , df , df ). A computational register containing n q = log ( 2 n f ) qubits along with two ancilla qubits will be required to encode the PDFs. The initial state to encode will be

$$\left | \phi ^ { 0 } \right \rangle = \left | 0 \right \$$

Collision. Before we express matrix version of Eq. (16), we define the following matrices:

$$\widetilde { W } & = \begin{pmatrix} - \mathbf I _$$

$$\widetilde { \mathbf F } _ { i } = \begin{pmatrix}$$

$$\begin{array} {$$

where I g rc denotes the identity matrix of size n g × n g , 2 g i = (2 , · · · , 2) and 0 g i = (0 , · · · , 0) are the vectors of size n g . The term ˜ β ric = χ ic β ric I g rc , where β ric is defined in Eq. (12), with r and c represent row and column indices, respectively. Based on the above set of matrices, we

$$\begin{array} { c c c c c c c c } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & && & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & - & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & _ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ; & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &   & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & . & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  ; & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & 
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &.. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &.
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &$$

Therefore, the collision operator is given by

| φ ∗ 〉 = ̂ F n e ̂ B n e · · · ̂ F 2 ̂ B F 2 ̂ 1 ̂ B D W 1 ̂ ̂ | φ 0 〉 (25) In Eq. (25), when the SV is multiplied by ̂ W yields the value 2 -ρ , which then multiplied by ̂ D results in the first order approximation of f i ρ , i.e., the column vector defined on the right side of Eq. (12). Upon successive multiplication by ̂ F B i ̂ i , yields the PDFs in the i th direction. Thus, for n e velocity directions, the collision operation is decomposed into 2 n e +2 operations. It is important to note that the state vector has been initialized with four sets of df . Among these, one will ultimately lead to the final outcome, while another is designated for storing temporary values. The remaining two sets remain constant, providing the initial state vector values during intermediate computations. Consequently, the matrices defined in Eq. (24) are not unique. In particular, the following is a valid choice:

$$( 2 1 ) & \quad \text{not unique. In particular, the following is a valid choice:} \\ \widehat { B } _ { i } = \begin{pmatrix} I & 0 & 0 & 0 \\ 0 & I & 0 \\ \widetilde { B } _ { i } & 0 & 0 \\ I \end{pmatrix} ; \ \widehat { \mathbf F } _ { i } = \begin{pmatrix} I & 0 & \widetilde { \mathbf F } _ { i } & 0 \\ 0 & 0 & I \\ I & 0 \\ I \end{pmatrix} \quad ( 2 6 ) \\ \quad \text{Streaming and Boundary Conditions.} \quad \text{The streaming} \\ \quad \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot \ \cdot$$

Streaming and Boundary Conditions. The streaming matrix S consists of n e × n e blocks with the size of each block n g × n g . The matrix S is square binary matrix where each row/column contains one entry equal to 1. Since every row index i can be decomposed as x + yn x + i e n n x y , and the corresponding column index j such that s ij = 1, can be computed by replacing the co-ordinate x by x -e i ∆ . In the present study, we utilized the wall boundary t condition, accordingly the column index will be computed as x ∆ + y ∆ n x + i -e n n x y , where ( x ∆ , y ∆ ) = x -e i and i -e is the reflection of the direction i e . Based on the matrix S , we define the streaming operator as

$$\begin{array} { c c c } S, \, \text{we define the streaming operator as} \\ n _ { g }, \\ \intertext { of } \begin{pmatrix} \widetilde { S } & \\ & I \\ & & I \\ \end{pmatrix} & ; \, \text{where } \widetilde { S } = \begin{pmatrix} S & \\ & 1 \end{pmatrix} \quad ( 2 7 ) \\ \intertext { \. }$$

From Eq. (27), it is clear that streaming is performed only on the computational qubits without any involvement from ancilla qubits. Thus, the PDF for the next time step can be obtained from the quantum state,

## Results and Discussion. -

$$| \phi ^ { n } \rangle = \widehat { S } \, | \phi ^ { * } \rangle \quad \quad ( 2 8 ) \quad \text{$\upsilon_{c}$} \\ \text{Discussion}. \, - \, \quad \ \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{\quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{\quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsi_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$} \quad \text{$\upsilon_{c}$,}$$

Verification - Case 1: 1D discontinuity flow. The classical shock tube discontinuity problem in one dimension is used to verify the present the algorithm. Because of the weak compressibility assumption, we refer it as a discontinuity rather than a shock. The one dimensional domain of size n g = 500 grid points and the D Q 1 3 lattice has been chosen. At the initial time, the velocity is set to zero everywhere and density is defined as,

$$\rho ( x, 0 ) = \begin{cases} 1. 0 + \Delta \rho & \text{if $x\leq \frac{n_{q}}{2}$} \\ 1. 0 & \text{otherwise} \end{cases} \quad ( 2 9 ) \quad a \text{fun} \\ \text{son o}$$

where ∆ ρ = 5 × 10 -5 is chosen to satisfy the weak incompressibility assumption made earlier. Bounce-back boundary condition is applied at the both ends of the domain. After simulating 200 time steps, we obtained a good comparison with the exact Riemann solution for normalized pressure ( p ∗ = ( p -p s ) /p s ; p s = ∆ ρc 2 s ) and velocity ( u ) [20] and is given in Fig. 1. The relative difference between exact and QLB results is to found to be 5%.

Fig. 1: Comparison of normalized pressure and velocity of 1D discontinuity flow obtained from QLB with the exact Riemann solution at 200 th time step.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00387v1-figure-aa345a1f-0d1d-4c9c-b55c-b9b8099a7175.jpg)

Case 2: 2D Kolmogorov flow. In order to compare the accuracy of present algorithm with that of Carleman linearized - LBM presented in [19], we performed the twodimensional simulation of Kolmogorov-like flow on a 32 × 32 grid. The D2Q9 lattice is used. Initially, the PDFs are defined as,

$$f _ { i } ( x, y ) = w _ { i } \left [ 1 + A _ { x } \cos \left ( \frac { 2 \pi k _ { x } } { N _ { y } } y \right ) \text{e} _ { i } \cdot \text{e} _ { x } + \text{the $rc$} \text{The r} \\ A _ { y } \cos \left ( \frac { 2 \pi k _ { y } } { N _ { x } } x \right ) \text{e} _ { i } \cdot \text{e} _ { y } \right ] \text{ } ( 3 0 ) \text{ the stat} \\ LB \text{ ope}$$

where e x = (1 0) and , e y = (0 1). , The parameters in Eq. (30) are taken as: A x = 0 3, . A y = 0 2, . k x = 1, and k y = 4. We choose different viscosity values ν ranging from 1 / 6 to 0 0088, and for each . ν we ran the simulation up to 100 time steps. After simulating 100 time steps, we compute the mean value of Root Mean Squared Error (RMSE) between the distribution functions obtained from the present algorithm ( f p i ) and the exact LBM simulation ( f e i ),

$$\stackrel { \nabla _ { i } \, \jmath } { s } - \bigvee _ { i } \nabla _ { i } \, \jmath, \\ \text{is} \, \ e \, \begin{array} { c } \nabla _ { i } \\ \nabla _ { i } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end {array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \, \left \langle \begin{array} { c } \nabla _ { j } \end{array} \right \rangle$$

Since CL transforms a finite set to an infinite dimensional set of equations, [19] analysed two different approaches, truncation (neglecting higher order terms) and closure (approximate the product of functions into products of a function and a constant). Fig. 2 shows the comparison of RSME with the results of Carleman second order truncation and closure approaches presented in [19]. Since the present approach does not require any order of truncation except the assumption on weak compressiblity, the &lt; RSME &gt; value is below 10 -5 compared to 10 -2 from Carleman closure approach.

Fig. 2: Comparison of root mean square error obtained for various viscosity values are compared with thre results of carleman second order truncation and closure approaches presented in [19].

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00387v1-figure-c9497803-5cd1-479e-8f69-df7004b24cd8.jpg)

Computational Complexity. In the present approach, the number of qubits required will be 2+log 2 ( n f ), whereas the Carleman truncated system of k th order would require 1+log ( 2 n f + n 2 f + . . . + n k f ) qubits. Fig. 3a shows the number of qubits required for the grid size ranging from 10 1 to 10 20 using D2Q9 lattice. For larger grid sizes, the CLsecond (CL2) and third (CL3) order truncation approach requires almost twice and thrice the qubit resource than the present approach.

The number of two-qubit gates is one of the ways to estimate the complexity of the present quantum algorithm. Since PDFs are encoded as amplitudes of the qubit states, the state preparation step is required before performing LB operations. According to [21], O (2 n q ) CNOT gates

are required for state preparation. Based on Eq. (25) and (28), the number of LB operators per time step is 2 n e +3, and each operator would require O (2 n q -1 (2 n q -1)) CNOT gates. In the second-order Carleman truncation, the qubit requirement will be n c q = 1 + log ( 2 n f + n 2 f ). Correspondingly, the O (2 n c q -1 (2 n c q -1)) CNOT gates will be required. Fig. 3b shows the number of CNOT gates required for different grid sizes. For smaller grids of size 10 3 , the gate counts of CL2 are a few orders higher than the present approach; however, for the larger grid, it is exponentially larger.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00387v1-figure-73d70119-d063-4a35-84c1-9627a8b68272.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00387v1-figure-b3cf773e-3d77-4b44-b644-71ef5bba8182.jpg)

nCNOT

Fig. 3: Comparison of (a) number of qubits and (b) number of CNOT gates required for number of grid points between the present algorithm vs. the Carleman second and third order truncation.

Summary. We have developed a quantum algorithm to deal with the non-linearity in the Lattice Boltzmann collision operator. By decomposing the collision operator, we significantly reduced the circuit width by half and circuit depth in exponential order when compared to the Carleman linearization technique. While there are still challenges to be addressed, such as the encoding and readout process for each time step and the construction of quantum gates, this work aims to minimize the quantum resources required by the CL technique and prevent the formation of an infinite system. Due to the weakly compressible assumption and the first-order approximation for 1 /ρ , this work is limited to moderate Reynolds numbers of O(100).

Acknowledgments. This research was kindly supported by the Quantum Computing Consortium that has been funded by the MAGNET program of the Israel Innovation Authority (IIA).

## REFERENCES

- [1] Harrow A. W., Hassidim A. and Lloyd S. , Physical Review Letters , 103 (2009) .
- [2] Cao Y., Anmer D., Frankel S. H. and Sabre K. , Molecular Physics , 110 (2012) 1675.
- [3] Bharadwaj S. S. and Sreenivasan K. R. , Proceedings of the National Academy of Sciences , 120 (2023) e2311014120.
- [4] Ingelmann J., Bharadwaj S. S., Pfeffer P., Sreenivasan K. R. and Schumacher J. , Computers &amp; Fluids , (2024) 106369.
- [5] Bharadwaj S. S. and Sreenivasan K. R. , arXiv preprint arXiv:2405.09767 , (2024) .
- [6] Kuznik F., Obrecht C., Rusaouen G. and Roux J.-J. , Computers &amp; Mathematics with Applications , 59 (2010) 2380.
- [7] Steijl R. , Applied Sciences 2023, Vol. 13, Page 529 , 13 (2022) 529.
- https://www.mdpi.com/2076-3417/13/1/529/htmhttps://www.mdpi
- [8] Steijl R. , , (2024) .
- [9] Mezzacapo A., Sanz M., Lamata L., Egusquiza I. L., Succi S. and Solano E. , Scientific Reports 2015 5:1 , 5 (2015) 1.

## https://www.nature.com/articles/srep13153

- [10] Budinski L. , Quantum Information Processing , 20 (2021) 57.
- [11] Budinski L. , Int. J. Quantum Information , 20 (2022) 2150039.
- [12] Schalkers M. A. and M¨ oller M. , J. Comp. Phys. , 502 (2024) 112816.
- [13] Dinesh Kumar E. and Frankel S. H. , arXiv:2405.08669 , (2024) .
- [14] Succi S., Itani W., Sreenivasan K. and Steijl R. , Europhysics Letters , 144 (2023) 10001.

https://dx.doi.org/10.1209/0295-5075/acfdc7

- [15] Carleman T. , , (1932) .
- [16] Itani W. and Succi S. , Fluids , 7 (2022) 24.
- [17] Itani W., Sreenivasan K. R. and Succi S. , Phys. Fluids , 36 (2024) 017112.

## https://doi.org/10.1063/5.0176569

- [18] Sanavio C., Scatamacchia R., de Falco C. and Succi S. , Physics of Fluids , 36 (2024) .
- [19] Sanavio C. and Succi S. , AVS Quantum Science , 6 (2024) 023802.

https://doi.org/10.1116/5.0195549

- [20] Toro E. F. , Riemann solvers and numerical methods for fluid dynamics: a practical introduction (Springer Science &amp; Business Media) 2013.
- [21] Shende V. V., Bullock S. S. and Markov I. L. , Synthesis of quantum logic circuits in proc. of Proceedings of the 2005 Asia and South Pacific Design Automation Conference 2005 pp. 272-275.
Dinesh Kumar E Steven H. Frankel
Decomposition of Nonlinear Collision Operator in Quantum Lattice Boltzmann Algorithm
We propose a quantum algorithm to tackle the quadratic nonlinearity in the
Lattice Boltzmann (LB) collision operator. The key idea is to build the quantum
gates based on the particle distribution functions (PDF) within the coherence
time for qubits. Thus, both the operator and a state vector are linear
functions of PDFs, and upon quantum state evolution, the resulting PDFs will
have quadraticity. To this end, we decompose the collision operator for a
$DmQn$ lattice model into a product of $2(n+1)$ operators, where $n$ is the
number of lattice velocity directions. After decomposition, the $(n+1)$
operators with constant entries remain unchanged throughout the simulation,
whereas the remaining $(n+1)$ will be built based on the statevector of the
previous time step. Also, we show that such a decomposition is not unique.
Compared to the second-order Carleman-linearized LB, the present approach
reduces the circuit width by half and circuit depth by exponential order. The
proposed algorithm has been verified through the one-dimensional flow
discontinuity and two-dimensional Kolmogrov-like flow test cases.