The simulation of large nonlinear dynamical systems, including systems generated by discretization of hyperbolic partial differential equations, can be computationally demanding. Such systems are important in both fluid and kinetic computational plasma physics. This motivates exploring whether a future error-corrected quantum computer could perform these simulations more efficiently than any classical computer. We describe a method for mapping any finite nonlinear dynamical system to an infinite linear dynamical system (embedding) and detail three specific cases of this method that correspond to previously studied mappings. Then we explore an approach for approximating the resulting infinite linear system with finite linear systems (truncation). Using a number of qubits only logarithmic in the number of variables of the nonlinear system, a quantum computer could simulate truncated systems to approximate output quantities if the nonlinearity is sufficiently weak. Other aspects of the computational efficiency of the three detailed embedding strategies are also discussed.

Plasma applications have historically been at the forefront of high-performance computing, and plasma physics computations are routinely performed using the largest available supercomputers.1,2 Even so, resolution, domain size, and dimensionality often limit the realism of plasma simulations. More computing power would allow for better understanding and more accurate predictions of the behavior of both natural and experimental plasmas. Quantum computers may eventually provide some of this computing power. Although they are at an early stage of development, error-corrected quantum computers have the potential to be much more powerful than classical computers, at least for some tasks.3 Plasma simulations are typically nonlinear, and the investigation of the power of quantum computation for nonlinear simulations, including plasma simulations, is an active research area.4–10 

Consider a kinetic plasma physics computation composed of the following steps. First, particle distribution functions fs(x,v,t), where s is the particle species, are initialized at time t =0 to specified functions of position x and velocity v. Second, the distribution functions are time-evolved numerically using the Vlasov equation or similar. Finally, some simple functional of fs(x,v,t) is extracted as the output of the computation. This process may also be repeated a small number of times to obtain different output quantities or to evaluate outputs at different t. If all these steps can be done efficiently, many important kinetic plasma physics problems can be solved. But in practice, the number of variables required to represent fs(x,v,t) makes this computation extremely expensive in general.

An analogous situation occurs when simulating many-particle, quantum-mechanical systems. Both the initial quantum-mechanical state and the desired output may have simple functional forms, but, since the number of variables needed to represent the state grows exponentially with the number of particles, this computation can be intractable. On the other hand, with a quantum computer, exponentially large states can be represented efficiently, and general quantum-mechanical simulations can be performed with exponentially lower computational costs than with a regular computer. This possibility is what led Feynman to originally propose the idea of quantum computation.11 

Moreover, universal quantum computers can obtain speedups for computations that are not quantum-mechanical in origin. For example, when the Vlasov–Maxwell system is linearized about a Maxwellian background, it can be transformed to have the form of a quantum-mechanical system.12 Then a quantum algorithm, i.e., an algorithm designed to run on a quantum computer, can potentially simulate this plasma physics system with a large speedup. In particular, we previously formulated an efficient quantum algorithm for the linear Landau damping problem and discussed extensions to six-dimensional phase space and full electromagnetics.12 Others who have studied the application of quantum algorithms to classical plasma physics problems include Parker and Joseph13 and Dodin and Startsev.5 

Quantum algorithms that perform general linear computations with a large speedup have also been designed. An important example is provided by the linear problem Ax=b, where A is an invertible matrix and b is a known vector. This can be solved by a quantum linear systems algorithm (QLSA) with costs only logarithmic in N, where N is the length of b.14,15 To go into a bit more detail, it is assumed that one can efficiently prepare a quantum state,

where the proportionality constant is determined by requiring |ψi to be normalized. The index j runs over the indices of b, and the state |j corresponds to a sequence of qubits with values matching the bits of the binary representation of j. Whether the |ψi state can be prepared efficiently depends on b, but it is possible in some general cases, such as if each bj can be computed efficiently and maxj|bj|/|b|=O(1/N), where |b|=jbj*bj denotes the complex vector magnitude.16 The QLSA then outputs a state,

where x=A1b. It is also assumed that the desired output quantity can be efficiently extracted from |ψf. A QLSA can then obtain an exponential speedup over a classical algorithm solving Ax=b subject to conditions on the matrix A, including that its condition number is not large.14 

Quantum linear systems algorithms can be used to perform various other linear computations, including the evolution of systems of linear differential equations.17 However, we are really interested in the evolution of nonlinear systems. A good example is the Vlasov–Poisson system for electrons,

where the t dependence of each variable has been omitted for brevity. A more compact expression for this system can be obtained by replacing the electric field E(x) with its value obtained through Coulomb's law,

(1)

where we have also rescaled quantities so that they are all dimensionless. The Debye length λDe is the distance unit, the plasma frequency ωpe is the inverse time unit, and ne/(λDeωpe)3 is the unit in which f(x,v) is expressed, where ne is the electron number density. To solve Eq. (1), f(x,v) can be represented on a grid, but a very large number of grid cells will be required in general.

One way to simulate nonlinear dynamics on a quantum computer is with the Koopman–von Neumann approach introduced by Joseph4 and studied further by Dodin and Startsev.5 We explore a different approach, with a goal of having the quantum computational costs scale only logarithmically with the number of variables N for simulation problems that have classical costs linear in N. For instance, this would mean costs logarithmic in the number of phase-space grid cells used to represent Eq. (1). But to apply quantum algorithms to a nonlinear simulation problem, we first map the nonlinear system to a linear one. The basic idea here can be understood as follows. Consider the extremely simple nonlinear system,

If we introduce a set of variables yr=xr for r=0,1,2,, then we obtain the infinite linear system

