Analysis of Carleman Linearization of Lattice Boltzmann

Abstract

We explore the Carleman linearization of the collision term of the lattice Boltzmann formulation, as a first step towards formulating a quantum lattice Boltzmann algorithm. Specifically, we deal with the case of a single, incompressible fluid with the Bhatnagar Gross and Krook equilibrium function. Under this assumption, the error in the velocities is proportional to the square of the Mach number. Then, we showcase the Carleman linearization technique for the system under study. We compute an upper bound to the number of variables as a function of the order of the Carleman linearization. We study both collision and streaming steps of the lattice Boltzmann formulation under Carleman linearization. We analytically show why linearizing the collision step sacrifices the exactness of streaming in lattice Boltzmann, while also contributing to the blow up in the number of Carleman variables in the classical algorithm. The error arising from Carleman linearization has been shown analytically and numerically to improve exponentially with the Carleman linearization order. This bodes well for the development of a corresponding quantum computing algorithm based on the lattice Boltzmann equation.

1. Introduction

Computational Fluid Dynamics (CFD) has accompanied computers as an application since their infancy, starting with von Neumann’s program to simulate the weather on the ENIAC machine around the 1950s. Even earlier, in 1922, Richardson described “human” computers computing the weather by hand, estimating that 64,000 of them, each calculating at a speed of 0.01 Flops/s, would be sufficient to predict the weather in real time [1]. Leaving aside human calculators, electronic ones have made a spectacular trek so far, from the few hundred Flops of ENIAC to the exaflop peak performance of the Sunway Oceanlite supercomputer [2]. This has spanned sixteen orders of magnitude in 70 years, close to a sustained Moore’s law rate, doubling every 1.5 years! Amazingly, CFD has been consistently on the forefront of the journey, and it continues to be to the present day. However, when it comes to quantum computing, CFD is yet to capture the limelight it deserves. In this paper, we present a brief survey of current ongoing research work in this direction and a preparatory technique, known as Carleman linearization, aimed at the development of a quantum computing algorithm for the lattice Boltzmann method for fluid flows.

1.1. Early Attempts for Quantum Simulation of Fluids

The earliest attempts at quantum simulation of fluids have been based on the lattice gas or lattice Boltzmann algorithms. The first quantum lattice Boltzmann scheme dates back to the 1990s [3,4]. Around the turn of the millenium, Yepez demonstrated fluid dynamic simulations on a special-purpose quantum computer based on nuclear magnetic resonance (NMR) [5], using the quantum lattice gas algorithm [6,7]. Leading the trail, Yepez has also investigated Burger’s equation [8,9], and entropic lattice Boltzmann models [10]. The latter has found its use in simulation of quantum fluid dynamics and other quantum systems [6]. We have recently seen a divergence towards citing Navier-Stokes as a future direction for papers discussing solving nonlinear differential equations on quantum computers, away from the early physically-motivated algorithms for fluid simulation, as the quantum lattice gas one mentioned above [5,7,8,9]. However, the attempts to revive the physically-motivated algorithms beyond quantum systems [11], by [12,13,14,15] and others, are promising. The work of [13] stands out as it presents itself as a method not of quantum computation per se, but of quantum simulation. The latter leverages the correspondence between the Dirac and lattice Boltzmann equations. We, thus, find a compelling reason [16] to explore lattice Boltzmann as the basis for quantum simulation of fluids, starting with its linearization, explored classically in this paper.

1.2. Carleman Linearization

Carleman linearization appears to cast linearization of a function through Taylor series expansion into matrix form, suitable for use in defining state estimator of a non-linear system of known dynamics as part of the Koopman operator approach, in what is known as Carleman-Koopman operator. The basic idea of Carleman linearization is to introduce powers of the original variable as a variables in the system. The recurrence through which the new equations in the system is defined leads to an infinite dimensional system. The latter is prone to admitting solutions different than those of the original system [17]. The linearization is achieved at the cost of infinite-dimensions. Upon truncation of the resulting infinite-dimensional upper-triangular matrix [18], the accuracy of the approximation also suffers, deteriorating with time, making the truncated system most suitable for asymptotically stable stationary solutions [19]. The error bound for a Carleman-linearized polynomial ODE reduced to quadratic form has been shown, through a power series approach [18], to depend on the initial condition as well as exponentially on time. They recommend discretizing the solution in time, and evaluating other basis functions. Their work has readily been extended by [20] for a quantum algorithm. Apart from discussing the complexity and error bounds of a quantum Carleman algorithm, [20] presented the results for a classical Carleman linearization of Burger’s equation.

2. Lattice Boltzmann

Between the Navier-Stokes equations which model the flow at a continuum level, and molecular dynamics which treat the microscale, LB stands out, representing the fluid as an ensemble of particles at the mesoscale. It stems from the minimal discretization of the Boltzmann kinetic equation. The lattice Boltzmann formulation is readily extensible for a range of physics, from the quantum to the relativistic, giving the method a versatility about which books could be written [21].

$$\frac{df}{dt}=\frac{\partial f}{\partial t}+\overrightarrow{v}·\overrightarrow{\nabla }f=\mathsf{\Omega }$$

(1)

The Boltzmann Equation (1) describes the probability density f of the fluid in the position-momentum space, driven by advection due to continuum particle velocity $\overrightarrow{v}$ and collision $\mathsf{\Omega }$ across space spanned by $\overrightarrow{x}$, and time t. To arrive at the lattice Boltzmann formulation, the probability density f from the Boltzmann equation Equation (1) is discretized into Q density distributions, each describing the fraction of fictitious particles: moving in a given D-dimensional lattice, with $\overrightarrow{v}$ restricted to speeds ${\overrightarrow{c}}_i$, ${\overrightarrow{c}}_i=c_i{\overrightarrow{e}}_i$, defined in the directions of the lattice vectors ${\overrightarrow{e}}_i$s.

