Quantum Simulator for Transport Phenomena in Fluid Flows.

Transport phenomena still stand as one of the most challenging problems in computational physics. By exploiting the analogies between Dirac and lattice Boltzmann equations, we develop a quantum simulator based on pseudospin-boson quantum systems, which is suitable for encoding fluid dynamics transport phenomena within a lattice kinetic formalism. It is shown that both the streaming and collision processes of lattice Boltzmann dynamics can be implemented with controlled quantum operations, using a heralded quantum protocol to encode non-unitary scattering processes. The proposed simulator is amenable to realization in controlled quantum platforms, such as ion-trap quantum computers or circuit quantum electrodynamics processors.

Transport phenomena in fluid flows play a crucial role for many applications in science and engineering. Indeed, a large variety of natural and industrial processes depend critically on the transport of mass, momentum and energy of chemical species by means of fluid flows across material media of assorted nature1. The numerical simulation of such transport phenomena still presents a major challenge to modern computational fluid dynamics. Among the reasons for this complexity stand out the presence of strong heterogeneities and huge scale separation in the basic mechanisms, namely advection, diffusion and chemical reactions23. In the last two decades, a novel concept for the solution of transport phenomena in fluid flows has emerged in the form of a minimal lattice Boltzmann (LB) kinetic equation. This approach is based on the statistical viewpoint typical of kinetic theory45. LB is currently used across a broad range of problems in fluid dynamics, from fully developed turbulence in complex geometries to micro and nanofluidics67, all the way down to lattice gas automata8 and quark-gluon applications9.

Recent improvements in ion trap and superconducting circuit experiments make these platforms ideal for challenging quantum information and simulation tasks. Trapped-ion experiments have demonstrated quantum information and simulation capabilities101112, including the quantum simulation of highly correlated fermionic systems13, fermionic-bosonic models1415 and lattice gauge theories16. Superconducting circuit setups can host nowadays top-end quantum information protocols, such as quantum teleportation17 and topological phase transitions18. These quantum devices are approaching the complexity required to simulate both classical and quantum nontrivial problems, as proposed by Feynman some decades ago19. Efforts in designing quantum algorithms for the implementation of fluid dynamics make use of quantum computer networks2021. In these works, the quantum degrees of freedom are used on the same ground as classical parameters, and the exponential gain of quantum computers is not properly exploited. In contrast, systems described by pseudospins coupled to bosonic modes, such as the aforementioned ion-trap and superconducting circuit platforms, can enjoy quantum superposition and have advantages with respect to pure-qubit quantum computers in simulating fluids.

In this article, we propose a quantum simulation of lattice Boltzmann dynamics, using coupled pseudospin-boson quantum platforms. Based on previously established analogies between Dirac and LB equations, we define here a full quantum mapping of transport equations in fluid flows. The LB dynamics is simulated sequentially by performing particle streaming and collision steps. The non-unitary collision process can be implemented with an heralded protocol, by sequential collapses of an ancillary qubit. The proposed mapping is amenable to realization in trapped-ion and superconducting circuit platforms.

Results

The lattice Boltzmann equation is a minimally discretized version of the original Boltzmann’s kinetic equation, in which the fluid is modeled as an ensemble of particles that move and collide within a uniform lattice. The lattice Boltzmann dynamics is described by the equation

Here, is the ith component of the particle fluid density associated with the lattice site at the time t, and with discrete velocity . The macroscopic fluid density at the site is retrieved as , while the fluid velocity is defined as the weighted sum of the discrete velocities, . The velocity components , with , satisfy mass-momentum-energy conservation laws and rotational symmetry. Typical lattices are D2Q9 or D3Q15 models, for the case of two dimensions with 9 speeds, and three dimensions with 15 speeds, respectively22.

