Quantum Computing for Simulation of Fluid Dynamics

Abstract

The implementation of quantum algorithms for the simulation of classical fluid dynamics poses a fundamental challenge due to the nonlinearity of the fluid equations. In this work, we provide a pedagogical introduction to quantum computing algorithms for simulating classical fluids, with a special focus on the Carleman-Lattice Boltzmann algorithm, which has captured significant attention in the last couple of years. While this algorithm demonstrates satisfactory convergence to analytical solutions for systems at low-to-moderate Reynolds numbers, it also shows an exponential depth of the corresponding quantum circuit. As a result much further analysis is needed to assess the availability of the Carleman-Lattice Boltzmann method on a quantum computer.

1. Introduction

Quantum computing has been vigorously pursued in the recent years, mostly in view of the spectacular speed-up it may provide as compared to the performance of electronic computers, for a number of scientific applications [1]. Quantum computing can be traced back to Richard Feynman’s epoch-making paper, in which he observed that physics “ain’t classical,” hence it ought to be simulated on quantum computers [2]. Following Feynman’s observations (with due credit to previous investigators such as T. Toffoli and E. Fredkin), early theoretical work on quantum computing was performed in the 1980s, for example, Deutsch’s work on the link between quantum theory, universal quantum computers and the Church–Turing principle [3]. Then, with the publication of Shor’s algorithm for integer factoring and Grover’s search algorithm in the middle of the 1990s, the research area gathered significant momentum in terms of theoretical work and quantum computing hardware as well. The research area of quantum computing has continued to grow ever since [4, 5, 6].

In terms of applications, the simulation of quantum many-body systems has received special attention, due to its scientific and industrial relevance, as well as due to its close link with quantum hardware, meaning the possibility of mapping quantum Hamiltonian directly into native quantum gates. In this chapter, we shall focus on a much less beaten track, namely the use of quantum computers for the simulation of classical fluids.

To this purpose, it is convenient to refer to the physics-computing plane defined by the following four quadrants: CC: Classical computing for classical physics; CQ: Classical computing for quantum physics; QC: Quantum computing for classical physics; QQ: Quantum computing for quantum physics as shown in Figure 1 (taken from [7]).

Figure 1.

The four quadrants in the physics-computing plane. The CC and CQ quadrants are the mainstay of current scientific computing. QQ is the quadrant invoked by Feynman, and QC is the quadrant relevant to this paper.

Feynman’s observation pertains to the CQ sector shown in Figure 1, where one is often confronted with exponential complexity barriers associated with the phase-space of quantum many-body problems [8, 9]. The basic idea is that such exponential barriers can be handled by the corresponding exponential capacity of qubit representations offered by the QQ quadrant. In this chapter, we shall focus on the opposite off-diagonal QC quadrant, where the power of quantum computing might be brought to the fruition of hard computational problems in classical physics. In particular, we shall focus on computational fluid dynamics (CFD), which presents a ceaseless quest for better and more efficient algorithms and computers, mostly on account of the problem of fluid turbulence [7, 10]. The physics of fluids presents two general features which are not shared by quantum physics: nonlinearity and dissipation. Both set a major challenge to quantum computing, which draws much of its power from the superposition principle and the unitary and Hamiltonian structure of quantum mechanics. In the following, we shall present a few potential strategies to deal with both issues above.

The chapter proceeds as follows. In Section 2, we introduce the main advantages and issues in dealing with the QC sector. To address some of these challenges, various hybrid classical-quantum approaches have been proposed. These strategies leverage the strengths of both classical and quantum computing, offloading certain tasks to the classical computer while utilizing the quantum computer for operations where it offers a significant speed advantage. In Section 3, we provide a concise review of key hybrid approaches from the literature. While not exhaustive, this review highlights the most prominent approaches. In Section 4, we review some of the fully quantum approaches, that is, algorithms where the computation is completely performed by the quantum computer, and only the extraction of the information is due to the classical computer. Sections 2, 3 and 4 are mostly based on the works of the authors published in Ref. [7]. In Section 5, we present and comment on one of the most promising among these approaches: the Quantum Carleman-Lattice Boltzmann algorithm, originally developed by the authors in Ref. [11]. In Section 6, we draw our conclusions and perspectives.

2. Challenges facing quantum computational fluid dynamics (QCFD)

As mentioned in the introduction, realizing the potential of quantum computing means leveraging distinctive features of quantum mechanics that, by definition, are not available on classical computers. However, it also follows that it is precisely these specific features that expose major challenges in realizing simulations with a quantum advantage.