This demonstrates one way that a nonlinear system can be mapped to an infinite linear one. In Sec. II, we consider the general problem of mapping nonlinear dynamical systems to infinite linear ones and review three specific mappings, with the one in Sec. II A generalizing the above example. Then in Sec. III, we suggest a way to truncate the infinite-dimensional linear systems to obtain linear systems of size poly(N), where N is the number of variables of the original nonlinear system, although these truncated systems might only provide good approximations to output quantities if the nonlinearity is sufficiently weak. Some analysis of the efficiency of quantum computation based on the truncated linear systems is provided in Sec. IV, and we discuss our findings in Sec. V.

First, we describe a method for expressing nonlinear dynamical systems as infinite-dimensional, linear dynamical systems. We write the original dynamical system as

(2)

where z(t) is the vector of variables and F(z) is a vector function of the components of z. We use N to denote the length of these vectors, i.e., the number of variables. Equation (2) is a system of ordinary differential equations; partial differential equations, such as the Vlasov equation, can be converted to this form through, e.g., spectral, finite element, or finite volume methods.

To map Eq. (2) to a linear system, we define a set of states,

(3)

where ŵ is a vector of operators and |0 is a fixed state. Up to an overall normalization factor, these states just generalize those considered by Kowalski.18 Additionally, we introduce a vector of operators ẑ that satisfies

(4)

The reason for denoting these operators by ẑ is revealed by the following evaluation:

(5)

and therefore

(6)

Note that, since the ẑ operators have eigenvalues of every variable value taken at any time, continuous evolution of z(t) implies that the space spanned by the |ψ(t) states is infinite-dimensional. Also, Eq. (6) implies that the ẑ operators commute within this space: ẑjẑk|ψ(t)=ẑkẑj|ψ(t) for every |ψ(t) state and linear combinations of them. One consequence is that the embedding does not introduce any quantum-mechanical effects, such as the uncertainty principle, into the evolution of the classical system.

Next, differentiating Eq. (3) with respect to time gives

(7)

where we assume that all components of F(z) are analytic functions. Then F(ẑ) is defined by replacement of z with ẑ in the power series representation of F(z). Equation (7) shows that the |ψ(t) states evolve linearly according to

(8)

where

(9)

The finite-time evolution can be expressed as

(10)

where |ψ(0)=ez(0)·ŵ|0 depends on the initial state of the classical system, z(0).

Viewing the linear evolution as occurring in a Hilbert space containing the |ψ(t) states, the inner product with some state |c yields an output quantity,

(11)

that depends on the final state of the nonlinear dynamical system. To go into more detail, we now consider a few specific instances of the embedding method. All of these have been previously studied in some form by other authors (e.g., Refs. 19–26) but we give an exposition of them here to illustrate various aspects of linear embedding and to collect results that will be used in Secs. III–V.

Let |n, where n is any length-N tuple of non-negative integers, be an orthonormal basis for our Hilbert space. Now suppose that the forms of ẑ and ŵ are

(12)

where

(13)

It is straightforward to check that these operators satisfy both parts of Eq. (4) with |0 being the n=0 basis state. Additionally, the ŵ operators satisfy

(14)

which can be used to simplify some expressions.

To concisely express the |ψ(t) states, we first write |n as

(15)

For instance, with N =2, we would have |n=|n1|n2. Now

(16)

Additional insight can be obtained after defining a number operator n̂ by

(17)

Using n to denote the n̂ eigenvalues, states can be broken up into components of different n. For example, the n =1 component of |ψ(t) is

(18)

and the n =2 component is

(19)

Another observation is that ẑ decreases n by one, except the n =0 component, |0, which it annihilates; and ŵ increases n by one. Those familiar with the technique of Carleman linearization22,27 should see that this particular linear embedding is just a representation of Carleman linearization in a Hilbert space. The development of an efficient quantum algorithm for dissipative nonlinear dynamical systems using Carleman linearization has been recently reported in Ref. 8.

Finally, we consider output quantities. Any variable can be obtained as an inner product,

(20)

Similarly, a linear combination of variables

(21)

is equivalent to c|ψ(t) with

(22)

Polynomials of the variables can also be obtained by adding n >1 components to |c. For instance, with |c=|ej+ek,c|ψ(t)=zj(t)zk(t). More generally, for any polynomial of the variables, there is a specific state |c for which c|ψ(t) evaluates that polynomial.

Now suppose that ẑ=â and ŵ=â, where â are standard bosonic lowering operators. This version of linear embedding has been extensively explored by Kowalski.28 In this case, Eq. (4) becomes

(23)

which amounts to the standard commutation relations and the statement that |0 is the ground state. Additionally, we can use |n to denote the occupation basis and express the number operator as

(24)

Now the eigenvalues n of n̂ represent numbers of bosonic particles.

In this version of linear embedding, the evolution operator is

(25)

and the states are

(26)

These are just coherent states; Eq. (6) becomes

(27)

As in Sec. II A, output quantities can be obtained with inner products. The only difference is that, for quadratic and higher-degree outputs, the 1/nj! factor in Eq. (26) can have an effect. For example, |c=|2ej yields c|ψ(t)=zj2(t)/2 rather than zj2(t).

While the states |ψ(t) evolve linearly, this linear evolution is generally nonunitary with both the Carleman and coherent states embeddings. One simple way to see this is to note that |ψ(t) [Eqs. (16) or (26)] changes in normalization when |zj| for one j changes. Nonunitary evolution presents a difficulty for efficient quantum computation, so a linear embedding that gives unitary evolution may be preferable.

This time, take ẑ=x̂ and ŵ=ip̂, where x̂ and p̂ are dimensionless versions of canonical position and momentum operators, respectively. Then

(28)

and both parts of Eq. (4) are met with |0=|x=0, the x=0 position eigenstate. This version of linear embedding was introduced by Koopman and von Neumann,19–21 and it has been studied in various forms by different authors.24–26,29 Since the eigenvalues of x̂ are real, the variables must be real, and to signify this, we switch to using x to denote the variables instead of z. To maintain the reality of x,F(·) must be a real function as well. The |ψ(t) states are position eigenstates, and they can be expressed as