Collisional properties are here expressed in scattering-relaxation form, making use of the local equilibrium distribution . The LB approach to compute the dynamics associated with Eq. (1) uses sequential computational steps. One initially performs a displacement (free-streaming) of each distribution component towards the nearest-neighbor lattice site pointed at by the discrete velocity . From there, the equilibrium distribution function is computed and the outcome of the collisional process is retrieved. Further iterations of these calculations allow the propagation of the lattice dynamics in time. We address the question of whether all these steps can be performed in a quantum simulator with practical quantum computing protocols.

The formal analogy between the Dirac and LB equations was first highlighted in423, where the velocity distribution of the particle is treated in a similar fashion as a relativistic spinor. This analogy is further exploited in the Majorana representation of the Dirac equation, by using real spinors24. The Dirac (Majorana) equation reads (ħ = 1 here and in the following)

where we have defined the Dirac (Majorana) streaming matrices , mass term β, and the imaginary prefactor i proper of quantum mechanical evolution.

Notice that the streaming matrices of the LB equation are diagonal, while the α, which generate a Clifford algebra, cannot be simultaneously diagonalized. Additionally, the mass matrix β is Hermitian, while standard collision matrices come in real symmetric form in the LB equation. Therefore, a complete codification of the LB scheme in quantum language requires the implementation of diagonal streaming matrices and of purely imaginary symmetric scattering matrices.

The components of the fluid density distribution function can be encoded in a set of quantum states defined on a proper Fock space. For example, in two dimensions, the distribution of the fluid density over the two coordinates can be described by a real quantum wavefunction that encodes the state of two bosonic modes, as depicted in Fig. 1. In the x-quadrature representation, it reads , where f(x1, x2) is a real distribution and the eigenstate of the quadrature of the first (second) bosonic mode. Several quantum distributions can be used by entangling the bosonic state to a multi-level system, such as a set of pseudospins, therefore the state of the complete system is given by , with η being real-valued coefficients. In order to keep a real-valued representation of , to be identified with a fluid density distribution function, one has to act only with purely imaginary interaction matrices.

Figure 1

(a) The distribution of the fluid density on a 2-dimensional lattice can be simulated, for example, via normal motional modes and internal levels of a set of trapped ions (b). (c) Superposition of two motional modes entangled with pseudo spin states can encode velocity distributions in different lattice directions.

The quantum simulation of the Dirac equation was originally proposed25 and afterwards realized in a trapped-ion experiment26. In general, streaming interactions involving matrices in the Dirac or Majorana representation can be implemented by using a pair of pseudospins coupled to one or more bosonic modes. In terms of creation and annihilation operators for the bosonic mode in the b direction, one can then consider and write Eq. (2) on the pseudospin-bosonic Hilbert space of ,

where K stands for the pseudospin-boson coupling and α act upon the pseudospin degrees of freedom.

Thus, the three streaming matrices are written in the Dirac representation as , in a pseudospin representation and the diagonal mass term as . These streaming matrices are diagonalized via the operators 24, which have to be physically implemented as quantum gates. Defining , the associated generators read , with and . In this way, a purely imaginary streaming step iβ∇ can be built, which mimics the diagonal streaming of the LB equation. The total wavefunction after the streaming steps can be retrieved with a sequential implementation, following the operator splitting method23. For example, in a 2-dimensional lattice, one has

The last collision step C, which scrambles particle distributions in different directions, is discussed below.

Standard collision operators in LB theory are represented by real symmetric matrices associated with non-unitary evolution operators. On the other hand, typical controlled quantum mechanics experiments produce unitary dynamics. Nevertheless, one can probabilistically encode non-unitary dynamics in a quantum device with a heralded protocol, by performing controlled operations conditioned on the state of an ancillary qubit, and then using the state of the latter as a flag for the success of the protocol. We consider a purely imaginary symmetric scattering matrix Ω, whose quantum evolution equation reads , providing a non-unitary evolution operator that describes lattice collisions .

The collision operator can be decomposed in a weighted sum of two commuting unitary operators, , with the constraint , assuming without loss of generality that γ > 0.