The main quantum mechanical concepts spawning the potential benefit of quantum algorithms are quantum superposition and quantum parallelism. The quantum state in a Q-qubit coherent register can be described by the Schrödinger wave function, defined by $2^Q$ complex amplitudes for $2^Q$ states in superposition. The square of each of these amplitudes defines the probability of finding the system in the corresponding state after quantum measurement. By encoding classical data in terms of these amplitudes an exponential saving in storage can be achieved when the number of qubits is compared to the required number of classical bits.

Let us illustrate the idea for the specific case of turbulent flow simulation. Turbulence features a ${\mathit{Re}}^3$ complexity, where the Reynolds number $\mathit{Re}$ represents the relative strength of convection (nonlinearity) versus dissipation. Most real life problems feature Reynolds numbers in the many-millions; for instance, an airplane features $\mathit{Re}\sim {10}^8$, implying $O\left({10}^{24}\right)$ floating-point operations per simulation. This is basically the best one can afford on a nearly-ideal Exascale classical computer. The simulation of regional atmospheric circulation flows takes us at least another two decades above in the Reynolds number, hence totally out of reach for any foreseeable classical computer [12, 13]. On the other hand, the minimum number of qubits $Q$ required to represent ${\mathit{Re}}^3$ complexity can be roughly estimated as $2^Q={\mathit{Re}}^3$, namely:

$$Q=3{\mathit{Log}}_2\mathit{\text{Re}}\sim 10{\mathit{log}}_{10}\mathit{Re}$$E1

This simple estimate shows that $80$ qubits match the requirement of full-scale airplane design, while $O\left(100\right)$ qubits would enable regional atmospheric models [12] or nuclear fusion applications [14]. However, several key challenges stand in the way of this task.

First, quantum measurement: Extracting classical information collapses the quantum wave vector into a single eigenstate (a house of cards). Hence, to get classical values for each of the amplitudes, multiple realizations of the quantum state vector are needed with an associated set of measurements.

Second, conditional operations: In the quantum circuit model, the “classical” information is not available to the quantum gate operations performed in the circuit. Specifically, gate operations can be conditional on the state of one or more control qubits, while specifying gate operations conditional on one or more of the complex amplitudes defining the wave function is impossible. Therefore, when classical data are encoded in terms of amplitudes, a rotation angle in a quantum gate operation, this information is required at the time the circuit is compiled. This angle cannot be changed at run time as a function of “classical” data encoded in quantum amplitudes.

Third, time marching: In an algorithm involving multiple iterations or time steps, the overhead associated with quantum measurements used to extract classical data and re-initialization of the quantum state for the next iteration, scale quadratically with the grid size [15] and consequently they severely tax the quantum CFD efficiency. It should be observed that the well-known Harrow-Hassidim-Lloyd (HHL) [16] algorithm for linear system solution assumes that the input and output data are encoded in terms of quantum amplitudes without including the cost of setting up the quantum state and extracting the classical solution.

Fourth, circuit depth: With the exception of quantum measurement operations, quantum mechanical operators are unitary, linear and reversible; hence, they can be directly implemented on a series of unitary quantum gates (quantum circuit model). Usually, in order to quantify the computational cost, producers and specialists take track of the number of two-qubit gates employed in the circuit. In fact, two-qubit gates, such as the CNOT gate, come in the hardware with an error that is several orders of magnitude (depending on the platform) larger than single-qubit gates errors. We know from the DiVincenzo theorem [17] that the number of two-qubits gates needed to implement a generic $Q$ qubits unitary operation scales as $4^Q$, thus voiding the exponential advantage of the embedding process. In order to sidestep this issue, the unitary operation that we want to perform in our circuit has to be combinatorially block diagonal [18]. Hence, this provides a stringent constraint on the possible algorithms that can actually be implemented to simulate classical fluids.

If we further assume that the velocity field is represented in terms of amplitude encoding [19], that is, the components of velocity vector at all grid points are represented in terms of the amplitudes defining the wave function, then the no-cloning theorem prevents the use of (temporary) copies of any of these amplitudes. So, evaluation of $u^2$ or $u\frac{\partial u}{\partial x}$ cannot be performed by storing a temporary copy $\mathit{\text{temp}}=u$, to perform the computation of the value of $u^2$ as $u\times \mathit{\text{temp}}$. Also, for data encoded in terms of the complex amplitudes of the Schrödinger wave function, there is a need for this state vector to have a unit norm, since these amplitudes represent probabilities of states. This means that for an operator attempting to compute the squares of these complex amplitudes, the resulting state vector loses unitarity. So, even without the no-cloning theorem complicating such a step, this points to a further problem with computing nonlinear terms. This loss of unit norm of the quantum state vector represents an example of a non-unitary operation, typically associated with a corresponding loss of information. A similar argument runs for orthogonality of the eigenstates, since a nonlinear propagator rotates the state by an angle which depends on the state itself.