(29)

Since the translation operator eix(t)·p̂ is unitary, the |ψ(t) states have the same normalization as |x=0.

The evolution operator is

(30)

To study this embedding in the occupation basis |n, we introduce

(31)

It is easily checked that these â obey

(32)

and thus they are standard bosonic lowering operators. Using the expressions for â and â in the occupation basis, the evolution operator

(33)

can also be expressed in an occupation basis.

To express |ψ(t) [Eq. (29)] in the occupation basis, we need to relate the position and occupation bases. However, since Eq. (31) is equivalent to the relation used in the quantum harmonic oscillator (QHO) problem with =mω=1, the required result is well known: the one-dimensional QHO eigenstates in the position basis are

(34)

where Hn(x) are the physicists' Hermite polynomials. The form of |ψ(t)|x(t) follows directly:

(35)

where we have chosen to drop the unimportant factors of π1/4.

Output quantities work somewhat differently in this version of linear embedding. The components of |ψ(t) in the nb subspace, where n denotes the eigenvalues of the total number operator n̂, are polynomials of the variables of degree b times a factor of ex·x/2. Therefore, if we write

(36)

for |c a state in the nb subspace, then p(x) is a polynomial of the variables with degree b. For example, |c=|ej yields p[x(t)]=2xj(t). A potentially simpler way to express an output quantity is as the expectation value of an observable: g:=ψ(t)|g(x̂)|ψ(t) for some analytic function g(x). However, only c|ψ(t) can be evaluated within the nb subspace, which is crucial to the analysis in Sec. III. The result [Eq. (60)] indicating that a truncated system can approximate the output for a sufficiently weak nonlinearity applies to c|ψ(t) values, not any g value. For this reason, we exclusively study the c|ψ(t) way of expressing outputs.

The ex·x/2 factor in Eq. (36) can complicate the representation of desired outputs, which are frequently polynomials of the variables. However, this issue is avoided if x·F(x)=0 for all x. This ensures that dt(x·x)=0, so the ex·x/2 factor in Eq. (36) is constant. Additionally, this has a consequence for the form of the evolution operator: for x·F(x)=0 to hold for all x, it must vanish identically, meaning that all terms are canceled. Therefore,

(37)

holds for r̂ any vector of commuting operators. Now, M̂ [Eq. (33)] can be expanded to give a sum of terms, each being proportional to a product of raising and lowering operators. This includes the terms that occur in the expansion of

(38)

which are made purely of raising or lowering operators. However, Eq. (37) can be applied to find that Eq. (38) evaluates to zero. All the remaining terms in the expansion of M̂ [Eq. (33)] have at least one raising and lowering operator. Also, supposing that F(x) is a polynomial in x of degree g, then each term has a maximum of g +1 operators. Therefore, M̂ couples between occupation basis components for which n, the n̂ eigenvalue, differs by at most g − 1.

The assumption that x·F(x)=0 for all x is nontrivial, but it holds in some important scenarios. For instance, it holds when

(39)

for any W(x), with the indices handled cyclically: if j is the last, then j +1 is the first. Equation (39) is a relevant form since it can be obtained when a first-order derivative in a partial differential equation is represented using centered differences and periodic boundary conditions.

Moreover, Eq. (39) can be generalized to a large degree. First, instead of the index j running over the usual vector components, it can run over any subset of the components, in any order, so long as they are handled cyclically. Components of F(x) that are not in the subset are taken to be zero. Second, F(x) can be any linear combination of terms of the described form, and each term can have a different W(x). As shown in  Appendix A, the electrostatic Vlasov equation [Eq. (1)] is an example of a partial differential equation that can be discretized to have this general form, implying that x·F(x)=0, where the variables x are the values of the distribution function on a grid.

Hereafter we assume that, when the position-space embedding is applied, the system satisfies the x·F(x)=0 condition. Therefore, the factor of ex·x/2 is a constant, and it is convenient to divide this out of the state. Then |ψ(t) becomes

(40)

and the output quantities [Eq. (36)] are updated to

(41)

where, as before, p(x) is a polynomial of the variables with a form determined by the state |c.

Evolution in the position-space embedding can also be described in terms of a “Hamiltonian” Ĥ=iM̂, and Ĥ can be written as

(42)

It is well known and easily verified that, for operators x̂ and p̂ satisfying [x̂,p̂]=i,

(43)

for f(x̂) any analytic function of x̂. Applying this, we find that

(44)

where divF(x̂) is defined by replacing x with x̂ in divF(x):=jFj(x)/xj. For many systems, divF(x) evaluates to zero. If it does not, it is still possible to eliminate that term by extending the original system. In particular, consider two systems, each identical to the original, including in their initial conditions. Whenever Fj(x) contains xj, we replace that xj occurrence with the corresponding variable from the other system. This has no impact on the dynamics, yet it results in the elimination of divF(x) for the doubled system. We prefer this strategy over the elimination of that term through a redefinition of the states (e.g., as done by Kowalski26) since the latter method introduces a factor of

(45)

which must be applied to obtain an output of the c|ψ(t) form, yet Eq. (45) cannot be evaluated in general without the full solution x(t) to the system. We assume that the described extension is applied when necessary to get a system for which divF(x) always vanishes. Then,

(46)

In the occupation basis, Eq. (46) becomes

(47)

Note that applying the standard rules for Hermitian conjugation to Ĥ now yields Ĥ=Ĥ. As a consequence, the embedded evolution will be unitary when restricted to any finite-dimensional subspace of the occupation basis.

Although a finite N will ultimately be required to perform computations, it is straightforward to formulate linear embedding for partial differential equations as well, e.g., as done by Kowalski.28 Consider a dynamical system of the form,

(48)

