Three Carleman routes to the quantum simulation of classical fluids

We discuss the Carleman approach to the quantum simulation of classical fluids, as applied to i) Lattice Boltzmann (CLB), ii) Navier-Stokes (CNS) and iii) Grad (CG) formulations of fluid dynamics. CLB shows excellent convergence properties, but it is plagued by nonlocality which results in an exponential depth of the corresponding circuit with the number of Carleman variables. The CNS offers a dramatic reduction of the number Carleman variables, which might lead to a viable depth, provided locality can be preserved and convergence can be achieved with a moderate number of iterates also at sizeable Reynolds numbers. Finally it is argued that CG might combine the best of CLB and CNS.


I. INTRODUCTION
One of the authors (SS) had the good fortune to be in friendly relations with Sreeni for over three decades, a great scientific honor and a human privilege alike.Yet, we never coauthored any paper until last year and on a topic which is no mainstream to either of us: quantum computing, and more specifically, quantum computing for fluids.Hence, on the occasion of his Festschrift, the same author resolved to write about this fascinating and over-challenging topic.
Quantum computing (QC) offers tantalizing prospects for cracking problems which lie far beyond reach of any foreseeable classical computer, typical examples in point being quantum many-body sytems as they occur in quantum material science and quantum chemistry 9,21,23 .
The main point rests with the very peculiar property of quantum system to live in a linear superposition of states, each of which carries independent information, along with the capability of processing it, independently.In a nutshell: built-in quantum parallelism.The potential is huge: the Hilbert space of a N-bit quantum system contains 2 N quantum states, which can be stored and processed in the form of a string of N qubits, i.e units of quantum information which, unlike classical bits, can take any value between 0 and 1.As a result, quantum computing offers exponential advantage over classical computing.
If this looks too good to be true, it is because it is indeed too good to be true: many obstacles stand on the way of the aforementioned blue-sky quantum computing scenario.First, decoherence: in order to realize the above potential qubits need to be entangled, meaning that any action on qubit A (Alice) is going to affect qubit B (Bob) as well, and vice-versa.Entanglement is a subtle form of correlation which proves very fragile against external perturbations, the environment, the result being that entanglement decays very fast, currently in a time lapse of hundreds of microseconds.At a processing speed of, say, one quantum update per nanoseconds, this means a few hundred thousands operations, before the qubit "dies out".
The second technological problem is quantum noise: even if the qubit is "alive", this does not mean that it computes errorfree.The same is true for classical bits, the difference being that the error rates is some fifteen orders of magnitude larger, something like 10 −3 against the 10 −18 rate of classical computers!Recovering such gigantic gap on mere technological grounds appears rather desperate, but error correction/mitigation algorithm may offer a robust helping hand to this purpose.This is why, despite the daunting barriers above, research on quantum computing is still burgeoning on both hardware and software fronts, over thirty years after the basic idea was first brought up 8 .
As mentioned above, the best candidates for quantum computing are quantum many-body systems (we focus on natural science, leaving aside paramount applications such as cryptography).Yet, it is only natural to wonder whether QC may contribute to solve also hard problems in classical physics, turbulence being a prominent example in point 3,12,29 (general gravitation would come next).That's exactly the topic where Sreeni and SS finally managed to publish our first paper 31 after three decades of unpublished scientific exchanges!