Fifth, dissipation: Dissipation breaks hermicity of the quantum propagator. This can be recovered by augmenting the system with its Hermitian conjugate, so that the doubled system is Hermitian by construction. One can dispense with doubling degrees of freedom by representing the non-unitary update as a weighted sum of unitaries. This introduces a failure probability which needs to be handled probabilistically, that is, by executing the update conditional on the successful measurement of an ancilla qubit.

To summarize, dealing with nonlinear and dissipative terms raises major extra challenges for quantum computational fluid dynamics (QCFD) besides the ones commonly encountered for the simulation of quantum systems.

3. Hybrid quantum/classical approaches

The challenges sketched above have resulted in the fact that most of the existing work in QCFD is based on a hybrid quantum/classical approach, with the quantum processor performing computations for which efficient quantum algorithms exist, while the output is then passed on to classical hardware to perform further computational tasks not (yet) amenable to quantum algorithms.

Figure 2 illustrates the idea. The quantum state ${\psi }_0$ is advanced to the next quantum state ${\psi }_1$ via a QQ algorithm. The quantum state ${\psi }_1$ is then measured to generate classical observables $C_1$ which are advanced to $C_2$ by a CC algorithm. The classical observables $C_2$ are then used to reconstruct the quantum state ${\psi }_2$, ready for the next QQ step. The Q2C conversion shown in Figure 2 involves quantum measurements and requires repeated measurements over a statistical sample of quantum states, since none of them can be re-used. The C2Q reconstruction requires the preparation of all the quantum eigenstates, hence a full reset of the quantum circuits from scratch. Both operations impose a substantial computational burden on the hybrid algorithm, typically scaling exponentially with the number of qubits.

Figure 2.

Sketch of a hybrid quantum-classical algorithm. The illustration shows the steps involved in a single time step or single iteration, including preparation of the subsequent step/iteration.

Examples of previous works using the hybrid quantum/classical approach include the works of Steijl [20], Gaitan [21] and Budinski [22]. The algorithm presented by Gaitan uses Kacewiz’s quantum amplitude estimation ODE algorithm [23] as applied to the set of nonlinear ODE’s resulting from standard discretization of the Navier-Stokes equations. As shown, for certain “non-smooth” problems (illustrated using the quasi-1D flow in converging-diverging duct with normal shock wave), the complexity analysis shows potential for exponential speed-up, so that the challenges associated with hybrid quantum/classical computing can potentially be overcome. For the linear advection-diffusion equation, Budinski [24] presented a quantum algorithm based on the lattice Boltzmann method [25, 26]. The algorithm can perform multiple successive time steps with no need for quantum measurements and re-initialization of qubit register between the time steps if suitable re-scaling of solution vector is applied to deal with non-unitarity. In the quantum circuit model, the fact that velocity field is unchanged between successive time steps implies that the operator implementing $u\frac{\partial u}{\partial x}$ can be re-used in multiple time steps. Budinski [22] then extended the work to Navier-Stokes equations in stream function-vorticity formulation. Then, velocity field updates during each time step means that quantum circuit implementation of convection terms cannot be re-used during multiple time steps and that the classical value of $u$ (as well as other flow field data) is needed to define quantum circuit implementation for the next time step. This shows that it is the nonlinearity that forces the use of a hybrid quantum/classical approach, similar to previous work where Navier-Stokes equations were discretized, for example, the hybrid quantum/classical algorithms for fluid simulations based on quantum-Poisson solvers [20].

In summary, key challenges for the hybrid quantum/classical approach are: (i) cost and complexity of (repeated) measurements; (ii) statistical noise due to sampling of the quantum solution; (iii) cost and complexity of (repeated) re-initialization.

The development of more efficient (re-)initialization techniques appears to be a key factor for the future success of quantum computing for fluids.

4. Quantum fluid dynamics strategies

For the sake of concreteness, let us refer to the Navier-Stokes equations for time-dependent incompressible flows:

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \frac{\partial \mathbf{u}}{\partial \mathbf{x}}=-\frac{\partial P}{\partial \mathbf{x}}+\nu \mathrm{\Delta }\mathbf{u}$$E2