where F can contain derivatives and integrals of the function f(t) with respect to the coordinates q. If we introduce operators f̂q and ĥq and a state |0 such that

(49)

then the states

(50)

are found to satisfy

(51)

and

(52)

where

(53)

Any of the specific versions of linear embedding can then be applied to partial differential equations. For example, coherent states embedding applied to the electrostatic Vlasov equation for electrons [Eq. (1)] gives

(54)

and

(55)

Quantum computers can perform some linear computations such as matrix inversion with costs only logarithmic in the system size,14 but that does not allow for handling a system of infinite size. Therefore, we seek to approximate a desired output quantity using linear evolution within some finite-dimensional subspace of the linear embedding space. Note that this is necessary even with a finite number of variables N. The infinite dimensionality of the linear embedding space is tied to the variables being represented with infinite precision, which is required to exactly represent continuous evolution in time.

The way in which we choose to truncate the linear embedding space makes a significant difference. For example, one can imagine a truncation of the position-space embedding (Sec. II C) that replaces the infinite-precision variables x with finite-precision variables and represents p̂ with a finite difference matrix in the position basis. In this case, the dimensionality of the truncated space scales as (1/ϵ)N, where ϵ is the chosen precision. Therefore, the space complexity, i.e., the required number of qubits, scales as Nlog2(1/ϵ). However, we will consider an alternative form of truncation and show that its space complexity can be only logarithmic in N. This potentially allows for the approximation of outputs with quantum computational complexity scaling only logarithmically with N, although accuracy is only ensured if the nonlinearity is sufficiently weak [Eq. (60)].

In what follows we assume that F(·) is a polynomial of degree g in the variables for some small integer g and that F(0)=0, i.e., F(·) does not have a constant term. Also, we assume that the desired computational output is a polynomial of degree b in the variables at the final time for some small integer b. If the position-space embedding is used, then we add the requirement that x·F(x)=0 for all x, as discussed in Sec. II C. With the stated assumptions, a couple of properties are shared by Carleman embedding, coherent states embedding, and position-space embedding. First, the output can be expressed as

(56)

for some state |c belonging to the nb subspace of the occupation basis. As usual, n is used to denote the eigenvalues of the total number operator n̂. Second, the evolution operator M̂ couples between occupation basis components for which n differs by at most g − 1.

Next, we rewrite the original dynamical system [Eq. (2)] as

(57)

where A is a matrix, η is a constant, and G(z) is purely nonlinear in z. In particular, G(z) is a degree g polynomial of the variables. Now the evolution operator can be decomposed as

(58)

By itself, M̂0 generates the linearized evolution of the original dynamical system, i.e., the evolution in the limit of η0. Since g =1 for linear evolution, the terms in M̂0 do not couple between different n. The need to consider n>b components comes from ηM̂1, which contains terms that change n by up to g − 1.

Now, we truncate the space such that only the

(59)

subspace is retained for some integer s0. Let c̃(t) denote the approximation to c(t) that is obtained by performing the embedded evolution in the truncated space and evaluating c|ψ(t) after some fixed, finite time t. Then,

(60)

We derive Eq. (60) in  Appendix B. Treating η as a dimensionless parameter characterizing the strength of the nonlinearity, we expect c(t)c̃(t) for small s if η is sufficiently small. Note that how small η needs to be for the approximation to be accurate can depend on t. If the dynamical system and output polynomial are actually linear, then the s =0 truncation is exact, and m =1. More generally, if the nonlinearity is sufficiently weak, which depends on the specific computation, then the output can be approximated with s, and thus also m, being of order one.

To find the dimensionality of the nm subspace, consider N +1 bins into which m particles are placed. The first N bins correspond to the numbers nj of Eq. (15) while the last bin holds any extra particles. The subspace dimensionality is the number of unique placements of the particles into the bins, which is given by the binomial coefficient,

(61)

Noting that

(62)

the number of qubits required to represent the nm subspace is

(63)

Therefore, the space complexity of representing the Eq. (59) subspace is

(64)

Some analysis of the efficiency of a quantum algorithm based on linear embedding is possible without specifying all the details. Let M̂ and |ψ(0) denote M̂ and |ψ(0), respectively, after restriction to the truncated subspace. Then the computation can be expressed as

(65)

where V̂ψ and V̂c are unitary operations that prepare states proportional to |ψ(0) and |c, respectively, from the computational starting state. After application of the Eq. (65) operation, the component along the computational starting state will be proportional to the output quantity, the details of which are determined by the chosen |c state. The technique of amplitude estimation can then be applied to estimate this component using O(1/ε) iterations of the Eq. (65) operation, where ε is an absolute accuracy for the output value.30 In particular, that will yield, up to a normalization factor, an estimate of |c̃(t)|. However, the full complex value can also be estimated using a simple algorithm extension.12 

The |ψ(0) state is a sum of components with particle counts up to m. To represent this using qubits, we can use m registers, each with log2(N+1) qubits. Here N +1 appears instead of N so that there is one extra basis state which can be used to indicate the absence of a particle: for the n < m components there are mn unused registers which should be set to the extra basis state. This representation achieves the Eq. (64) scaling. The |c state is a sum of components (some of which may vanish) with particle counts up to bm [Eq. (59)]. Therefore, |c can be represented in the same way as |ψ(0).

Overall efficiency requires efficient preparation of states proportional to |ψ(0) and |c. In particular, we want these operations to cost poly(logN). We are less concerned with the scaling with m, since we do not expect to achieve overall efficiency when m is large. One basic state preparation strategy is to start by making a superposition over an ancilla register of log2(m+1) qubits with the values of this register representing particle counts. Then, for each particle count nm, the corresponding component of |ψ(0) [e.g., Eq. (19) with t =0] or |c [e.g., Eq. (22)] is prepared. Whether that can be done with poly(logN) complexity will depend on the details of the initial state and output quantity, but it is possible in many cases, including when the maximum |bj|2 is a factor of O(1) larger than the average |bj|2, where bj is the prepared amplitude for state component j, assuming that each bj can be computed efficiently using the index j.16,31