II. QUANTUM COMPUTING FOR FLUIDS
It is often heard that turbulent flows raise a ceaseless demand for increasingly more powerful computers and improved computational methods.The main culprit being the fact that the number of active degrees of freedom scales like the cube of the Reynolds number, and since Reynolds number in Nature easily exceeds millions (automobiles), billions (regional weather forecast) and trillions (astrophysics), no foreseeable classical computer can meet this demand.Quantum computing has potential to put this quest at rest 2 .
Indeed, the blue-sky scenario potential of quantum computing for fluids is mind-boggling: given that a turbulent fluid at Reynolds Re features Re 3 dynamic degrees of freedom, the 0 0.4 0. number of qubits required to represent it is given by 25,27 q = 3Log 2 Re ∼ 10Log 10 Re This means that a full airplane simulation (Re ∼ 10 8 ) takes just q ∼ 80 qubits, well within the nominal capabilities of current quantum hardware 11 .At the top of the line (as of 2024), q ∼ 500, would formally enable simulations at Re ∼ 10 500 which is far beyond any conceivable need in the physics of fluids!Again, too good to be true: the 500 nominal (physical) qubits must be mapped back to the actual number of effective (logical) qubits, whose ratio is often estimated at about 1 : 1000.Even taking a more optimistic 1 : 100 ratio, the Exascale bar, which we place conventionally around Re ∼ 10 8 , would require 8000 physical qubits, about an order of magnitude beyond the current quantum hardware capabilities.
In the previous section we have mentioned the two major technological stumbling blocks for quantum computers: decoherence and quantum noise.When we move from computers to computing, other major limitations take stage.The first is that not every problem can be formulated in terms of an efficient quantum algorithm, meaning by this an algorithm which can be (efficiently) mapped into the circuits of a quantum computer.Such circuits consisted of a collections of one and two qubit gates performing linear and unitary operations, in full compliance with the linear and unitary nature of quantum mechanics.This alone signals two major elephants in the room for quantum computing of fluids systems: quantum mechanics is linear, unitary and hamiltonian, the physics of fluids is generally none of the three.
Several strategies can be conceived to turn around these problems, but the one that Sreeni and SS have been exploring together is the so-called Carleman linearization (or embedding), CL for short.

III. CARLEMAN LINEARIZATION
The general idea of the CL procedure is to turn a nonlinear finite-dimensional system into an infinite-dimensional linear one.Let us illustrate the point with specific reference to the simplest and yet representative example of the logistic equation: where R, the competition rate, is a measure of the nonlinearity.The exact solution is readily derived and reaches its stable time asymptotic value x ∞ = 1/R for t ≫ 1.The Carleman procedure consists in renaming x (1) ≡ x and x (2) ≡ x 2 , so that the logistic equation takes the linear form Iterating the procedure to the k-th order delivers: This is an endless hierarchy, which is then truncated at a given order K by setting x K = 0, in the hope that the truncated solution captures the essential behaviour of the exact one.In practice, it can be shown that the Carleman hierarchy is just another way of representing the exact solution as an infinite power series in the saturating term Rx 0 (e t − 1).Each successive iteration prolongs the time-horizon of convergence, see Fig. 1, beyond which convergence is abruptly lost.
It should be observed that convergence is controlled by the parameter Rx 0 , namely the ratio x 0 /x ∞ between the initial and time asymptotic values.
Despite its simplicity, the logistic equation delivers a few useful hints for the Carleman linearization of fluids.First, it is a quadratic non-linearity and second it shows that the timehorizon of the Carleman series depends on the strength of the nonlinearity via the ratio x 0 /x ∞ .Third, the logistic equation bears a closed similarity to the collision operator of kinetic equations.This observation prompted out the earliest attempts to apply the Carleman procedure to the Lattice formulation of fluids, as we describe next.