Given a specific diagonalizable collision operator C and weight γ, one can then find its decomposition in terms of unitaries. In order to find a decomposition in terms of unitaries, C must first be diagonalized as . This reduces the problem of finding U and U down to an eigenvalue equation, δ = α + γβ, with δ, α and β being the ith eigenvalues of the collision and unitary operators respectively. Notice that, due to the properties of the scattering matrix, . Taking into account the normalization conditions, one has the system of equations

The eigenvalues α, β can now be written as a function of the initial collision operator and weight γ,

The unitary operators U are reconstructed via . The real domain of Eqs (6) provides the range of validity of the method developed here. Simple algebra leads to the set of inequalities

By defining δ and δ as the maximal and minimal eigenvalues of the spectrum of C, the system of inequalities in Eq. (7) can be reduced to one of the two inequalities or , respectively when or . If longer evolution times t are considered, the spectral range of C changes accordingly. The weighted γ-sum derived here can be implemented with quantum computing algorithms, using ancillary qubits and controlled U and U gates27. By measuring the ancilla state, one can determine whether the desired operation has been performed or not. The success of the protocol depends on the weighted sum of unitary operators, with a failure probability .

As P is an increasing function of γ, choosing maximizes the probability of success. This directly connects the simulation time of the scattering process C with the best choice for γ. To propagate the dynamics of a given collision process C, one can split the step time Δt into N time intervals Δt/N and perform the heralded protocol at each step, such that . At each step, one has a collision operator , with an optimal γ0. In this way, as the step size gets smaller, the success probabilities for each step increase, while the total success probability accumulates single success rates from the individual steps. In Fig. 2a, we plot the success probability of the simulation of the single step, as a function of N, for random symmetric purely imaginary matrices. As expected, the success probability per step increases as the size for the single time step gets smaller. The success of the whole protocol is constant and does not depend on N. In Fig. 2b is shown that the optimal protocol is performed at γ = γ0.

Figure 2

Probabilities of success P for the implementation of the scattering process.

(a) Probability of success P per time step of simulating real symmetric random matrices as a function of the number of ancillary measurements N (solid lines), together with accumulated probabilities for the whole protocol (dashed lines). Each curve represents a different instance of a random matrix. (b) Probability of success of a single step as a function of γ/γ0, when N = 10.

Discussion

The scheme proposed can be adapted to a variety of transport fluid problems. As an example, we consider the implementation of an advection-diffusion process in two spatial dimensions. The dynamics of the transported species, e.g. pollutants or bacteria, is described by the equation

where is the scalar field transported by a fluid with space-dependent velocity and constant diffusivity D.

The problem in Eq. (8) can be recast in LB form, as in Eq. (1). The corresponding equilibrium distribution function is defined as

with w = 1/4, . Note that, by definition, the space-time dependence of the local equilibria is entirely carried by the macroscopic fields ρ and .

The scattering matrix reads , where , Ai(2) = cix ≡ , and are the four eigen-vectors.

The first three corresponding eigenvalues are given by

which follows from mass conservation and the expression of the diffusion constant , respectively. By choosing different values for ω2 and ω3, one can implement anisotropic diffusivities along the x and y directions. Finally, ω4 dictates the decay of higher order fields of no direct macroscopic significance, hence it is chosen as ω4 = 1 so as to annihilate the corresponding field in a single collision step.

The relative scattering matrix Ω is defined by . By choosing a Couette flow, e.g. , where L is the typical size of the fluid domain, one has , with and . Here, velocities are numbered 1 ÷ 4 counterclockwise starting from the +x direction.

The latter defines the quantum scattering matrix as composed of three contributions, namely , where is proportional to the position quadrature of the bosonic mode associated with the y direction. The three contributions to the scattering matrix represent classical linear wave propagation and damping, mass conservation and macroscopic advection, respectively. They can be implemented with the quantum simulation protocol previously introduced. The bounds to γ can be obtained, e.g., for the first contribution to the scattering matrix −A, by computing the spectrum of for different time steps Δt, for D = 0.05, in units of 1/ω4. The result is shown in Fig. 3.