In a quantum algorithm, prepared states must be normalized. Consequently, the output obtained from Eq. (65) will be based on the normalized versions of |ψ(0) and |c. A factor χζ, where χ:=ψ(0)|ψ(0) and ζ:=c|c, must then be applied to the output to obtain the actual result. Therefore, to obtain the result to within some error tolerance δ, the absolute accuracy ε of the output value must be

(66)

and to extract the result to within error δ using amplitude estimation, O(χζ/δ) iterations of the Eq. (65) operation are required.30 

With position-space embedding the truncated system evolution is unitary, so it can potentially be performed efficiently using a Hamiltonian simulation algorithm such as the algorithm by Low and Chuang.32 However, analysis of initial state normalization reveals a complication. The even-degree Hermite polynomials do not vanish at the origin, and therefore the initial state normalization receives large contributions from low-degree components. For example, some of the Eq. (40) components are

(67)

and those contribute

(68)

to ψ|ψ. Equation (68) is minimized by xj=1/4 for all j, and the value at the minimum is N/4. Consequently, an m =2 truncation has χ2>N/4. More generally, the scaling of χ2 with N is 1/χ2=O[Nm/2]. This is a problem for the overall efficiency of an algorithm based on position-space embedding since performing amplitude estimation results in a cost factor of O(χζ/δ).

With coherent states embedding, the initial state has a normalization of

(69)

Truncation of the space will mean that only a finite number of the components are kept, resulting in a lower initial state normalization. Specifically,

(70)

To prevent χ from growing as some power of N, we need |z(0)| to not scale as a power of N. We can achieve that through a rescaling of the variables. In particular, we can switch from z to

(71)

with a constant

(72)

so that γ|z(0)|=O(1). The transformed system is

(73)

and the truncated initial state normalization becomes

(74)

With Carleman embedding, a similar situation occurs. In this case, the initial state normalization before truncation is

(75)

when |zj|<1 for all j and infinite otherwise, although the truncated state normalization will always be finite. Again, we can rescale the variables [Eq. (71)], this time with

(76)

where the factor of 2 is somewhat arbitrary. Next, we can apply

(77)

which holds for non-negative α and β with α+β<1, to bound the initial state normalization as

(78)

Thus, χ=O(1) is achieved.

Applying rescaling [Eq. (71)] to Carleman embedding or coherent states embedding results in a factor of γn being applied to components of |ψ(0) for each particle count nm. At the same time, to maintain the original output quantity, which is assumed to be a degree b polynomial of the variables, the n-particle component of |c is rescaled by γn for each nb. The general form of the output quantity is

(79)

where Cn is a rank-n tensor. Suppose that, prior to any rescaling of the variables, the variable values are bounded as N, while the entries of Cn scale as O(Nn). For example, the variables can represent a bounded distribution function over a grid of size N, and the output quantity can approximate some integral over the distribution function by summing over the values on the grid. Then |z(0)|=O(N), and a rescaling with γ1/N will result in χ=O(1) and Cj1,,jnn=O(Nn/2). Now, the components of |c in the occupation basis are proportional to the entries of the C tensors (in a manner independent of N). Therefore, the scaling of ζ2 with N can be obtained by summing all squared absolute entries of Cn for each n, which yields ζ=O(1) with respect to N in this case.

So Carleman embedding and coherent states embedding can potentially avoid having the cost factor χζ grow as a power of N, but with these embeddings the linear evolution given by eM̂t is generally not unitary. We can perform this nonunitary evolution using a QLSA, e.g., as done by Berry et al.17 However, the more nonunitary the evolution is, the higher the costs will be. Of particular concern is that the condition number κ of the evolution may grow exponentially with the simulation time t, forcing the costs to also grow exponentially with t.

But there is reason to suspect that, in some cases, an exp(t) cost scaling can be avoided. With coherent states embedding, if the system is given by Eq. (57) with the matrix A being anti-Hermitian, then M̂0 [Eq. (58)] is formally anti-Hermitian,

(80)

Additionally, if we order the basis states based on their particle count n, from low to high, then M̂ takes on a block-upper-triangular structure, with one block for each n. This occurs because M̂1 always decreases n, while M̂0 leaves n constant. All the blocks along the diagonal are anti-Hermitian since they are components of M̂0. Then a block-diagonal, unitary transformation can be applied to diagonalize M̂0 while maintaining the overall block-upper-triangular structure. That will result in an upper-triangular matrix with purely imaginary entries along the diagonal. The same holds for the matrix M̂ obtained by truncating M̂ to the nm subspace, which implies that the eigenvalues of M̂ are purely imaginary. As a consequence, the condition number κ of the evolution eM̂t will remain bounded as t if M̂ is diagonalizable. If M̂ happens to be nondiagonalizable, the block structure of M̂ limits the size of any Jordan chain to m, which leads to a bound of κ=O(t2(m1)). A worse scaling with t could still occur if the truncation degree m needed to obtain accurate results increases with t, which will depend on the details of the classical dynamical system.

Another important question is how the costs of implementing the eM̂t evolution scale with N. With the same assumptions as in Sec. III, the general form of the classical system is

(81)

where Ar is a rank-(r+1) tensor. Without any loss of generality, we take Ar to be symmetric in its last r indices. Now suppose that all the Ar are q-sparse in their first index. In other words, for fixed j1,,jr, the number of j0 indices such that Aj0,,jrr0 is at most q. For many systems, q will be small, including the electrostatic Vlasov system [Eq. (1)] discretized on a grid in x and v. The evolution operator M̂ associated with Eq. (81) is not generally sparse when expressed in the occupation basis, even for small q, but it does have a structure that may allow for approximating the eM̂t evolution efficiently; this is explored further in  Appendix C. Still, the costs grow with the size of the entries of M̂, and thus with the entries of Ar. In particular, there is a cost factor (explained in  Appendix C) of