$$\nabla \cdot \mathbf{u}=0$$E3

where $\mathbf{u}$ is the velocity vector, $P$ is the pressure, defined for location $\mathbf{x}$ as function of time $t$. The kinematic viscosity is defined by $\nu$ and density has been conventionally set to a unit constant value. Eq. (3) enforces mass conservation, while Eq. (2) is based on momentum conservation in each coordinate direction. Eq. (2) and Eq. (3) highlight that it is the convection term that represents the nonlinearity, that is, the second term on the left-hand side of Eq. (2). Writing the Navier-Stokes equations in non-dimensional form, such that $\mathbf{x}$ and $\mathbf{u}$ are scaled by reference length $L_{\mathit{ref}}$ and $U_{\mathit{ref}}$, respectively, it follows that in Eq. (2) the term $\nu$ is replaced by $1/\mathit{Re}$, where Reynolds number $\mathit{Re}=U_{\mathit{ref}}{L}_{\mathit{ref}}/\nu$. For Stokes flow, that is, with Reynolds approaching $0$, nonlinear terms are vanishingly small, but still not to be neglected since they are responsible for non-trivial long-range correlations especially important in biological flows [46]. For high Reynolds number (turbulent) flows, obviously the nonlinear terms play the leading role.

Different strategies have been adopted to simulate fluid dynamics on quantum computer, mostly involving hybrid algorithms. Since they have been reviewed in the recent perspective [7], PNAS [27], and their time complexity has been analyzed thoroughly in [28], here we simply list them, directing the reader to the original literature.

4.1 Nonlinear quantum ODE solvers

This approach consists in treating the discretized Navier-Stokes equations as a set of nonlinear ODEs and advancing them in time by developing bespoken quantum nonlinear ODE solvers. This approach has been pioneered by Gaitan [21] and demonstrated for the case of a Laval nozzle on grids with $O\left(30-60\right)$ grid points over $1400$ time steps. The critical issue, though, is that the concrete implementation of the quantum oracle is left to a classical computer.

4.2 Nonlinear variational quantum eigenvalue solvers

It has recently been proposed that variational quantum eigenvalue (VQE) solvers, a major tool of the QQ sector, might be extended to the fluid equations as well [29]. The basic idea is to use quantum computers for the efficient generation of variational eigenfunctions by suitable circuit parametrization. The associated energy functional is then minimized through a classical procedure. The appeal of these methods is the possibility to import algorithmic know-how from quantum mechanical applications, especially quantum chemistry. To date, however, there is no evidence that the procedure can meet the standards of accuracy required by computational fluid dynamics.

4.3 Functional approach

It is long known that any nonlinear problem can be mapped into a linear one in a space of higher dimensions, a technique also known as Carleman embedding or Carleman linearization [30]. The Carleman embedding trades nonlinearity for infinite dimensions and nonlocality. An alternative procedure to recover linearity is to take a probabilistic approach and formulate a functional kinetic (Liouville) equation for the probability distribution function (PDF) of the fluid velocity field. Formally [31]:

$${\partial }_{t}P\left[\mathbf{u}\right]+\frac{\delta }{\delta \mathbf{u}}\left(f\left(\mathbf{u}\right)P\left[\mathbf{u}\right]\right)=0$$E4

where $P=P\left[\mathbf{u}\right]$ is the functional PDF and $\dot{\mathbf{u}}=f\left(\mathbf{u}\right)$ is the equation of motion (hence $f\left(\mathbf{u}\right)$ is a nonlinear and non-local operator in function space). This is linear by construction, regardless of the nonlinearity of the dynamics. Once the fluid equations are discretized on a grid with $G$ lattice nodes, the Liouville distribution $P_G$ becomes a $3G$-variate PDF $P_G\left({\mathbf{u}}_1\dots {\mathbf{u}}_G\right)$ which lives in a $O\left(3G\right)$-dimensional space. For large-scale flow applications, $G$ can reach values up to multi-billions, hence fully under the “curse of dimensionality.”

Fortunately, for most practical purposes, low-order marginals prove often capable to retain the essential physical information. Marginalization is known to introduce coupling terms between different marginals, thus raising the need for appropriate closures expressing the high-order marginals as functionals of the low-order ones. This is a classical topic in non-equilibrium statistical physics, which may draw significant benefits from modern developments in tensor networks theory [32, 33].

The idea can be taken one level further by noting that Eq. (4) is linear, but still not unitary. In Ref. [34] the authors presented a method to relate the classical evolution and the unitary evolution of a quantum system, making use of the Carleman-Koopman embedding, which transforms the Liouville equation into an equivalent many-body Schrödinger equation.