A. Carleman Lattice Boltzmann
The Carleman lattice Boltzmann (CLB) procedure has been first advocated in 13 and lately further explored by a number of authors 5,6,16,26 .The LB equation 1 takes the following form where f i ≡ f i (⃗ x,t) is the probability of finding a fluid parcel ("population") at position ⃗ x and time t with discrete velocity ⃗ c i and i = 0, . . ., n v , chosen from a suitably discrete lattice.In the above f eq i is a local-quadratic function of the flow field J a (⃗ x,t) = ∑ i c ia f i (⃗ x,t)/ρ, where ρ(⃗ x,t) = ∑ i f i (⃗ x;t) is the local fluid density.In equations: where latin indices label spatial components and the discrete velocities are rescaled with the sound speed ⃗ c i → ⃗ c i /c s .The parameter ω is inversely proportional to the relaxation time of the fluid and to the viscosity ν, which can be obtained through the formula .
In order to deal with the nonpolynomial term ρ −1 J a J b we apply the weakly-compressible fluid approximation, ρ ∼ 1, substituting ρ −1 ≈ 2 − ρ.This makes the equilibrium functions f eq i cubic polynomials of the LB distribution functions.The LB method has met with spectacular success for the simulation of fluid problems across a broad spectrum of applications, scales and regimes 19,28 .Its main virtues are as follows: First, free-streaming proceeds along discrete characteristics ⃗ ∆x i = ⃗ c i ∆t and it is floating-point free, for it amounts to shifting the populations from the source site ⃗ x to the corresponding lattice destination ⃗ x i = ⃗ x +⃗ c i ∆t.Since streaming always lands the population on lattice sites, it is exact on a computer and free of most lattice artifacts.Second, equilibria are nonlinear (quadratic) but local, the result being that dissipation is an emergent property which does not require second order spatial derivatives.
When it comes to quantum computing, a third point pops out: the nonlinearity of the LB method is governed by local quadratic terms J 2 versus the linear ones, J, thus implying that the nonlinearity is formally controlled by the Mach number instead of the Reynolds number, which makes a huge difference.
As a result, it has been argued that CLB (Carleman-Lattice Boltzmann) might work better than CNS (Carleman-Navier Stokes) 16 .Indeed, both 16 and 26 report excellent convergence up to Reynolds numbers of the order 10 − 100, with just two or three Carleman iterates.
We tested CLB on a 2D lattice of 32 × 32 points in the grid G, where the initial conditions were set to for each of the gridpoints ⃗ x ∈ G, with a, b = 1, 2 and a ̸ = b.
The values of the distribution functions f i with are obtained accordingly to standard methods.This is shown in Figs.2(a) and 2(b), where we plot the macroscopic density J 1 ( x,t) for a random point x, found by the exact LB method and the corresponding approximate solutions provided by the Carleman linearization at orders 1 and 2. The two figures refer to different choices of the parameters ω, namely ω = 1 in Fig. 2(a) and ω = 1.5 in Fig. 2(b).The other parameters and the choice for the initial conditions are detailed in the caption.In Figure 2(c) we plot the average deviation ⟨RMSE⟩ (Root Mean Squared Error) between the results obtained by the exact LB calculation and the ones obtained by the Carleman linearization of different orders.It is defined as Fig. 2(c) shows the results for ω = 1 and ω = 1.5 respectively.It's worth to notice that although the Mach number should be the relevant parameter for accounting the non-linearity of the system in the LB framework, larger deviations are obtained for higher Reynolds (lower viscosity).
Unfortunately, the free-streaming also brings an undesired non-local coupling of the Carleman variables.This is easily understood by considering for instance the local Carleman pair f i j (⃗ x,⃗ x;t) ≡ f i (⃗ x,t) f j (⃗ x,t).At time t + 1 (∆t = 1), this flies into the nonlocal pair A sketch of this feature is shown in Fig. 3. Due to this nonlocality, the number of CLB variables in a block-time algorithm marching from t = 0 all the way to the end time t, scales like (n v G) k , which is totally unviable despite the low-k convergence.Ultimately this amounts to an exponential depth problem of the corresponding quantum circuit, which is found to lie close to the theoretical upper bound De ∼ 4 q = N 226 .Of course, a single step algorithm from t to t + ∆t would feature (n v ) k scaling.However, this faces with the notorious problem of measurement.
Hence, despite its wonderful mathematical structure and fast convergence, the CLB procedure is blocked off by the depth barrier.In the end, this due to the fact that then CLB matrix is nonlocal hence it cannot be expressed as a compact product of the 2 qubit matrices of the native quantum gates.
In light of the lesson learned with CLB, it is worth revisiting CNS in the hope of mitigating the depth problem.