(82)

where a rescaling of the form described earlier [Eqs. (71) and (73)] is assumed to have been applied. Whether Eq. (82) scales as a power of N still depends on the classical system, and on the details of how the system is being varied as N is increased. For systems obtained by discretizing a local, nonlinear partial differential equation, this scaling may be improved by switching to Fourier space (within which the system is no longer local).

Numerical simulations of nonlinear dynamical systems on regular computers have costs at least linear in the number of variables N. Kinetic plasma simulation is especially expensive since very many variables (e.g., phase-space grid points, Fourier modes, or basis functions) are needed to accurately represent the six-dimensional particle distribution function. This motivates investigating whether a quantum computer might be able to perform the same computation with costs sublinear in N. In particular, we have focused on computations that may be possible with costs growing only logarithmically with N.

We have investigated the following approach. The nonlinear dynamical system is first mapped to an infinite-dimensional linear dynamical system using linear embedding. The linear system is then expressed in the occupation basis and truncated to a system of size O[(N+1)m] for a small integer m. If a quantum algorithm can obtain an exponential speedup for the resulting linear computation, then the costs will scale with N as O[mln(N+1)] on a quantum computer. However, whether this logarithmic scaling with N is achievable will depend on the classical dynamical system being studied.

There is no guarantee that a truncated linear system will give a good approximation to the exact nonlinear evolution in general. We noted that the truncation will give correct results up to an order smb in a parameter η associated with the strength of the nonlinearity. This implies that a truncation will be accurate if η is sufficiently small. However, whether that is the case is a difficult question that depends on the details of the problem being studied. Similarly, when η is larger, at what order the terms become negligible, if ever, will be problem dependent. Whether a low-m truncation can approximate the output of any difficult, nonlinear kinetic plasma problem is also an open question.

Additionally, it could be that the nonunitarity of the evolution of the truncated linear system will make quantum computation inefficient. In particular, if the costs end up being logarithmic in N but exponential in the simulation time t, that will typically not give a speedup over the O(Nt) costs for doing the simulation classically. We conjecture that there are certain classes of nonlinear systems for which an exp(t) scaling can be avoided. The structure and conservation properties of physical systems such as the Vlasov–Poisson system might be relevant here, and this is an interesting topic that we leave for future work.

The research was supported in part by the U.S. Department of Energy under Grant No. DE-SC0020393.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

Here we show that a very simple discretization of the electrostatic Vlasov equation [Eq. (1)] satisfies the x·F(x)=0 condition discussed in Sec. II C. Let the variable fj,k represent the value of the distribution function at x=jΔx+x0 and v=kΔv+v0 for some choices of Δx,Δv,x0, and v0. The components of j and k are integers between zero and some maximum, which can be different for each component. Also, for each expression indexed by j and k, there is a corresponding vector representation obtained by collapsing the j and k indices into a single index in some fixed manner. Then we seek to show that f·F(f)=0, which we can do by breaking F(f) up into a number of terms, F(f)=iFi(f), and showing that f·Fi(f)=0 for each i.

Using centered differences, the viif(x,v) term takes the form

(A1)

with ei defined as in Eq. (13). We assume periodic boundary conditions so that the indices are handled cyclically, i.e., ji+1 and ji1 wraparound when they fall outside the range of valid indices. The evaluation of f·Fi(f) for this term gives

(A2)

where the last equality holds because each fj,kfj+ei,k appears twice in the sum with opposite signs and the same prefactor. The key to this cancelation, besides the cyclic handling of the indices, is that the prefactor (kΔv+v0)i/(2Δx) does not depend on the position index ji, which is the only index that differs between fjei,kfj,k and fj,kfj+ei,k.

The same sort of cancelation occurs for the nonlinear terms of Eq. (1), which take the form,

(A3)

This time the prefactor is a function of the variables and is proportional to the electric field component Ei(x). We can express f·Fi(f) for the Eq. (A3) term as

(A4)

where all terms cancel because the electric field does not depend on the velocity index ki. Therefore, f·Fi(f)=0 holds for all the terms in our discretized form of Eq. (1).

We now derive Eq. (60), using the definitions and assumptions from Sec. III. Let M̂0 and M̂1 denote the operators obtained by restricting M̂0 and M̂1, respectively, to the nm subspace. Since M̂0 and M̂1 are finite-dimensional operators with finite entries, they have finite spectral norms, which we denote as Λ0 and Λ1, respectively. The output approximation is

(B1)

where |ψ(0) is the initial state restricted to the nm subspace. c̃j(t,η) is the same as c̃(t) from Sec. III, now just with the η dependence made explicit. Alternatively, we can write c̃(t,η) as a power series in η,

(B2)

If we collect the Eq. (B1) terms that have j factors of η and apply the subadditive and submultiplicative properties of spectral norms, we obtain a bound,

(B3)

from which it follows that the Eq. (B2) series is convergent for all finite t and η. Therefore,

(B4)

holds for finite t and any integer s0. Meanwhile, the exact output quantity c(t,η) [equal to c(t) from Sec. III] can be expressed as the right side of Eq. (B1) without the primes or as

(B5)

but in this case, the prior argument for convergence does not apply since M̂0 and M̂1 do not have finite spectral norms. Instead, we apply a result from the theory of first-order, ordinary differential equations: for an initial value problem given by

(B6)

where F(z,η) is infinitely differentiable with respect to η, if there is a unique η = 0 solution z0(t) for t[0,T] with T finite, then for any integer s0, there exist functions zj(t) such that

(B7)