This is formally very appealing, but still subject to the dimensional curse; even in quantum chemistry nobody knows how to handle molecules with billions atoms with ab initio methods. In addition, Giannakis et al. [34] show that for chaotic systems, the spectrum of the Hamiltonian contains continuous modes that are pretty hard to tame.

Always in the spirit of functional methods, in Ref. [28], the authors develop another form of Koopman-von Neumann approach based on the level-set method as applied to classical nonlinear field theories [28]. The formalism is applied to hyperbolic PDE’s and Hamilton-Jacobi equations, but its extension to the Navier-Stokes equations is addressed only marginally, making it difficult to draw firm conclusions.

4.4 Carleman embedding

Introducing indices $m$ and $n$ to denote interacting neighbors, and using Einstein summation, the fluid equations can be recast in terms of a linear operator $\mathcal{L}$ and quadratic operator $\mathcal{Q}$, as follows,

$$\frac{{\mathit{du}}_l}{\mathit{dt}}={\mathcal{L}}_{\mathit{lm}}{u}_m+\mathit{Re}{\mathcal{Q}}_{\mathit{lmn}}{u}_m{u}_n$$E5

The pressure term was removed for simplicity, although this all but a minor item, since pressure-flow-stress coupling is going to affect the structure of the Carleman matrix.

Discretizing velocity derivatives introduces a key challenge: Carleman linearization creates a growing hierarchy of variables ($U_{\mathit{lm}}=u_l{u}_m$). Each level’s tensor has $O\left(G{\kappa }^{k-1}\right)$ components, $\kappa$ being the sparsity of the $\mathcal{Q}$ matrix, and spreads across a larger spatial neighborhood. This nonlinear uplift incurs dimensionality and locality costs.

Despite this, Liu et al. [35] presents a promising quantum Carleman linearization algorithm for the 1D Burgers equation showing a polylog scaling with the number of grid points, that is an exponential improvement over classical approach. However, for a given time span $T$, the algorithm shows a complexity $T^2\mathit{\text{Poly}}\left(\mathit{\text{logT}}\right)$. The authors simulated shock formation at surprisingly high Reynolds numbers (up to 40), exceeding theoretical limits and warranting further analysis.

The key question then is: Is there a maximum Reynolds number above which the Carleman procedure fails to converge? Also, when it does converge, how many Carleman levels are required to reach the desired accuracy at a given Reynolds number?

The answer may well depend on the chosen representation of the fluid equations and in the following we turn to a promising candidate in this respect, namely the Lattice Boltzmann method [25, 36].

5. Carleman-lattice Boltzmann

In a recent paper [37], Cheung and co-workers argued that Lattice Boltzmann might be more convenient for Carleman linearization, since the fluid nonlinearity is formally encoded in the Mach number instead of the Reynolds number, which makes of course a huge difference, as the Mach number is typically many orders of magnitude smaller than the Reynolds number (for an airplane $\mathrm{Ma}\sim 1$ vs. $Re\sim {10}^8$). They performed a classical Taylor-Green Vortex (TGV) simulation based on a Carleman-Lattice Boltzmann scheme, showing excellent agreement at moderate Mach number, with just three Carleman levels [37]. It is an unsteady flow of a decaying vortex which has an exact closed form solution of the incompressible Navier-Stokes equations. The initial conditions for the three velocity components are:

$$u_x=A_{x}cos\left(k_{x}x\right)sin\left(k_{y}y\right)sin\left(k_{z}z\right),$$E6

$$u_y=A_{y}sin\left(k_{x}x\right)cos\left(k_{y}y\right)sin\left(k_{z}z\right),$$E7

$$u_z=A_{z}sin\left(k_{x}x\right)sin\left(k_{y}y\right)cos\left(k_{z}z\right).$$E8

The TGV is a crucial benchmark in CFD because it allows a direct comparison with a known analytical solution. More recently, in Ref. [11], the authors showed the convergence of the Carleman method to the analytical solution of a Kolmogorov-like flow, a close 2D analogue of the TGV with periodic initial conditions:

$$u_x=A_{x}cos\left(k_{x}y\right)$$E9

$$u_y=A_{y}cos\left(k_{y}x\right).$$E10

It was shown that the Carleman system truncated at second order for Moderate-to-low Reynolds number, up to 100 (see Figure 3) provided excellent agreement up to few hundred time steps.

Figure 3.