The discretized lattice Botlzmann equation takes the form

$$\frac{1}{\Delta t}(f_i(\overrightarrow{x}+{\overrightarrow{c}}_i\Delta t,t+\Delta t)-f_i(\overrightarrow{x},t))=-\frac{1}{\tau }(f_i(\overrightarrow{x},t)-f_i^{eq}(\overrightarrow{x},t))$$

(2)

for the simple case of a single-phase fluid, which components are permeable to each other such that they can be considered well-mixed and homogeneous [22]. The discrete probability density is the fundamental variable of the lattice Boltzmann approach. It describes the probability of finding a fictitious fluid particle at a given point defined by the position vector $\overrightarrow{x}$, at a given instant of time t, with a particular speed [23]. Implicit, in the discretization of the Boltzmann equation, is the discretization of the phase-space, involving discrete position and velocity values. The lattice Boltzmann proceeds by updating the discrete probability densities at each cell in two steps: collision and advection/streaming shown in Figure 1. This is the collision step:

$$f_i(\overrightarrow{x},t+\Delta t)=f_i(\overrightarrow{x},t)-\frac{\Delta t}{\tau }(f_i(\overrightarrow{x},t)-f_i^{eq}(\overrightarrow{x},t))$$

(3)

The streaming step is:

$$f_i(\overrightarrow{x},t+\Delta t)\stackrel{}{\to }{f}_i(\overrightarrow{x}+{\overrightarrow{c}}_i\Delta t,t+\Delta t)$$

(4)

Collision is practically a relaxation towards equilibrium, a nonlinear operation dependent on terms local to each cell. On the other hand, streaming involves the transfer of the discrete densities to nearby cells, a nonlocal linear operation. The equilibrium distribution $f_i^{eq}$ is written as a function of $c_i=c$ the lattice speeds, the fluid density $\rho$, the lattice vectors ${\overrightarrow{e}}_i$, and the flow velocity $\overrightarrow{u}$.

$$\rho ={\mathsf{\Sigma }}_{i=1}^Q{f}_i$$

(5)

flow velocity $\overrightarrow{u}$:

$$\overrightarrow{u}=\frac{c}{\rho }{\mathsf{\Sigma }}_{i=1}^Q{f}_i{\overrightarrow{e}}_i$$

(6)

weights w which sum to unity. A common model for the equilibrium function is:

$$f_i^{eq}(\overrightarrow{x},t)=w_i\rho (1+(\frac{3}{c}{\overrightarrow{e}}_i·\overrightarrow{u}+\frac{9}{2c^2}{({\overrightarrow{e}}_i·\overrightarrow{u})}^2-\frac{3}{2c^2}\overrightarrow{u}·\overrightarrow{u}))$$

(7)

Replacing the expression for the equilibrium expression for the incompressible case, we have:

$$\begin{array}{cc}\hfill (f_i(\overrightarrow{x}+{\overrightarrow{c}}_i\Delta t,t+\Delta t)=& (1-\frac{\Delta t}{\tau })f_i(\overrightarrow{x},t))\hfill \\ & +\frac{\Delta t{w}_i}{\tau }(1+3{\overrightarrow{e}}_i·f_j{\overrightarrow{e}}_j+\frac{9}{2}{({\overrightarrow{e}}_i·f_j{\overrightarrow{e}}_j)}^2-\frac{3}{2}{f}_j{f}_k{\overrightarrow{e}}_j·{\overrightarrow{e}}_k)\hfill \end{array}$$

(8)

We note that the lattice vectors ${\overrightarrow{e}}_i$, do not correspond to unit vectors, and are generally not orthonormal. In a single dimension, the left and right vectors have nonzero inner product of $-1$. In higher dimensions, we can identify such a pair for each direction, in addition to nonzero inner products with and between diagonal lattice vectors shown in Figure 2. For example, in D1Q3, we have:

$$e_i=\left\{\begin{array}{cc}-1\hfill & i=1\hfill \\ 0\hfill & i=2\hfill \\ 1\hfill & i=3\hfill \end{array}\right.$$

(9)

thereby:

$$\overrightarrow{\mathsf{\Omega }}=-\frac{dt}{\tau }\left(\begin{array}{c}\frac{1}{2}(f_1+f_3-{(f_1-f_3)}^2-\frac{1}{3})\\ (f_2+{(f_1-f_3)}^2-\frac{2}{3})\\ \frac{1}{2}(f_1+f_3-{(f_1-f_3)}^2-\frac{1}{3})\end{array}\right)$$

(10)

Continuing with the consideration of the D1Q3 case, by introducing the vectors of first-order and second-order products of the discrete densities appearing in the lattice Boltzmann formulation, denoted ${\overrightarrow{f}}^{\left(1\right)}$ and ${\overrightarrow{f}}^{\left(2\right)}$ respectively:

$$\begin{array}{c}{\overrightarrow{f}}^{\left(1\right)}={\left(\begin{array}{ccc}{f}_1& f_2& f_3\end{array}\right)}^T\\ {\overrightarrow{f}}^{\left(2\right)}={\left(\begin{array}{ccc}{f}_1^2& f_3^2& f_3{f}_1\end{array}\right)}^T\end{array}$$

(11)

we may write $\overrightarrow{\mathsf{\Omega }}$ as a second-order “polynomial” with matrix-valued coefficients, similar to the form introduced by [18].

$$\overrightarrow{\mathsf{\Omega }}=F_1\overrightarrow{f}\left(1\right)+F_2{\overrightarrow{f}}^{\left(2\right)}+{\overrightarrow{F}}_0$$

(12)

where $F_1$ and $F_2$, for the D1Q3 are:

$$F_1=-\frac{\Delta t}{\tau }\left(\begin{array}{ccc}\frac{1}{2}& 0& \frac{1}{2}\\ 0& 1& 0\\ \frac{1}{2}& 0& \frac{1}{2}\end{array}\right)$$

(13)

and:

$$F_2=-\frac{\Delta t}{\tau }\left(\begin{array}{ccc}-\frac{1}{2}& -\frac{1}{2}& 1\\ 1& 1& -2\\ -\frac{1}{2}& -\frac{1}{2}& 1\end{array}\right)$$

(14)

The constant vector ${\overrightarrow{F}}_0$ appearing could be absorbed with a fixed variable, say $f_4=1$ to account for constants. It could easily be defined with a zero derivative adding an empty row to the bottom of $F_1$ and $F_2$. The constants appearing in any polynomial entry of the vector could be, thus, defined as 1 or $1^2$ multiplied by some coefficient. Therefore, they should not affect the norm of $F_1$ or $F_2$ considered below.

$${\overrightarrow{F}}_0=\frac{\Delta t}{\tau }\left(\begin{array}{c}\frac{1}{6}\\ \frac{2}{3}\\ \frac{1}{6}\end{array}\right)$$

(15)

The fact that the rows and columns of $F_1$ sum up to unity, whereas those $F_2$ sum to zero is at heart of the favourable error bounds obtained in Section 3.4, and that in presence of uniform initial conditions or absence of streaming, the method becomes exact. It is well-known that Equation (8) is inconsistent with the incompressible assumption, $\rho =1$, and density fluctuations around unity are expected to arise during evolution [24]. Neglecting these fluctuations amounts to an error in the velocity components $\left|\right|\overrightarrow{u}\left|\right|$ proportional to the square of the Mach number $O\left(M{a}^2\right)$. This is due to the fact that this standard lattice Boltzmann formulation recovers the compressible Navier-Stokes upon expansion. More sophisticated formulations have been developed for incompressible and nearly incompressible flows [25]. The particles with random motion, are restricted to the lattice nodes with microscopic velocities $c_i$ defined over lattice directions, allowing us to model the collision of particles and their streaming in seperate, uncoupled steps. The latter forces the nonlinearity of fluid flow, captured in the collision term, to be local, whereas the non-local streaming terms remain linear. Moreover, streaming is exact.

Nonlinearity Ratio

We define the nonlinearity ratio R as a measure of how much Equation (8) deviates from a complete linear behavior. It is defined as

$$R=\left|\right|f(t=0)\left|\right|\frac{\left|\right|F_2\left|\right|}{\left|\mathrm{Real}\right({\lambda }_{MAX}\left(F_1\right)\left)\right|},$$

(16)

where $F_1$ is the matrix of linear coefficients of Equation (8) and $F_2$ is the matrix of the second order terms. Here we have considered the supremum norm of a matrix $\left|\right|A\left|\right|$ or a vector $\left|\right|\overrightarrow{x}\left|\right|$ as

$$\left|\right|A\left|\right|=\underset{ij}{max}|A_{ij}|,||\overrightarrow{x}\left|\right|=\underset{i}{max}\left|x_i\right|$$

(17)

and ${\lambda }_{MAX}\left(F_1\right)$ as the largest eigenvalue in modulo of $F_1$. Note that for a linear system $\left|\right|F_2\left|\right|=0$ such that one has no nonlinearity. This makes R “qualitatively” similar to the Reynolds number which is a ratio of nonlinear convective forces to linear viscous forces. When we rescale a matrix by a constant, its eigenvalues get rescaled as well as the norm. As such, we see that the $\frac{\Delta t}{\tau }$ factor cancels out from the nonlinearity ratio. This is not surprising given that R is practically a measure of the relative weight of nonlinear terms to linear terms. However, this independence of R on a critical parameter, $\tau$, should be taken with cautiously, as the it role is more nuanced in the lattice Boltzmann formulation.

We may proceed by calculating it for the case of $\frac{\Delta t}{\tau }=1$, and it remains valid otherwise.

In the case of D1Q3, the two different norms taken for $F_1$ and $F_2$ both give unity, such that:

$$R_{D1Q3}=\left|\right|f(t=0)\left|\right|\le 1$$

(18)

where it is less than or equal to 1 by definition of a probability distribution, and where the equality holds in the unlikely case that the distribution is initialized such that a single discrete density is unity. The fact that lattice Boltzmann’s variables are formulated as probability densities is indeed advantageous. While a similar feat could be achieved by rescaling variables for a well-posed problem, as suggested by [20], their relative values and numerical accuracy used limit the effectiveness of such an approach. Moreover, $Re\left({\lambda }_{MAX}\right)=-1<0$, this means that both conditions set by [20] for arbitrary time-convergence are met. Another advantage of the lattice Boltzmann method is that it is by definition discretized, such that the error arising from discretizing the equation, the forward Euler for time derivative for example as discussed in [20], is irrelevant.

3. Carleman Linearization for Lattice Boltzmann

The basic tenant of Carleman linearization is a change of variables done such that the variables of the original nonlinear system, $\overrightarrow{f}\left(t\right)$, are replaced by a larger set of variables $\overrightarrow{V}\left(t\right)$. The original variables form a subset of the larger set of Carleman variables. The additional variables are monomials up to order $O_c$ of the original $f_i$.

Denoting the vector of Carlemann variables of kth order as $V^{\left(k\right)}$, and $P^{\left(k\right)}$ a vector of some polynomial functions of kth order in Carlemann variables, using the numbering scheme shown in Figure 3, we may write:

$$\begin{array}{cc}\hfill {\overrightarrow{V}}^{\left(1\right)}={\left(\begin{array}{ccc}{f}_1& f_2& f_3\end{array}\right)}^T,& {\partial }_t{\overrightarrow{V}}^{\left(1\right)}={\overrightarrow{P}}^{\left(2\right)}({\overrightarrow{V}}^{\left(1\right)},{\overrightarrow{V}}^{\left(2\right)})\hfill \\ \hfill {\overrightarrow{V}}^{\left(2\right)}={\left(\begin{array}{cccccc}{f}_1^2& f_1{f}_2& f_2^2& f_2{f}_3& f_3^2& f_3{f}_1\end{array}\right)}^T,& {\partial }_t{\overrightarrow{V}}^{\left(2\right)}={\overrightarrow{P}}^{\left(3\right)}({\overrightarrow{V}}^{\left(1\right)},{\overrightarrow{V}}^{\left(2\right)},{\overrightarrow{V}}^{\left(3\right)})\hfill \\ \hfill \cdots \end{array}$$

(19)

For lattice Boltzmann with the BGK equilibrium function, we see that the time evolution of each set of monomial of discrete density variables of given order could be written as a function of variables of same order, all preceding orders, and an order above. In Equation (19), we see that those of first order monomials are written in terms of first and second order monomials, and those of second order monomials are written in terms of first, second and third order monomials.

$$\begin{array}{cc}\hfill {\overrightarrow{V}}^{\left(1\right)}={\left(\begin{array}{ccc}{f}_1& f_2& f_3\end{array}\right)}^T,& {\partial }_t{\overrightarrow{V}}^{\left(1\right)}={\overrightarrow{P}}^{\left(2\right)}({\overrightarrow{V}}^{\left(1\right)},{\overrightarrow{V}}^{\left(2\right)})\hfill \\ \hfill {\overrightarrow{V}}^{\left(2\right)}={\left(\begin{array}{cccccc}{f}_1^2& f_1{f}_2& f_2^2& f_2{f}_3& f_3^2& f_3{f}_1\end{array}\right)}^T,& {\partial }_t{\overrightarrow{V}}^{\left(2\right)}={\overrightarrow{P}}^{\left(3\right)}({\overrightarrow{V}}^{\left(1\right)},{\overrightarrow{V}}^{\left(2\right)},\sout{{\overrightarrow{V}}^{\left(3\right)}})\hfill \end{array}$$

(20)

Carleman embedding involves, then, choosing an order of linearization and neglecting higher-order terms from the driving functions of ODEs. This is what we refer to as truncation, and it is demonstrated in Equation (20) to second order, where third-order terms otherwise required for the time evolution of second order terms are dropped.

$$\begin{array}{c}\overrightarrow{V}=\left(\begin{array}{c}{\overrightarrow{V}}^{\left(1\right)}\\ {\overrightarrow{V}}^{\left(2\right)}\end{array}\right)=\left(\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_1^2\\ f_1{f}_2\\ f_2^2\\ f_2{f}_3\\ f_3^2\\ f_3{f}_1\end{array}\right)=\left(\begin{array}{c}{V}_1\\ V_2\\ V_3\\ V_4\\ V_5\\ V_6\\ V_7\\ V_8\\ V_9\end{array}\right)\end{array}$$

(21)

The Carlemann variables $\overrightarrow{V}$ are then introduced to replace the remaining monomials, as shown in Equation (21).

An example of additional variables of second order is shown in

Table 1. Example of D1Q3 expanded to a second order truncation in Carleman linearization, with $N=9$ Carleman variables. Note that the dummy variable $V_1=1$ is defined to simplify the form of Equation (23).

Carleman Variables	Lattice Variables
$V_1$	1
$V_2$	$f_1^2$
$V_3$	$f_3^2$
$V_4$	$f_1{f}_2$
$V_5$	$f_1{f}_3$
$V_6$	$f_2{f}_3$
$V_7$	$f_1$
$V_8$	$f_2$
$V_9$	$f_3$

for a D1Q3 lattice, taking into account the model in Equation (8).

The dynamics of the extended system are then derived from the original system of a single phase, homogeneous, fluid, with terms beyond a chosen order $O_c$ dropped.

$$\frac{\partial f_i\left(t\right)}{\partial t}={\mathsf{\Omega }}_i\left(\overrightarrow{f}\left(t\right)\right)$$

(22)

which is linearized into the total derivative of the Carleman variables vector $\overrightarrow{V}\left(t\right)$ equating a constant coefficient matrix, the Carleman linearization matrix C, multiplying the Carleman variables vector:

$$\frac{\partial \overrightarrow{V}\left(t\right)}{\partial t}=C\overrightarrow{V}\left(t\right)$$

(23)

where the Carleman linearization matrix C is obtained by deriving the following system of equations:

$$\frac{\partial V_n}{\partial t}={\mathsf{\Sigma }}_{i=1}^Q\frac{\partial V_n\left(t\right)}{\partial f_i}{\mathsf{\Omega }}_i\left(f\left(t\right)\right)$$

(24)

and identifying the variables $V_i$ on the right hand side of Equation (24).

3.1. Number of Variables

An upper bound for the number of Carleman variables for a desired order $O_c$, when considering N original lattice variables, is

$$N=\frac{(O_c+(Q+1)-1)!}{\left(O_c\right)!((Q+1)-1)!}=\left(\begin{array}{c}{O}_c+Q\\ Q+1\end{array}\right)$$

(25)

This is an upper bound to the number of Carleman variables used in the linearized system because in practice some terms of order $O_c$ do not appear in the linearized equations. For example, from Table 1 $f_2^2$, the second order term for the rest particles density, does not appear in the D1Q3 formulation. In general, the exact number of variables is given after specifying a lattice structure and an equilibrium function.

3.2. Carleman Linearization of Collision Step

The collision operator Equation (3) can be expressed as a function of the Carleman variables, $V_i$s. To do that, we go back to Equation (1), and consider collision sans streaming. The total derivative of the discrete densities is described by the nonlinear collision operator on the right-hand side.

This assumes that the obtained solution for the discrete density distribution at each site is streamed exactly, and the nonlinear terms recalculated to evolve the system. While this allows us to isolate the error from the linearization of the collision step, it is undesirable in practice. On another note, the linearization of the collision step allows one to explore other discretization scheme to arrive at the LB formulation from the Boltzmann equation Equation (1). In particular, we conjecture that an implicit scheme for the time discretization would improve the error bounds at the cost of the sparsity of the resulting matrix.

3.3. Carleman Linearization of Streaming Step

Instead, we need to consider Equation (1). Another limitation of considering only a traditional one-dimensional model, such as Burger’s is that the problem of streaming the coupled terms of the Carleman linearization is avoided [20]. The forward and backward directions of the velocity are accounted for in one velocity variable with its negative and positive values, and derivatives of the velocity are discretized in terms of that one variable as well.

We have thus far considered linearization of the collision step solely. However, for a self-contained quantum fluid simulation algorithm, Carleman streaming must be considered, which leads us to a modified Boltzmann equation for the Carleman variables:

$$\frac{dV(\overrightarrow{x},t)}{dt}=\frac{\partial V(\overrightarrow{x},t)}{\partial t}+{\overrightarrow{c}}_i·{\mathsf{\Sigma }}_{i=1}^Q\frac{\partial V(\overrightarrow{x},t)}{\partial x_i}=CV$$

(26)

If the partial derivatives are left undiscretized, we have to include additional variables to account for the new terms appearing, contributing to the blowup of variables, and subjecting a description of streaming to truncation. On the other hand, following the same discretization scheme typical of LB, in similarity to the derivation of Equation (2), we arrive at the a modified LB equation where the coefficients of the partial derivatives of the original variables $\overrightarrow{f}$ now show dependence on the variables themselves after replacing the expressions of the Carleman variables $\overrightarrow{V}$ and computing their partial derivatives. Equation (27) shows an example:

$$\frac{\partial V_n(\overrightarrow{x},t)}{\partial f_j(\overrightarrow{x},t)}\frac{\partial f_j(\overrightarrow{x},t)}{\partial t}+{\overrightarrow{c}}_i·\left(\frac{\partial V_n(\overrightarrow{x},t)}{\partial f_j(\overrightarrow{x},t)}\frac{\partial f_j(\overrightarrow{x},t)}{\partial x_i}\right)=C_{n}V$$

(27)

For D1Q3 second order linearization, with $n=4$, we have $V_4=f_1{f}_2$:

$$\begin{array}{cc}\hfill & f_2(\overrightarrow{x},t)(f_1(\overrightarrow{x}+\Delta {\overrightarrow{x}}_1,t)-f_1(\overrightarrow{x},t))+f_1(\overrightarrow{x},t)(f_2(\overrightarrow{x}+\Delta {\overrightarrow{x}}_2,t)-f_2(\overrightarrow{x},t))=0\hfill \\ & f_2(\overrightarrow{x},t)f_1(\overrightarrow{x}+\Delta {\overrightarrow{x}}_1,t)-f_2(\overrightarrow{x},t)f_1(\overrightarrow{x},t)+f_1(\overrightarrow{x},t)f_2(\overrightarrow{x}+\Delta {\overrightarrow{x}}_2,t)-f_1(\overrightarrow{x},t)f_2(\overrightarrow{x},t)=0\hfill \\ & f_2(\overrightarrow{x},t)f_1(\overrightarrow{x}+\Delta {\overrightarrow{x}}_1,t)-2V_4(\overrightarrow{x},t)+f_1(\overrightarrow{x},t)f_2(\overrightarrow{x}+\Delta {\overrightarrow{x}}_2,t)=0\hfill \end{array}$$

(28)

where we see new terms unaccounted for in the Carleman variables appear combining both nonlocality, and nonlinearity.

Another way to explain the above is that the discrete densities of the particles are weighted by their contribution to the Carleman variables, the new variables of the system $\overrightarrow{V}$, when streamed. When discretizing the equation describing the evolution of the Carleman variables, a coupling between terms at different locations appears due to the different partial derivatives appearing in the term. Therefore, classical Carleman linearization of the lattice Boltzmann formulation exchanges the local nonlinearity of the collision step, with nonlocal linearity of the streaming step to which the linearized collision term is coupled. We note that additional V variables must introduced for the nonlocal coupled terms, i.e., $f_1({\overrightarrow{x}}_1+\Delta {\overrightarrow{x}}_1,t)f_2({\overrightarrow{x}}_2,t)$ to keep the system linear. This is indeed used for the Burger’s equation in previous literature [20]. However, this further exacerbates the blowup in variable count for the classical Carleman scheme.

This leads to the fact that streaming is also described by an infinite differential system that must also be truncated. That is, the exactness of streaming, a major advantage of the lattice Boltzmann method, is lost. On a classical computer, to study the collision step, it is possible to recalculate the nonlinear terms to achieve exact streaming. Remarkably, while Carleman linearization slashes out the exact streaming advantage of lattice Boltzmann, we are able to retrieve it by going into the quantum paradigm Itani et al. [26].

3.4. Error Bound

In Equation (8), we see that LB fits into a generalized quadratic ODE considered by [18,20], and it is a multi-population extension of the logistic equation suggested by [20] for treatment.

For a given system of differential equations in terms of a set of variables represented by a vector $\overrightarrow{f}$, of which LB is an example, we denote the solution of the exact system as $\overrightarrow{f}$, and that of the Carleman linearized system as ${\overrightarrow{f}}_C$. We use the same definition of the error $\epsilon$ from [18,20], in terms of the max supremum over time of the vector difference of exact and approximate solution normalized by their respective supremum norms:

$$\epsilon \left(t\right)=\parallel \frac{\overrightarrow{f}\left(t\right)}{\parallel \overrightarrow{f}\left(t\right)\parallel }-\frac{{\overrightarrow{f}}_C\left(t\right)}{\parallel {\overrightarrow{f}}_C\left(t\right)\parallel }\parallel$$

(29)

As we integrate the system of equations, we define the maximum error made as

$${\epsilon }_{max}=\underset{t\in T}{max}\epsilon \left(t\right)$$

(30)

When using Carleman linearization and the Euler time discretization as in Equation (2), they require [20] the minimum number of Carleman variables N will be a function of ${\epsilon }_{max}$ as be:

$$N=-\frac{lo{g}_2\left(R\right)}{lo{g}_2\left(2(1+\frac{1}{{\epsilon }_{max}})\right)}$$

(31)

and the largest time-step $\Delta t$:

$$\Delta t=\frac{1}{N\parallel F_1\parallel }=\left(-\frac{lo{g}_2\left(R\right)}{lo{g}_2\left(2(1+\frac{1}{\epsilon })\right)}\right)\parallel F_1\parallel {)}^{-1}$$