holds for t[0,T]; this is a particular case of what Hoppensteadt33 calls the regular perturbation theorem. The need to consider a finite-time interval arises because the series in η might not converge in the limit of t. These conditions are met: F(z,η):=Az+ηG(z) is infinitely differentiable in η, the η = 0 solution is just z0(t)=eAtw, and we only consider finite simulation times.

Next, since c(t,η) is assumed to be a specified polynomial of z(t,η), we can plug Eq. (B7) into this polynomial to obtain

(B8)

where the functions cj(t) must be the same as in Eq. (B5) due to the uniqueness of asymptotic expansions. Finally, we use that c̃j(t)=cj(t) for js with the value of s from Eq. (59). This is true because it takes s +1 applications of M̂1 to couple from |c to any component with n > m. Therefore, all terms in

(B9)

that are affected by the truncation [Eq. (59)] have at least s +1 factors of η, which implies that c̃j(t)=cj(t) for js. Combining this result with Eqs. (B4) and (B8) yields Eq. (60). Of course, in practice, we need the approximation to be accurate for particular values of η and t, while Eq. (60) only ensures accuracy for sufficiently small η, in a manner that can depend on t.

Here we outline a strategy for implementing eM̂t for classical systems of the form given in Eq. (81) with each tensor being q-sparse in its first index. Let M denote a matrix representing M̂ in the occupation basis. We consider only Carleman embedding and coherent states embedding. In both cases,

(C1)

ẑj annihilates |n when nj = 0. Additionally, the ẑ operators occur on the right in M̂, and every term has at least one ẑ factor,

(C2)

This allows us to bound the number of different occupation basis components in M̂|n, where |n lies within the truncated subspace, meaning that jnjm. The number of ways that a number of particles between 1 and g can be removed from |n is bounded by

(C3)

For each removal of particles counted in Eq. (C3), there is a corresponding product of ẑ operators which can appear in the terms of M̂ with up to q different ŵ factors. Therefore, the number of occupation basis components in M̂|n, which equals the number of nonzero entries in a column of M, is bounded by

(C4)

Since N does not appear in Eq. (C4), we consider M to be sparse along its columns. This is not generally true along its rows, i.e., there can be poly(N) nonzero entries in a row of M.

Evolution by M can be converted into a matrix inversion problem using the technique of Berry et al.17 The matrix L to be inverted has up to two more nonzero entries in each column than the evolution matrix.17 To perform the matrix inversion using a QLSA, the Hermitian matrix

(C5)

is inverted instead.14 While L is sparse along its columns, L̃ is not sparse. However, if Hamiltonian simulation of L̃ can be performed efficiently, then that can be used to implement an efficient QLSA.14,15

The Hamiltonian simulation algorithm by Low and Chuang34 can be efficient for some nonsparse Hamiltonians. This includes when the Hamiltonian H can be expressed as

(C6)

where the states |χj and |φk can be efficiently prepared, and Γ is not too large: the costs grow linearly with Γ.12,35 We now examine how that can be applied to the simulation of L̃.

Let W denote the size of L and dq(m+1)g+2 the maximum number of nonzero entries in any column of L. We introduce a unitary operation Ô such that

(C7)

where j is the row in L with the lth occurrence of a nonzero entry along column k; if there are not l occurrences or kW, then j can be arbitrary. In Eq. (C7) and below, we take the registers to have dimensions of 2, d, 2W, and W, from left to right. Next, we introduce a unitary operation P̂ that acts as

(C8)

for some constant Λ, where Ljki is the lith nonzero entry along column ki of L, and |ϕj is any superposition of |l|k|x components with k < W. The sets {li} and {ki} depend implicitly on j; in particular, {ki}={k|Ljk0}. Then

(C9)

where the |ϕj states satisfy Ô|1|ϕj=|1|ϕj.

Now we introduce states

(C10)

Also, we define |φj as the state obtained by applying the operation +W(mod2W) to the third register of |χj. Then |χj and |φk satisfy Eq. (C6) with H=L̃ and Γ=dΛ. Moreover, they are straightforward to prepare using an implementation of the ÔP̂ operation.

The P̂ operation requires preparing a superposition containing amplitudes that are proportional to the entries of the classical system tensors Ar. Consequently, the cost to implement the P̂ operation will depend on the details of the classical system. However, there is one result that holds generally: the unitarity of P̂ gives a lower bound on the constant Λ. It cannot be less than the maximum normalization of any row in L. This is the source of the Eq. (82) cost factor. Factors that grow as powers of m and depend on whether Carleman embedding or coherent state embedding is applied have been dropped from Eq. (82) for simplicity.