(a) The deviation $\epsilon$ in time of the second order Carleman solution from the original LB solution for the case of a Kolmogorov-like flow at $\mathit{Re}\sim 10$, on a $32\times 32$ grid. (b) The mean deviation $⟨\epsilon ⟩$ as a function of the Reynolds number for the same flow.

Let us remind that the Lattice Boltzmann (LB) method is based on the following scheme: In equations:

$$f_i\left(\mathbf{x}+{\mathbf{v}}_i\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)-f_i\left(\mathbf{x},t\right)=-\frac{f_i-f_i^{\mathit{eq}}}{\tau }$$E11

where $f_i\left(\mathbf{x},t\right)$, is the probability to find a representative fluid parcel with discrete velocity ${\mathbf{v}}_i$ at position $\mathbf{x}$ and time $t$. In the above, $f_i^{\mathit{eq}}$ is the corresponding discrete local equilibrium and $\tau$ is a local relaxation time, fixing the viscosity of the lattice fluid. Importantly, the local equilibrium is a quadratic function of the Mach number $\mathit{Ma}=u/c_s$, $c_s$ being the speed of sound. A defining feature of the LB method is that it allows to use a fixed, generally small, number of discrete velocities ${\mathbf{v}}_i$. For instance, it turns out that just nine discrete velocities are sufficient to fully recover the physical properties of a two-dimensional flow defined on a square lattice. The nine velocities are depicted in Figure 4 and described in Table 1. The model is called D2Q9. A key feature of LB, not shared by the macroscopic representation of fluid dynamics, is that streaming proceeds along straight lines, defined by the discrete velocities. This operation is exact, in that it does not involve any floating-point calculation but just a memory shift from the lattice position $\mathbf{x}$ to another lattice position $\mathbf{x}+{\mathbf{v}}_i\mathrm{\Delta }t$. Also to be observed that collisions are completely local, hence perfectly amenable to parallel computing. Because of this, dissipation is an emergent property which does not require any second-order spatial derivative, as it is the case for the Navier-Stokes equations. Besides their mathematical elegance, these properties are at the roots of the computational efficiency of the method and its amenability to parallel computing.

Figure 4.

The D2Q9 model, defined on a two-dimensional square lattice, with the corresponding set of nine velocity vectors labelled by index $i$, cf. Table 1.

i	0	1	2	3	4	5	6	7	8
$c_i$	$\left(0,0\right)$	$\left(1,0\right)$	$\left(0,1\right)$	$\left(-1,0\right)$	$\left(0,-1\right)$	$\left(1,1\right)$	$\left(-1,1\right)$	$\left(-1,-1\right)$	$\left(1,-1\right)$

Table 1.

The velocity set of the D2Q9 model.

A key point of using the discrete-velocity Boltzmann formalism instead of Navier-Stokes is that, owing to the double dimensional phase-space, in the latter nonlocality (streaming) is linear while nonlinearity (collision) is local, while in Navier-Stokes the two merge into a single $\mathbf{u}\nabla \mathbf{u}$ convective term.

On classical computers, this disentanglement proves extremely beneficial because of the local nature of the collision term, which makes the system highly amenable to parallel computation. On quantum computers instead it is not simple to relate together the streaming and the collision steps, as we are going to show next. When writing the Lattice Boltzmann functions $f_i\left(x,t\right)$ in a quantum state we can employ amplitude encoding, using a quantum register $\mid \cdot ⟩{}_v$ for embedding the discrete velocities $v_i$ and a second quantum register for encoding the lattice sites $\mid \cdot ⟩{}_x$. The fluid is then described by the following quantum state

$$\mid \psi \left(t\right)⟩=\sum _i\sum _x{f}_i\left(x,t\right)\mid i⟩\mid x⟩.$$E12

Since the collision step is local and depends only on the velocities of the fluid in the lattice site, one may apply an operator $\widehat{\mathrm{\Omega }}$ which performs the collision step on the quantum register adapted to the velocities and another operator $\widehat{S}$ encoding the streaming process to be applied to all quantum registers, being it both non-local and velocity dependent. This ideal circuit is shown in Figure 5(a). It turns out that since $\widehat{\mathrm{\Omega }}$ is a nonlinear and non-unitary operator, this circuit cannot be realized, unless specific symmetry conditions hold, such as those used in Ref. [38]. Unfortunately, these do not conform to the universal quantum computation paradigm.

Figure 5.