B. Carleman Navier-Stokes
For the sake of simplicity let us consider the compressible Navier-Stokes equations which, upon spatial discretization and advanced in time by means of a simple Forward-Euler scheme: where the right hand side is evaluated at (⃗ x,t).In the above, σ ab is the dissipative tensor and D a denotes the Centered Finite Difference operator described in section V, which is the  The discretization shown above, based on Forward-Euler time-stepping and Centered Finite Differences, may appear quite naive as it is well known to suffer from various forms of numerical instability; on the one hand, in the high Reynolds regime, such instabilities are usually cured by applying some form of upwinding or by adding an artificial excess diffusion 15,22 , on the other hand, at low Reynolds, issues related to local incompressibility regions in the flow warrant either the use of different representations for the dependent variables 4,24,33 or the use of suitable projection techniques to filter out spurious modes 7 .Nonetheless, the simple discretization approach presented is attractive in the current context due to its locality properties that are discussed below (see Appendix V), and does therefore deserve to be tested and assessed.
The second level of the Carleman approximation requires an equation for the second order tensor ρ −1 J a J b , which we approximate assuming weak compressibility of the fluid.
Yet another problem would arise with a non-ideal equation of state P = P(ρ) but for the present purposes we shall stipulate an linear ideal-gas relation P = ρc 2 s .With this position, we need equations for J ab ≡ J a J b and the cubic term ρJ ab .The Euler forward method in Carleman form delivers the follow-ing set of equations written in matrix form as where the greek indices run from 0 to d and we identified J 0 with ρ.The matrices A, B and C can be derived from Eqs. ( 13) and ( 14 and the other terms are 0. By far, A is the most populated matrix.In the case of a two-dimensional system its explicit form is The tensor B has components and the tensor C [... ]0 has non-null components The evolution of the second order Carleman variables is obtained by multiplying the Eq. ( 15) by J β , which delivers several terms of order two to six.
The tensor product between the vectors J α and J β leads to nonlocal terms with components J αβ (x, y) for any grid point ⃗ x,⃗ y 17 , and with a total dimension (n v G) 2 .However, a much smaller vector can be obtained if we consider only the local variables J αβ (⃗ x) ≡ J αβ (⃗ x,⃗ x) which lead to a second order vector of dimension just n 2 v G. Thus, the potential advantage of CNS over CLB relates to the count of Carleman variables, based not only on the number of variables per site but potentially also on better locality.In fact, the conservative form of Eqs. ( 13),( 14) is preserved also at higher Carleman orders.All the terms can in fact be written under divergence and the derivative operator acts relating local variables at different locations.Nonetheless, the Carleman variables remain function of a single location.For instance, the discrete derivative applied to the first term in the square brackets in Eq. ( 14) in direction ⃗ e a assumes the form: Note that in the LBM formulation, this role is taken by the streaming step, which requires the calculation of the non local variables J ab (⃗ x,⃗ y) at each grid point, see Fig. 3 and analogously for the higher order Carleman terms 26 .This feature is further discussed in the Appendix.
At Carleman level K = 1, we have V 1 ≡ [ρ, J a ] which makes 4 in d = 3 against 27 for LB, meaning two qubits for CNS and five for CLB.With a 4 q depth, this means 4 2 = 16 for CNS and 4 5 = 1024 for CLB, which is already a significant gain, even leaving nonlocality apart.
In order to make a fair comparison between the two methods, we analysed the Kolmogorov-like flow with the same initial conditions as defined in Eq. ( 9).We exploit Eq. ( 8) to derive the corresponding viscosity ν, implicit in the NSE (14) in the stress-tensor term σ ab .We show in Fig. 4 the solution obtained by Euler forward marching method for the same point plotted in Fig. 2, with ν = 1/6 in Fig. 4 From Figs. 4(a) and (b), we see that the convergence of CNS is valid only for small time frames up to O(10) timesteps.Increasing the Carleman order leads to a better approximation of the dynamics, even though the Carleman approximation remains valid for a very small time frame.We observe that the time frame of the CNS convergence is limited to less than O(100) time steps.We are aware of the fact that Euler timemarching with centered finite-differences is unconditionally unstable.However with time step currently used in the simulation (∆t = 0.01), such instability is not expected to become apparent up to t ∼ 1, since it is driven by a negative viscosity of order U 2 ∆t 2 ≃ 5 • 10 −3 .Furthermore, we wish to point out that in the time frame shown in Fig. 4 the density remains constant up to second digit.Finally, we are led to conclude that convergence of CNS is significantly poorer than CLB.This might be due to the Reynolds versus Mach argument.
Hence, contrary to what happens for CLB, CNS has convergence issues, which cannot be traced to the lack of stability of the Euler method.Much more work is needed to handle the behavior shown above.Incidentally, one may also consider resorting to different and possibly more effcient linearization strategies than Carleman (see, e.g, 24 ) IV. CARLEMAN-GRAD PROCEDURE Before closing, we discuss a third alternative which may combine the best of CLB and CNS, namely the application of the Carleman procedure to the Grad formulation of generalized hydrodynamics 10 .The basic idea is to take progressive moments of the Boltzmann probability distribution and inspect the resulting open hierarchy of hyperbolic PDE's.At order three in the Grad expansion, we obtain where ρ = f dv is the fluid density, J a = f v a dv is the fluid current, P ab = f v a v b dv is the momentum flux tensor and Q abc = f v a v b v c dv is the energy flux tensor.
As is well known, this is a hyperbolic superset of the Navier-Stokes equations, which are recovered in the limit of weak departure from local equilibrium.Under such limit, the third equation can be closed, by assuming adiabatic relaxation of P ab to its equilibrium, namely: where τ = 1/ω and the equilibrium expressions, and The appeal of the Grad formulation is its hyperbolic character, which reflects in the conservative nature of the equations.
Euler time marching delivers: (27)  P ab (t + ∆t) = P ab (t) − ∆tD c Q eq abc − ω∆t(P ab − P eq ab ) (28) This scheme preserves locality because P ab is an independent variable under divergence.
The first order Carleman step is to set the nonlinear terms J ab = 0 and J abc ≡ J a J bc /ρ = 0.  13) and ( 14), and comparison with CNS up to the fourth order for ν = 1/6 in (a) and ν = 1/18 in (b).The parameters of the Kolmogorov flow at t = 0 are