1.
S.
Murty
, “
Engineering computations at the National Magnetic Fusion Energy Computer Center
,”
Nucl. Technol.-Fusion
4
,
25
32
(
1983
).
2.
R.
Service
, “
Design for U.S. exascale computer takes shape
,”
Science
359
,
617
(
2018
).
3.
J.
Preskill
, “
Quantum computing in the NISQ era and beyond
,”
Quantum
2
,
79
(
2018
).
4.
I.
Joseph
, “
Koopman–von Neumann approach to quantum simulation of nonlinear classical dynamics
,”
Phys. Rev. Res.
2
,
043102
(
2020
).
5.
I. Y.
Dodin
and
E. A.
Startsev
, “
On applications of quantum computing to plasma simulations
,”
Plasma Res. Express
(submitted).
6.
M.
Lubasch
,
J.
Joo
,
P.
Moinier
,
M.
Kiffner
, and
D.
Jaksch
, “
Variational quantum algorithms for nonlinear problems
,”
Phys. Rev. A
101
,
010301
(
2020
).
7.
R.
Steijl
,
Quantum Algorithms for Nonlinear Equations in Fluid Mechanics
(
IntechOpen
,
2020
).
8.
J.-P.
Liu
,
H.
Øie Kolden
,
H. K.
Krovi
,
N. F.
Loureiro
,
K.
Trivisa
, and
A. M.
Childs
, “
Efficient quantum algorithm for dissipative nonlinear differential equations
,” arXiv:2011.03185 (
2020
).
9.
S.
Lloyd
,
G.
De Palma
,
C.
Gokler
,
B.
Kiani
,
Z.-W.
Liu
,
M.
Marvian
,
F.
Tennie
, and
T.
Palmer
, “
Quantum algorithm for nonlinear differential equations
,” arXiv:2011.06571 (
2020
).
10.
Y.
Shi
,
A. R.
Castelli
,
X.
Wu
,
I.
Joseph
,
V.
Geyko
,
F. R.
Graziani
,
S. B.
Libby
,
J. B.
Parker
,
Y. J.
Rosen
,
L. A.
Martinez
, and
J. L.
DuBois
, “
Simulating nonnative cubic interactions on noisy quantum machines
,”
Phys. Rev. Lett.
(submitted).
11.
R. P.
Feynman
, “
Simulating physics with computers
,”
Int. J. Theor. Phys.
21
,
467
488
(
1982
).
12.
A.
Engel
,
G.
Smith
, and
S. E.
Parker
, “
Quantum algorithm for the Vlasov equation
,”
Phys. Rev. A
100
,
062315
(
2019
).
13.
J. B.
Parker
and
I.
Joseph
, “
Quantum phase estimation for a class of generalized eigenvalue problems
,”
Phys. Rev. A
102
,
022422
(
2020
).
14.
A. W.
Harrow
,
A.
Hassidim
, and
S.
Lloyd
, “
Quantum algorithm for linear systems of equations
,”
Phys. Rev. Lett.
103
,
150502
(
2009
).
15.
A.
Childs
,
R.
Kothari
, and
R.
Somma
, “
Quantum algorithm for systems of linear equations with exponentially improved dependence on precision
,”
SIAM J. Comput.
46
,
1920
1950
(
2017
).
16.
A. N.
Soklakov
and
R.
Schack
, “
Efficient state preparation for a register of quantum bits
,”
Phys. Rev. A
73
,
012307
(
2006
).
17.
D. W.
Berry
,
A. M.
Childs
,
A.
Ostrander
, and
G.
Wang
, “
Quantum algorithm for linear differential equations with exponentially improved dependence on precision
,”
Commun. Math. Phys.
356
,
1057
1081
(
2017
).
18.
K.
Kowalski
,
Methods of Hilbert Spaces in the Theory of Nonlinear Dynamical Systems
(
World Scientific
,
Singapore
,
1994
), p.
1
.
19.
B. O.
Koopman
, “
Hamiltonian systems and transformation in Hilbert space
,”
Proc. Natl. Acad. Sci.
17
,
315
318
(
1931
).
20.
J.
von Neumann
, “
Zur operatorenmethode in der klassischen mechanik
,”
Ann. Math.
33
,
587
642
(
1932
).
21.
J.
von Neumann
, “
Zusatze zur arbeit ‘zur operatorenmethode…
,’”
Ann. Math.
33
,
789
791
(
1932
).
22.
T.
Carleman
, “
Application de la théorie des équations intégrales linéaires aux systèmes d'équations différentielles non linéaires
,”
Acta Math.
59
,
63
87
(
1932
).
23.
W.-H.
Steeb
, “
Embedding of nonlinear finite dimensional systems in linear infinite dimensional systems and Bose operators
,”
Hadronic J.
6
,
68
76
(
1983
).
24.
B.
Chirikov
,
F.
Izrailev
, and
D.
Shepelyansky
, “
Quantum chaos: Localization vs. ergodicity
,”
Physica D
33
,
77
88
(
1988
).
25.
T.
Alanson
, “
A ‘quantal’ Hilbert space formulation for nonlinear dynamical systems in terms of probability amplitudes
,”
Phys. Lett. A
163
,
41
45
(
1992
).
26.
K.
Kowalski
, “
Nonlinear dynamical systems and classical orthogonal polynomials
,”
J. Math. Phys.
38
,
2483
2505
(
1997
).
27.
K.
Kowalski
and
W.-H.
Steeb
,
Nonlinear Dynamical Systems and Carleman Linearization
(
World Scientific
,
Singapore
,
1991
).
28.
K.
Kowalski
,
Methods of Hilbert Spaces in the Theory of Nonlinear Dynamical Systems
(
World Scientific
,
Singapore
,
1994
).
29.
V. S.
Varadarajan
,
Geometry of Quantum Theory
(
Van Nostrand Reinhold
,
New York
,
1970
), Vol.
II
.
30.
G.
Brassard
,
P.
Høyer
,
M.
Mosca
, and
A.
Tapp
, “
Quantum amplitude amplification and estimation
,” in
Quantum Computation and Information
, Contemporary Mathematics Vol.
305
(
American Mathematical Society
,
Providence, RI
,
2002
), pp.
53
74
.
31.
L.
Grover
and
T.
Rudolph
, “
Creating superpositions that correspond to efficiently integrable probability distributions
,” arXiv:Quant-ph/0208112 (
2002
).
32.
G. H.
Low
and
I. L.
Chuang
, “
Optimal Hamiltonian simulation by quantum signal processing
,”
Phys. Rev. Lett.
118
,
010501
(
2017
).
33.
F. C.
Hoppensteadt
,
Analysis and Simulation of Chaotic Systems
, 2nd ed. (
Springer
,
New York
,
2000
), pp.
152
153
.
34.
G. H.
Low
and
I. L.
Chuang
, “
Hamiltonian simulation by qubitization
,”
Quantum
3
,
163
(
2019
).
35.
G. H.
Low
and
I. L.
Chuang
, “
Hamiltonian simulation by uniform spectral amplification
,” arXiv:1707.05391 (
2017
).