The sketch of three quantum circuits for different lattice Boltzmann implementations for a grid with $G$ lattice sites for the D2Q9 model. In (a) the collision is applied on the velocity register, whereas the streaming is a global operation applied on both the position and the velocity registers. In (b) the evolution is obtained as a series of streaming in the eight velocities of the D2Q9 model. In (c) the quantum register $\mid k⟩$ holds the information about the truncation order, whereas $\mid i⟩$ and $\mid j⟩$ are the velocities associated to the Carleman variables.

We may then opt for hybrid strategies whereby only streaming (collisions) would be performed on a quantum computer, leaving collisions (streaming) to a classical one. Interestingly, this leads to different circuits, as sketched in Figure 5. Figure 5b and c shows the circuits for collision-free streaming and streaming-free collision processes, that we are going to discuss next.

5.1 Collision-free streaming

The quantum computing algorithm for streaming is based on the approach used by Steijl and co-workers [39]. The streaming is a non-local process shifting the particles in the next neighboring site, depending on the corresponding velocity, without affecting their value. Thus, this step can be executed by moving the probability function $f_i\left(\mathbf{x},t\right)$ to $f_i\left(\mathbf{x}+{\mathbf{v}}_i,t+1\right)$ at the next time step.

The corresponding operation on the quantum state reads as follows:

$$\mid \psi \left(t\right)⟩=\sum _i\sum _{\mathbf{x}}{f}_i\left(\mathbf{x},t\right)\mid i⟩\mid \mathbf{x}⟩\stackrel{\widehat{S}}{\to }\mid \psi \left(t+1\right)⟩=\sum _i\sum _{\mathbf{x}}{f}_i\left(\mathbf{x},t\right)\mid i⟩\mid \mathbf{x}+{\mathbf{v}}_i⟩$$E13

That is, only the positions register is affected by the streaming process, whereas the velocity register acts as a control. The streaming operator is obtained as a combination of the streaming in the different directions, for the D2Q9 model with eight velocities (plus the null velocity) is $\widehat{S}={\oplus }_{i=1}^8{\widehat{S}}_i,$ that are therefore applied depending on the value of the qubit in the velocity register. The streaming can be obtained by applying a series of multiqubit controlled operations on the position register. We show in Figure 6 an example for the ${\widehat{S}}_1$ streaming operator when applied on a $8\times 8$ lattice.

Figure 6.

The quantum circuits for the operators $S_1$ and $S_2$, cf. with Table 1, for a D2Q9 model define on a $8\times 8$ lattice. The first three qubits are used to encode the position in the $x$ direction and the latter three for the position in the $y$ direction. They are applied to the quantum state only if the velocity register has value 1 or 2 respectively.

5.2 Streaming-free collision

A completely different approach is employed to perform the collision step: Carleman linearization. Based on the Lattice Boltzmann method, Itani et al. [40] developed an approach termed Carleman for second-quantized Lattice Boltzmann. This terminology stems from the fact that the Boltzmann operator, defined by:

$$\mathcal{B}{f}_i=-{\mathbf{v}}_i\cdot \nabla f_i-\frac{f_i-f_i^{\mathit{eq}}}{\tau },$$E14

can be expressed in terms of the second-quantization annihilation and generation operators via the relation $\nabla =\left(\widehat{a}-{\widehat{a}}^{+}\right)$. Collisions follow the bosonic encoding first proposed by Mezzacapo et al. [41], whereby dissipative effects are represented as a weighted sum of two unitary operators. The scheme is as “Feynmansque” as it can possibly get, since it builds on a one-to-one analogy between LB and the Dirac equation, as first proposed in Ref. [42]. As a result, the formal solution $f_t=e^{\mathcal{B}t}{f}_0$ can be computed in analogy with quantum mechanics. However, it is subject to a number of additional questions: primarily truncation effects due to the finite number of bosonic excitations and the long-time behavior of non-unitarity errors [43].

On a similar but operationally different line, Sanavio et al. [11] developed a quantum circuit implementing the collision step of a Carleman-Lattice Boltzmann algorithm. Since the number of variables with the Carleman method increases exponentially with the truncation order $k$, the circuit needs more qubits to embed the relevant information. The collision step requires $k+1$ additional quantum registers, embedding the truncation order and the products between local Boltzmann functions with different velocities, say $f_i\left(x\right),f_i\left(x\right)f_j\left(y\right),f_i\left(x\right)f_j\left(y\right)f_k\left(z\right)$ and so on. With this setting, the collision quantum circuit is genuinely local, as it depends only on the velocity registers and it can be implemented using only a fixed number of two-qubit gates, regardless of the number of lattice sites. This meets the promise of quantum advantage but unfortunately it is inconsistent with the streaming process, which must therefore be directed to a classical machine. This may seem an ideal option for hybrid quantum computing, since streaming is a floating-point free operation. However, accessing data is by no means a cost-free operation and it is indeed known to take a significant fraction of the computer time spent on collisions [44].