A. Picard iterations
However, further Carleman steps require knowledge of J ab to compute P eq ab .The dynamic equation for J ab can be obtained the same way as for CNS, at the cost of introducing six additional fields, for a total of 16, namely still 4 qubits.
A possible alternative is to resort to Picard iteration by computing the above terms using the values of ρ and J a from the previous iteration.In other words, by storing a separate copy of the Carleman variable at each iteration level l, one would compute J p+1 ab = (J p a J p+1 b + J p+1 a J p b )/(ρ p + ρ p+1 ), which is linear in the Carleman variables J p+1 a .With P Picard iterations this brings about another 4P Carleman variables.Assuming fast convergence, say P = 3 takes the count to 10 + 12 = 22, namely 5 qubits per site and a circuit depth 2 10 = 1024, possibly doable in the near future.
Clearly, this iterative procedure could be applied to the CNS framework as well, but the Grad picture offers a number of specific advantages.
First, like in LB, dissipation occurs via adiabatic relaxation of P ab to its equilibrium value, a fully local process which does not entail any communication in space, hence no Laplace operator as required in the NS picture.Second, again like in LB, the nonlinearity appears to be controlled by the Mach number.Third, since P ab obeys its own conservative evolution equation, weak-compressibility effects are naturally incorporated with no need of ad-hoc numerics.
As compared to CLB, the advantage is a lesser number of variables and a non-directional streaming, which favors locality.This said, the iterative procedure adds another layer of complexity, hence its benefits must be carefully weighed against the corresponding costs.Work is current underway to provide a quantitative assessment of the CLB procedure.

B. Summary of Carleman linearization
Summarizing, it appears like Carleman linearization for fluids presents a dual picture: on the one side, CLB is very compact and elegant, but it involves an exponential growth of variable not only because it needs many variables per site but because free-streaming couples variables across sites.
On the other side, the CNS framework offers a drastically reduced number of Carleman variables, first because there are less variables per site and potentially more locality, as long as conservative difference schemes can be developed.
Finally, the CG framework may offer an optimal trade-off between the two.
Much further (hard) work is needed to assess whether any of the three Carleman routes discussed in this paper may finally open the way to the quantum simulation of classical fluids.

V. APPENDIX: LOCAL AND NON-LOCAL FINITE DIFFERENCE TERMS
Using centered differences , with ∆x = ∆y = 1, the divergence of a generic second order tensor T ab reads as follows: These expressions are local since they do not involve any coupling between T ab at different sites, all we need to store is T ab (i, j).

Figure 2 .
Figure 2. The dynamics of the Kolmogorov-like flow on a 32 × 32 grid, plotted for a point of the grid, evolved with the LB method and comparison with CLB at first (CLB 1 ) and second (CLB 2 ) order for ω = 1 in (a) and ω = 1.5 in (b).In (c) the average Root Mean Squared Error ⟨RMSE⟩ between Carleman at second order and exact lattice Boltzmann.The parameters of the Kolmogorov flow at t = 0 are u 1 = u 2 = 0.1, k 1 = k 2 = 1.

2 [
T xx (i + 1, j) − T xx (i − 1, j) + T xy (i, j + 1) −T xy (i, j − 1(i + 1, j) − T yx (i − 1, j) + T yy (i, j + 1)−T yy (i, j − 1)] Figure 1.Convergence of the Carleman Linearization for the logistic equation at increasing orders of the truncation.As one can see, increasing the order of truncation extends the convergence time horizon, beyond which the solution drastically departs from the exact one.
Figure 4.The dynamics of the Kolmogorov-like flow on a 32 × 32 grid, plotted for a point of the grid, evolved with the Euler method from the NS equations (