Figure 3

Spectrum of a collision operator (solid red line) for advection-diffusion process of a four-speed lattice as a function of the evolution time step Δt, in units of 1/ω4.

The allowed region for γ is bounded by dashed blue lines using Eq. (7) and shadowed in the picture.

Natural quantum platforms for prospective implementation of the proposed scheme could be ions trapped in linear Paul traps or superconducting circuit setups, in which the sequential streaming and collision steps in Eq. (4) can be realized. The pseudospin-bosonic state can be encoded, in the case of ion traps, in the internal level and motion modes of the ions28, while in a superconducting architectures, one can use the first levels of charge-like qubits, e.g. transmon qubits, and microwave resonators29. One may consider opening similar avenues in other quantum technologies as is the case of quantum photonics30 and Bose-Einstein condensates31.

A practical implementation of the protocol proposed can make use of many-body interactions, involving couplings with bosonic modes. These type of gates have been considered in superconducting architectures32 or in ion-trap platforms33. For a four-speed lattice, the diagonal streaming processes can be realized with a combination of a qubit-boson interaction and two entangling gates among the qubits. For example, the corresponding evolution operator for the streaming in the X direction can be written as

where we have defined an entangling operation between the two qubits and local rotations of the i-th qubit about the j-th axis, . The U interaction can then be diagonalized in the qubit space via the realization of two S matrix, , which can be achieved by a combination of entangling two-qubit gates and a phase gate. In the case of more internal degrees of freedom of the lattice, the two-body entangling gate can be substituted by a collective interaction . Similar reasoning applies to the streaming in the Y direction, considering a different bosonic mode a. The unitary matrices that implement the collision process can be implemented in a controlled way27 by using an additional ancillary qubit Ψ and performing the quantum gate . These gates can be decomposed in general with a Lie-Trotter-Suzuki decomposition, in terms of many-body interactions of the type , , with . These many-body interactions can be obtained by sequential implementation of collective gates and single qubit rotations3233. In the case of a Couette flow, the term with the linear spatial dependence of the scattering matrix can be implemented by considering a single qubit rotation entangled with a bosonic displacement, similar to Eq. (11). The quantum resources necessary to implement Lie-Trotter-Suzuki decompositions scale polynomially in the internal degrees of freedom of the lattice, and sub-polynomially in the digital error34. Notice also that the quantum resources needed are invariant with respect to the size of the simulated lattice. The latter will depend on the accessibility and readability of distributions over Fock spaces in practical implementations, e.g. the ability to characterize distributions over current quadrature in superconducting architectures3536.

Note that the above scheme readily extends to the case of reactive flow, by augmenting the collision operator with a local source term proportional to the chemical reaction rate. Such kind of advection-diffusion-reaction phenomena in complex geometries, say catalytic reactors, represent a very active area of applications of the LB scheme. Further developments may include the implementation of hydrodynamic non-linearities to model the Navier-Stokes fluid dynamic equations. This requires the inclusion of quadratic terms in the LB equilibrium distribution. Such nonlinear behavior can be provided in a quantum mechanical experiment by preparing multiple copies of the system37, feedback mechanisms38, or non-unitary operations induced by measurements.

We have developed a protocol to reproduce the dynamics of fluid transport phenomena in a quantum mechanical experiment, by using pseudospins coupled to bosonic modes that can be implemented in different quantum platforms. This proposal paves the way to quantum simulation and retrieval of complex classical fluid dynamics in controlled quantum systems.

Additional Information

How to cite this article: Mezzacapo, A. et al. Quantum Simulator for Transport Phenomena in Fluid Flows. Sci. Rep. 5, 13153; doi: 10.1038/srep13153 (2015).