5.3 Fully quantum algorithm: collision and streaming

In Ref. [11], the authors managed to generalize the circuit for the streaming step to Carleman variables of higher orders, thus allowing the subsequent application of streaming and collision step on Carleman variables with a single embedding and a single readout of the results after evolution. This results in a fully quantum algorithm. In order to embed different nonlinear terms one needs $2\left(k+1\right)$ additional quantum registers, one encoding the information about the degree of non-linearity up to order $k$, and the other $2k$ to embed position and velocity of the multiproduct functions.

This method was found to meet with an exponential depth of the quantum circuit as a function of the number of qubits. The reason is that the collision operator cannot be written in the form of a sparse and combinatorially block-diagonal matrix, nor in terms of a compact combination of Pauli matrices. Hence, to date, it cannot be quantum compiled.

Possible ways out are currently under investigation: one option is the application of the Carleman linearization method directly to the Navier-Stokes equations. The advantage of the Carleman-Navier-Stokes procedure rests in a much lesser number of Carleman variables, with a corresponding reduction of the circuit depth. The flip side, however, is a more complex structure of the Carleman matrix, which, combined with the Reynolds versus Mach number as a measure of nonlinearity, may impair Carleman convergence. A detailed analysis along this direction is currently underway.

6. Conclusions

In summary, we have presented a survey of the few leading approaches to the quantum simulation of classical fluids. Various obstacles stand in the way of the efficient simulation of fluid flows on quantum computers, especially at high Reynolds numbers, mostly on account of the strong nonlinear effects.

In this chapter, we focused on the Quantum Carleman-Lattice Boltzmann algorithm. Our analysis, informed by prior works [7, 11, 37, 40, 43], indicates that Carleman linearization offers a promising route for modeling low-to-moderate Reynolds number flows, achieving errors on the order of ${10}^{-3}$ for $Re$ in the hundreds up to hundred of time steps. However, a significant limitation lies in the exponential depth of the corresponding quantum circuit. Future research should focus on optimizing quantum CLB’s resource requirements, potentially through alternative linearization strategies [45] or tailored quantum circuit design. This work highlights the potential of quantum algorithms for fluid simulations, while underscoring the crucial need for addressing their computational demands to unlock real-world computational fluid dynamics applications.

Before closing, it should be observed that while, everybody would expect quantum computers to handle turbulence, the physics of fluids is littered with interesting problems at low-Reynolds, especially in microfluidics, soft matter and biology [36, 47, 48], but also in high-energy physics, typically quark-gluon plasma experiments. Hence, even in case turbulent flows would prove beyond reach to quantum computing, there would nevertheless be plenty of interesting fluid applications in the low-moderate Reynolds regime.

For instance, it could be of great interest to devise a quantum multi-scale application, coupling quantum algorithms for biomolecules in a water solvent described by a quantum algorithm for low-Reynolds fluid flow. Given that low-Reynolds flows are non-local, perhaps the inherent nonlocality of quantum mechanics could prove helpful in representing the classical nonlocality of low-Reynolds flows.

On a more philosophical ground, should the physics of fluids show a case for “classical advantage,” there would still be interesting lessons to learn. Indeed, as observed earlier on, while it is true that nature is not classical, it is equally true that nature has a very strong tendency to become classical at macroscopic scales/high temperatures. Withstanding such a tendency, or perhaps even taking advantage of it, is indeed the prime struggle of quantum computers. Looking into ways to achieve this very hard goal is not only of practical but also of major foundational interest, since it is likely to shed new light into the still open problem of “classicalization,” that is, the emergence of classical behavior out of underlying quantum mechanical microscopic physics. This is only one (outstanding) example of the fact that quantum computers offer the unique chance to ask questions that the founding fathers of quantum mechanics could only formulate as “Gedanken Experiments.”

Acknowledgments

The authors have benefited from valuable discussions with many colleagues, particularly P. Coveney, N. Defenu, G. Galli, M. Grossi, B. Huang, W. Itani, A. Mezzacapo, S. Ruffo, A. Solfanelli K. Sreenivasan, R. Steijl and T. Weaving. The authors also acknowledge financial support from the Italian National Centre for HPC, Big Data and Quantum Computing (CN00000013).