(32)

given that all the eigenvalues of $F_1$ are negative real, and $R<1$.

The above formulas are derived for a single variable quadratic ODE system, $\overrightarrow{f}\in {\mathcal{R}}^1$, for which $Q=1$ and $N=O_c+1$. Furthermore, in [20], it is shown that the dependence of N and $\Delta t$ with the error ${\epsilon }_{max}$ is of the form Equations (31) and (32) even when $R>1$ for the Burgers equation. Note that for the LB problem the bounds of Equations (2), (31) and (32) are not valid as the quadratic system is always multivariate-$Q>1$-and $R>1$ for typical choices of the equilibrium function, as for a single-phase fluid in Equation (8).

Now we prove that the leading order in the error improves exponentially with the Carleman order, for 1,2 and 3D systems.

Let $\left(r\right)$ denote a Carleman variable $V_i$ of rth order, $V_i^{\left(r\right)}$, such that $C_{ij}^{\left(r\right)\left(k\right)}$ describes the matrix coefficient describing the contribution of $V_j^{\left(k\right)}$ to $V_i^{\left(r\right)}$. We have:

$$\begin{array}{cc}\hfill V_i^{\left(1\right)}(\overrightarrow{x},t)=& V_i^{\left(1\right)}(\overrightarrow{x},t-\Delta t)+{\mathsf{\Sigma }}_j{C}_{ij}^{\left(1\right)\left(1\right)}{V}_j^{\left(1\right)}(\overrightarrow{x},t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_k{C}_{ik}^{\left(1\right)\left(2\right)}{V}_k^{\left(2\right)}(\overrightarrow{x},t-\Delta t)\hfill \end{array}$$

(33)

and:

$$\begin{array}{cc}\hfill V_i^{\left(2\right)}(\overrightarrow{x},t)=& V_i^{\left(2\right)}(\overrightarrow{x},t-\Delta t)+{\mathsf{\Sigma }}_j{C}_{ij}^{\left(2\right)\left(1\right)}{V}_j^{\left(1\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_j,t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_k{C}_{ik}^{\left(2\right)\left(2\right)}{V}_k^{\left(2\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_k,t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_l{C}_{il}^{\left(2\right)\left(3\right)}{V}_l^{\left(3\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_l,t-\Delta t)\hfill \\ \hfill =& E_i^{\left(3\right)}(\overrightarrow{x},t-\Delta t)+V_i^{\left(2\right)}(\overrightarrow{x},t-\Delta t)\hfill \end{array}$$

(34)

such that when streaming is considered:

$$\begin{array}{cc}\hfill V_i^{\left(1\right)}(\overrightarrow{x},t)=& V_i^{\left(1\right)}(\overrightarrow{x},t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_j{C}_{ij}^{\left(1\right)\left(1\right)}{V}_j^{\left(1\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_j,t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_k{C}_{ik}^{\left(1\right)\left(2\right)}{V}_k^{\left(2\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_k,t-\Delta t)\hfill \end{array}$$

(35)

Replacing Equation (34) into Equation (35), we have:

$$\begin{array}{cc}\hfill V_i^{\left(1\right)}(\overrightarrow{x},t)=& V_i^{\left(1\right)}(\overrightarrow{x},t-\Delta t)+{\mathsf{\Sigma }}_j{C}_{ij}^{\left(1\right)\left(1\right)}{V}_j^{\left(1\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_j,t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_k{C}_{ik}^{\left(1\right)\left(2\right)}(E_k^{\left(3\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_k,t-2\Delta t)+V_k^{\left(2\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_k,t-2\Delta t))\hfill \end{array}$$

(36)

For the collision-only setup considered for the classically linearized collision term, $\Delta x_i$ = 0, such that one is able to verify:

$${\mathsf{\Sigma }}_j{C}_{ij}^{\left(1\right)\left(2\right)}{E}_j^{\left(2\right)}(\overrightarrow{x},t)=0$$

(37)

for D1Q3, D2Q9 and D3Q27, “full” lattices, for which the expressions become linear in Carleman variables up to second order. Equation (35) then reduces to:

$$\begin{array}{cc}\hfill V_i^{\left(1\right)}(\overrightarrow{x},t)=& V_i^{\left(1\right)}(\overrightarrow{x},t-\Delta t)+{\mathsf{\Sigma }}_j{C}_{ij}^{\left(1\right)\left(1\right)}{V}_j^{\left(1\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_j,t-\Delta t)\hfill \\ & +{\mathsf{\Sigma }}_k{C}_{ik}^{\left(1\right)\left(2\right)}\left(V_k^{\left(2\right)}(\overrightarrow{x}+\Delta {\overrightarrow{x}}_k,t-2\Delta t)\right)\hfill \end{array}$$

(38)

Dropping the E term, and with further replacements similar to above, one can see that the first order terms depend only on the initial conditions of the second order terms in the whole domain, not only neighbouring cells, and no third order terms are needed. This is to say that turning off streaming resolves the nonlinearity of the problem, as expected.

With streaming, a simple Taylor expansion of Equation (35) shows that Carleman linearization of order yields a solution with error of the order:

$${\epsilon }_{max}=O(\Delta t\Delta x^{O_c}).$$

(39)

If we initialize the flow to be uniform, the inclusion of streaming does not introduce errors as above, as Equation (37) still holds for the initial conditions at the cells are identical, and so is the collision step, such that $E(\overrightarrow{x},t)=E(\overrightarrow{x}+\Delta {\overrightarrow{x}}_1,t)=\dots =E(\overrightarrow{x}+\Delta {\overrightarrow{x}}_Q,t)$. At the boundaries, this still holds if they were periodic, but the latter amounts to the trivial case where the kinetic energy of the flow relaxes to zero. The argument for identical evolution for a uniform initial field under periodic boundaries is visualized in Figure 4. Periodic boundaries refer to a fully periodic domain, in all its dimensions, i.e., triply periodic in $3D$, which is typically useful for fundamental studies of homogeneous turbulence.

If (one of) boundaries are not periodic, the error is first generated in the collision step at the boundary, and propagated to the interior of the domain. As illustrated in Figure 5, with each time cycle, the error grows, and it propagates further inside, such that we may speak of a numerical error boundary layer. For example, in a pipe with periodic flow, where the error is first generated at the walls of the pipe, rather than the periodic inlet and outlet.

4. Numerical Results

4.1. Logistic Equation

To demonstrate the utility of Carleman linearization, we consider the logistic equation, as suggested by [20].

$$\begin{array}{cc}\frac{df}{dt}=Kf(1-f)& \forall f\in [0,1]\end{array}$$

(40)

We note that the K factor appears for both first and second order terms, and, thus, cancels out in the calculation of R which remains dependent on the initial conditions solely. In the results, we take $K=1$. We still see an abrupt cut-off in the utility of the linearization for the resulting analytical solution. This is explained by the fact that even though $R\le 1$, the coefficient of the first-order term is definite positive, unity, thereby fulfilling neither $Re\left({\lambda }_1\right)<0$ [20] nor $\mu \left(F_1\right)<0$ [18], the conditions set to gurantee error convergence. According to [18,20], this explains the blow-up in the error shown in the analytical solution. Namely, we restate upper-bounds for time T using the power-series method of [18]:

$$T=\frac{1}{\parallel F_1\parallel }ln(1+\frac{\parallel F_1\parallel }{\parallel f(t=0)\parallel \parallel F_2\parallel })=ln(1+\frac{1}{\parallel f(t=0)\parallel })$$

(41)

which predicts the time of validity for the linearization as shown in

Table 2. The maximum endtime validity for Carleman linearization of the logistic equation as predicted by Equation (41).

$\mathit{f}(\mathit{t}=0)$	T
0	∞
0.1	2.40
0.2	1.79
0.3	1.47
0.4	1.25
0.5	1.10
0.6	0.98
0.7	0.89
0.8	0.81
0.9	0.75
1	0.69

, and with which the results agree, as the analytical solutions presented in Figure 6.

However, for the numerical solution computed through time-discretization of the equation, we note the error shows slower evolution, giving a longer effective time-period to work with, as could be seen with the numerical solutions extending well-beyond their analytical counterparts before blowing up in Figure 6. Namely, we note that $\left|\right|F_2\left|\right|=|Re\left({\lambda }_{MAX}\left(F_1\right)\right)|=K$ such that their ratio is independent of K, but when the equation is discretized $f(T+\Delta t)$ becomes dependent on $f\left(t\right)$ and the original driving function scaled by $\Delta t$, $\Delta K\left(f\right(1-f\left)\right)$. Thus, the choice of K and $\Delta t$ weigh on the contribution of nonlinear terms in the discretized equation, such that we expect more accurate solutions for smaller $\Delta t$ and/or K. This is in line with the findings for the validity of the linearization of the Burger’s equation with $R\approx 40$, and an invitation for more applied, rather than analytical, work in the field. In terms of lattice Boltzmann, which is defined as a discretized problem, this reveals a shortcoming of using the analytical bounds presented for a general quadratic equation. However, given that Carlemann linearization is exact for a collision-only scheme, and the exacerbation of the blowup of variables when streaming is considerd, as discussed in Section 3.3, this is not investigated on a classical computer, and the results presented for a D1Q3 lattice are limited to a collision-only scheme.

4.2. D1Q3

We now concern ourselves with the results of linearizing the collision term of a D1Q3 lattice Boltzmann formulation. As mentioned in methodology section, exact streaming of the linearized system is only possible on a classical computer with the explicit computation of the nonlinear terms, which defies the purpose of a linearization scheme. Therefore, we restrict ourselves to the collision step only. The test case involves two steps, the first is initializing the discrete densities of a single cell, and the second is performing successive collisions.

$${\overrightarrow{f}}^{\left(1\right)}(t+\Delta t)=(I+F_1){\overrightarrow{f}}^{\left(1\right)}+F_2{\overrightarrow{f}}^{\left(2\right)}+{\overrightarrow{F}}_0$$

(42)

where $\tau$ plays a role since we have the contribution of the identity matrix arising from the discretization which does not scale by $\frac{\Delta t}{\tau }$ whereas $F_2$, as well as $F_1$ and ${\overrightarrow{F}}_0$ do for that matter.

The discrete densities are initialized to differ from the weights of the model ${\left(\begin{array}{ccc}\frac{1}{6}& \frac{2}{3}& \frac{1}{6}\end{array}\right)}^T$ by the velocity, chosen to be $u=0.1$. In the results shown in Figure 7, the difference has been distributed amongst the left and right discrete densities, such that their difference is u, ${\left(\begin{array}{ccc}\frac{1}{6}-\frac{u}{2}& \frac{2}{3}& \frac{1}{6}+\frac{u}{2}\end{array}\right)}^T$. We note that this symmetric initialization yields the constant first-order solution shown in blue, but it is not necessary, choosing ${\left(\begin{array}{ccc}\frac{1}{6}& \frac{2}{3}& \frac{1}{6}+u\end{array}\right)}^T$, for example, has also yielded exact solution from second order during our runs. Apart from the initial velocity, we also experiment with different timestep-relaxation time ratios, varying them between $0.1$ and 2, and we see that the exactness holds. The choice of 2 as a maximum ratio is due to the limitations of the lattice Boltzmann scheme under considerations, as would be discussed below. In the absence of streaming, we see in Figure 7 that the solution is exact for all orders of linearization starting from the second.

In terms of computing the Carleman linearization matrices, a different approach is used for the lattice Boltzmann scheme than for the logistic equation. In the presence of a single variable, we know apriori the variables that would appear in the linearization, and their derivatives are easily defined either recursively, in relation to one another, or explictly. However, in the case of the lattice Boltzmann formulation, we have found it more efficient to symbolically express the derivatives of each order, starting with the first, identify the variables that do appear in the expression, and then compute the derivatives of the new ones amongst them. As such, we avoid calculating the derivatives of all monomials, some of which never appear in the formulations, such that we have been able to compute the linearization up to 25th order within a handful of minutes on a regular Intel i7-8750H laptop. Apart from the details of the implementation, setting up the Carlemann linearization matrices allows for parallel in time computation, such that one is able to exploit algorithms to calculate powers of matrices, to precompute the matrix for each timestep, irrespective of the initial data.

Moreover, as could be inferred from Figure 8, with increasing linearization order, the Carleman matrix increases in sparsity and approaches a diagonal matrix, becoming more desirable in terms of numerical stability, and even quantum matrix-inversion algorithms.

Dynamic similarity holds for the lattice Boltzmann model of flow under consideration, such that the Reynold’s number is equal in lattice and physical units. In lattice units, we may write it as:

$$Re=\frac{LU}{\nu }$$

(43)

For the single cell considered, without streaming, there is no characteristic length scale. We may instead define it using the characteristic velocity U and the relaxation constant time $\tau$. The problem setup is that of a flow of decaying velocity, such that we may use the initial velocity u as the characteristic velocity U. Moreover, the kinematic viscosity is defined in lattice units as:

$$\nu =c_S^2(\tau -\frac{1}{2})$$

(44)

It becomes immediately obvious that $\tau >\frac{1}{2}$ is a limitation of the model, or in lattice units, $\frac{\Delta t}{\tau }<2$. Even in the vicinity of 2, the model is considered unstable. However, we note that the Carleman linearized system still tracks the original system exactly beyond the range of stability of the lattice Boltzmann model. If we define the Reynolds number for the grid [27]:

$$R{e}_g=\frac{U\Delta x}{c_S^2(\tau -\frac{1}{2})}=\frac{U}{c_S^2(\tau -\frac{1}{2})}$$

and consider the low-Mach number requirement for the incompressibility assumption, typically $Ma<0.3$, we may arrive at the fact that the grid Reynolds number should be of order $O\left(10\right)$, which in fact, for a given $\tau$, is another way of saying that $\Delta x$ should be chosen in physical units as small enough to capture the long-range structure that arises [27]. As for $\tau$, the recommendation is that it is slightly larger than unity. Put together, they indicate that the lattice Boltzmann scheme under consideration is fit for $ReO(10-100)$. We take the chance here to reiterate the findings that Carleman linearization still exactly tracks the unstable lattice Boltzmann system. Therefore, the limitations on Reynolds number arise from the inherent limitations of the lattice Boltzmann scheme considered. While we expect our findings to extend to other lattice Boltzmann models, this is conditional upon the equations of the extended models and the resulting structure of the coefficient matrices. This paper should be taken as a motivation to investigate other flavours of the lattice Boltzmann method beyond the vanilla flavour considered here.

In terms of Carleman linearization, a $\tau$ larger than unity in lattice units spells good news for the method. It means that the nonlinearity ratio remains bounded by unity even when considering the time discretization:

5. Conclusions

The classical algorithm suffers from a blowup in variable count and sacrifices the exactness of streaming. However, we have shown that the error of the classical Carleman technique could be mitigated, even effaced, in specific applications. On the bright side, we have shown that, at least for the case explored in this paper, the error of the Carleman linearization decreases exponentially with the order of the linearization. This bodes well for the development of a quantum LB algorithm based on Carleman linearization [26].