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.

## I. INTRODUCTION

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 Joseph^{13} 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 $|\psi i\u27e9$ to be normalized. The index *j* runs over the indices of **b**, and the state $|j\u27e9$ corresponds to a sequence of qubits with values matching the bits of the binary representation of *j*. Whether the $|\psi i\u27e9$ state can be prepared efficiently depends on **b**, but it is possible in some general cases, such as if each *b _{j}* can be computed efficiently and $maxj|bj|/|b|=O(1/N)$, where $|b|=\u2211jbj*bj$ denotes the complex vector magnitude.

^{16}The QLSA then outputs a state,

where $x=A\u22121b$. It is also assumed that the desired output quantity can be efficiently extracted from $|\psi f\u27e9$. 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,

where we have also rescaled quantities so that they are all dimensionless. The Debye length *λ _{De}* is the distance unit, the plasma frequency

*ω*is the inverse time unit, and $ne/(\lambda De\omega pe)3$ is the unit in which $f(x,v)$ is expressed, where

_{pe}*n*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.

_{e}One way to simulate nonlinear dynamics on a quantum computer is with the Koopman–von Neumann approach introduced by Joseph^{4} 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,\u2026$, 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.

## II. LINEAR EMBEDDING

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

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,

where $w\u0302$ is a vector of operators and $|0\u27e9$ 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 $z\u0302$ that satisfies

The reason for denoting these operators by $z\u0302$ is revealed by the following evaluation:

and therefore

Note that, since the $z\u0302$ operators have eigenvalues of every variable value taken at any time, continuous evolution of $z(t)$ implies that the space spanned by the $|\psi (t)\u27e9$ states is infinite-dimensional. Also, Eq. (6) implies that the $z\u0302$ operators commute within this space: $z\u0302jz\u0302k|\psi (t)\u27e9=z\u0302kz\u0302j|\psi (t)\u27e9$ for every $|\psi (t)\u27e9$ 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

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

where

The finite-time evolution can be expressed as

where $|\psi (0)\u27e9=ez(0)\xb7w\u0302|0\u27e9$ depends on the initial state of the classical system, $z(0)$.

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

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.

### A. Carleman embedding

Let $|n\u27e9$, 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 $z\u0302$ and $w\u0302$ are

where

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

which can be used to simplify some expressions.

To concisely express the $|\psi (t)\u27e9$ states, we first write $|n\u27e9$ as

For instance, with *N *=* *2, we would have $|n\u27e9=|n1\u27e9\u2297|n2\u27e9$. Now

Additional insight can be obtained after defining a number operator $n\u0302$ by

Using *n* to denote the $n\u0302$ eigenvalues, states can be broken up into components of different *n*. For example, the *n *=* *1 component of $|\psi (t)\u27e9$ is

and the *n *=* *2 component is

Another observation is that $z\u0302$ decreases *n* by one, except the *n *=* *0 component, $|0\u27e9$, which it annihilates; and $w\u0302$ increases *n* by one. Those familiar with the technique of Carleman linearization^{22,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,

Similarly, a linear combination of variables

is equivalent to $\u27e8c|\psi (t)\u27e9$ with

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

### B. Coherent states embedding

Now suppose that $z\u0302=a\u0302$ and $w\u0302=a\u0302\u2020$, where $a\u0302$ are standard bosonic lowering operators. This version of linear embedding has been extensively explored by Kowalski.^{28} In this case, Eq. (4) becomes

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

Now the eigenvalues *n* of $n\u0302$ represent numbers of bosonic particles.

In this version of linear embedding, the evolution operator is

and the states are

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

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\u27e9=|2ej\u27e9$ yields $\u27e8c|\psi (t)\u27e9=zj2(t)/2$ rather than $zj2(t)$.

While the states $|\psi (t)\u27e9$ 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 $|\psi (t)\u27e9$ [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.

### C. Position-space embedding

This time, take $z\u0302=x\u0302$ and $w\u0302=\u2212ip\u0302$, where $x\u0302$ and $p\u0302$ are dimensionless versions of canonical position and momentum operators, respectively. Then

and both parts of Eq. (4) are met with $|0\u27e9=|x=0\u27e9$, 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\u0302$ 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,\u2009F(\xb7)$ must be a real function as well. The $|\psi (t)\u27e9$ states are position eigenstates, and they can be expressed as

Since the translation operator $e\u2212ix(t)\xb7p\u0302$ is unitary, the $|\psi (t)\u27e9$ states have the same normalization as $|x=0\u27e9$.

The evolution operator is

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

It is easily checked that these $a\u0302$ obey

and thus they are standard bosonic lowering operators. Using the expressions for $a\u0302$ and $a\u0302\u2020$ in the occupation basis, the evolution operator

can also be expressed in an occupation basis.

To express $|\psi (t)\u27e9$ [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 $\u210f=m\omega =1$, the required result is well known: the one-dimensional QHO eigenstates in the position basis are

where $Hn(x)$ are the physicists' Hermite polynomials. The form of $|\psi (t)\u27e9\u221d|x(t)\u27e9$ follows directly:

where we have chosen to drop the unimportant factors of $\pi \u22121/4$.

Output quantities work somewhat differently in this version of linear embedding. The components of $|\psi (t)\u27e9$ in the $n\u2264b$ subspace, where *n* denotes the eigenvalues of the total number operator $n\u0302$, are polynomials of the variables of degree $\u2264b$ times a factor of $e\u2212x\xb7x/2$. Therefore, if we write

for $|c\u27e9$ a state in the $n\u2264b$ subspace, then $p(x)$ is a polynomial of the variables with degree $\u2264b$. For example, $|c\u27e9=|ej\u27e9$ yields $p[x(t)]=2\u2009xj(t)$. A potentially simpler way to express an output quantity is as the expectation value of an observable: $\u27e8g\u27e9:=\u27e8\psi (t)|g(x\u0302)|\psi (t)\u27e9$ for some analytic function $g(x)$. However, only $\u27e8c|\psi (t)\u27e9$ can be evaluated within the $n\u2264b$ 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 $\u27e8c|\psi (t)\u27e9$ values, not any $\u27e8g\u27e9$ value. For this reason, we exclusively study the $\u27e8c|\psi (t)\u27e9$ way of expressing outputs.

The $e\u2212x\xb7x/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\xb7F(x)=0$ for all **x**. This ensures that $dt(x\xb7x)=0$, so the $e\u2212x\xb7x/2$ factor in Eq. (36) is constant. Additionally, this has a consequence for the form of the evolution operator: for $x\xb7F(x)=0$ to hold for all **x**, it must vanish identically, meaning that all terms are canceled. Therefore,

holds for $r\u0302$ any vector of commuting operators. Now, $M\u0302$ [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

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\u0302$ [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\u0302$ couples between occupation basis components for which *n*, the $n\u0302$ eigenvalue, differs by at most *g* − 1.

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

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\xb7F(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\xb7F(x)=0$ condition. Therefore, the factor of $e\u2212x\xb7x/2$ is a constant, and it is convenient to divide this out of the state. Then $|\psi (t)\u27e9$ becomes

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

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

Evolution in the position-space embedding can also be described in terms of a “Hamiltonian” $H\u0302=iM\u0302$, and $H\u0302$ can be written as

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

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

where $div\u2009F(x\u0302)$ is defined by replacing **x** with $x\u0302$ in $div\u2009F(x):=\u2211j\u2202Fj(x)/\u2202xj$. For many systems, $div\u2009F(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 *x _{j}*, we replace that

*x*occurrence with the corresponding variable from the other system. This has no impact on the dynamics, yet it results in the elimination of $div\u2009F(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 Kowalski

_{j}^{26}) since the latter method introduces a factor of

which must be applied to obtain an output of the $\u27e8c|\psi (t)\u27e9$ 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 $div\u2009F(x)$ always vanishes. Then,

In the occupation basis, Eq. (46) becomes

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

### D. Continuum limit

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,

where *F* can contain derivatives and integrals of the function *f*(*t*) with respect to the coordinates $q$. If we introduce operators $f\u0302q$ and $h\u0302q$ and a state $|0\u27e9$ such that

then the states

are found to satisfy

and

where

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

and

## III. TRUNCATION OF THE SPACE

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\u0302$ with a finite difference matrix in the position basis. In this case, the dimensionality of the truncated space scales as $(1/\u03f5)N$, where *ϵ* is the chosen precision. Therefore, the space complexity, i.e., the required number of qubits, scales as $N\u2009log2(1/\u03f5)$. 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(\xb7)$ is a polynomial of degree *g* in the variables for some small integer *g* and that $F(0)=0$, i.e., $F(\xb7)$ 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\xb7F(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

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

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

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

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

Now, we truncate the space such that only the

subspace is retained for some integer $s\u22650$. Let $c\u0303(t)$ denote the approximation to *c*(*t*) that is obtained by performing the embedded evolution in the truncated space and evaluating $\u27e8c|\psi (t)\u27e9$ after some fixed, finite time *t*. Then,

We derive Eq. (60) in Appendix B. Treating *η* as a dimensionless parameter characterizing the strength of the nonlinearity, we expect $c(t)\u2248c\u0303(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 $n\u2264m$ subspace, consider *N *+* *1 bins into which *m* particles are placed. The first *N* bins correspond to the numbers *n _{j}* 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,

Noting that

the number of qubits required to represent the $n\u2264m$ subspace is

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

## IV. EFFICIENCY

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

where $V\u0302\psi \u2032$ and $V\u0302c$ are unitary operations that prepare states proportional to $|\psi (0)\u2032\u27e9$ and $|c\u27e9$, 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\u27e9$ state. The technique of amplitude estimation can then be applied to estimate this component using $O(1/\epsilon )$ 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\u0303(t)|$. However, the full complex value can also be estimated using a simple algorithm extension.^{12}

The $|\psi (0)\u2032\u27e9$ state is a sum of components with particle counts up to *m*. To represent this using qubits, we can use *m* registers, each with $\u2308\u2009log2(N+1)\u2309$ 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 *m* − *n* unused registers which should be set to the extra basis state. This representation achieves the Eq. (64) scaling. The $|c\u27e9$ state is a sum of components (some of which may vanish) with particle counts up to $b\u2264m$ [Eq. (59)]. Therefore, $|c\u27e9$ can be represented in the same way as $|\psi (0)\u2032\u27e9$.

Overall efficiency requires efficient preparation of states proportional to $|\psi (0)\u2032\u27e9$ and $|c\u27e9$. In particular, we want these operations to cost $poly(log\u2009N)$. 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 $\u2308\u2009log2(m+1)\u2309$ qubits with the values of this register representing particle counts. Then, for each particle count $n\u2264m$, the corresponding component of $|\psi (0)\u2032\u27e9$ [e.g., Eq. (19) with *t *=* *0] or $|c\u27e9$ [e.g., Eq. (22)] is prepared. Whether that can be done with $poly(log\u2009N)$ 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 *b _{j}* is the prepared amplitude for state component

*j*, assuming that each

*b*can be computed efficiently using the index

_{j}*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 $|\psi (0)\u2032\u27e9$ and $|c\u27e9$. A factor $\chi \zeta $, where $\chi :=\u27e8\psi (0)\u2032|\psi (0)\u2032\u27e9$ and $\zeta :=\u27e8c|c\u27e9$, 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

and to extract the result to within error *δ* using amplitude estimation, $O(\chi \zeta /\delta )$ 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

and those contribute

to $\u27e8\psi |\psi \u27e9$. 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 $\chi 2>N/4$. More generally, the scaling of $\chi 2$ with *N* is $1/\chi 2=O[N\u2212\u230am/2\u230b]$. 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(\chi \zeta /\delta )$.

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

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

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

with a constant

so that $\gamma |z(0)|=O(1)$. The transformed system is

and the truncated initial state normalization becomes

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

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

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

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

Thus, $\chi =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 $|\psi (0)\u2032\u27e9$ for each particle count $n\u2264m$. 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\u27e9$ is rescaled by $\gamma \u2212n$ for each $n\u2264b$. The general form of the output quantity is

where *C ^{n}* is a rank-

*n*tensor. Suppose that, prior to any rescaling of the variables, the variable values are bounded as $N\u2192\u221e$, while the entries of

*C*scale as $O(N\u2212n)$. For example, the variables can represent a bounded distribution function over a grid of size

^{n}*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 $\gamma \u221d1/N$ will result in $\chi =O(1)$ and $Cj1,\u2026,jnn=O(N\u2212n/2)$. Now, the components of $|c\u27e9$ in the occupation basis are proportional to the entries of the

*C*tensors (in a manner independent of

*N*). Therefore, the scaling of $\zeta 2$ with

*N*can be obtained by summing all squared absolute entries of

*C*for each

^{n}*n*, which yields $\zeta =O(1)$ with respect to

*N*in this case.

So Carleman embedding and coherent states embedding can potentially avoid having the cost factor $\chi \zeta $ grow as a power of *N*, but with these embeddings the linear evolution given by $eM\u0302\u2032t$ 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\u2009(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\u03020$ [Eq. (58)] is formally anti-Hermitian,

Additionally, if we order the basis states based on their particle count *n*, from low to high, then $M\u0302$ takes on a block-upper-triangular structure, with one block for each *n*. This occurs because $M\u03021$ always decreases *n*, while $M\u03020$ leaves *n* constant. All the blocks along the diagonal are anti-Hermitian since they are components of $M\u03020$. Then a block-diagonal, unitary transformation can be applied to diagonalize $M\u03020$ 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\u0302\u2032$ obtained by truncating $M\u0302$ to the $n\u2264m$ subspace, which implies that the eigenvalues of $M\u0302\u2032$ are purely imaginary. As a consequence, the condition number *κ* of the evolution $eM\u0302\u2032t$ will remain bounded as $t\u2192\u221e$ if $M\u0302\u2032$ is diagonalizable. If $M\u0302\u2032$ happens to be nondiagonalizable, the block structure of $M\u0302\u2032$ limits the size of any Jordan chain to *m*, which leads to a bound of $\kappa =O(t2(m\u22121))$. 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\u0302\u2032t$ evolution scale with *N*. With the same assumptions as in Sec. III, the general form of the classical system is

where *A ^{r}* is a rank-$(r+1)$ tensor. Without any loss of generality, we take

*A*to be symmetric in its last

^{r}*r*indices. Now suppose that all the

*A*are

^{r}*q*-sparse in their first index. In other words, for fixed $j1,\u2026,jr$, the number of

*j*

_{0}indices such that $Aj0,\u2026,jrr\u22600$ 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\u0302\u2032$ 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\u0302\u2032t$ evolution efficiently; this is explored further in Appendix C. Still, the costs grow with the size of the entries of $M\u0302\u2032$, and thus with the entries of

*A*. In particular, there is a cost factor (explained in Appendix C) of

^{r}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).

## V. DISCUSSION

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[m\u2009ln\u2009(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 $s\u221dm\u2212b$ 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\u2009(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.

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

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

### APPENDIX A: VLASOV DISCRETIZATION

Here we show that a very simple discretization of the electrostatic Vlasov equation [Eq. (1)] satisfies the $x\xb7F(x)=0$ condition discussed in Sec. II C. Let the variable $fj,k$ represent the value of the distribution function at $x=j\Delta x+x0$ and $v=k\Delta v+v0$ for some choices of $\Delta x,\u2009\Delta v,\u2009x0$, 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\xb7F(f)=0$, which we can do by breaking $F(f)$ up into a number of terms, $F(f)=\u2211iFi(f)$, and showing that $f\xb7Fi(f)=0$ for each *i*.

Using centered differences, the $\u2212vi\u2207if(x,v)$ term takes the form

with **e**_{i} defined as in Eq. (13). We assume periodic boundary conditions so that the indices are handled cyclically, i.e., $ji+1$ and $ji\u22121$ wraparound when they fall outside the range of valid indices. The evaluation of $f\xb7Fi(f)$ for this term gives

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\Delta v+v0)i/(2\Delta x)$ does not depend on the position index *j _{i}*, which is the only index that differs between $fj\u2212ei,kfj,k$ and $fj,kfj+ei,k$.

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

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

where all terms cancel because the electric field does not depend on the velocity index *k _{i}*. Therefore, $f\xb7Fi(f)=0$ holds for all the terms in our discretized form of Eq. (1).

### APPENDIX B: NONLINEAR EXPANSION

We now derive Eq. (60), using the definitions and assumptions from Sec. III. Let $M\u0302\u20320$ and $M\u0302\u20321$ denote the operators obtained by restricting $M\u03020$ and $M\u03021$, respectively, to the $n\u2264m$ subspace. Since $M\u0302\u20320$ and $M\u0302\u20321$ are finite-dimensional operators with finite entries, they have finite spectral norms, which we denote as $\Lambda 0$ and $\Lambda 1$, respectively. The output approximation is

where $|\psi \u2032(0)\u27e9$ is the initial state restricted to the $n\u2264m$ subspace. $c\u0303j(t,\eta )$ is the same as $c\u0303(t)$ from Sec. III, now just with the *η* dependence made explicit. Alternatively, we can write $c\u0303(t,\eta )$ as a power series in *η*,

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,

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

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

but in this case, the prior argument for convergence does not apply since $M\u03020$ and $M\u03021$ 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

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

holds for $t\u2208[0,T]$; this is a particular case of what Hoppensteadt^{33} 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\u2192\u221e$. These conditions are met: $F(z,\eta ):=Az+\eta 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,\eta )$ is assumed to be a specified polynomial of $z(t,\eta )$, we can plug Eq. (B7) into this polynomial to obtain

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\u0303j(t)=cj(t)$ for $j\u2264s$ with the value of *s* from Eq. (59). This is true because it takes *s *+* *1 applications of $M\u03021$ to couple from $|c\u27e9$ to any component with *n *>* m*. Therefore, all terms in

that are affected by the truncation [Eq. (59)] have at least *s *+* *1 factors of *η*, which implies that $c\u0303j(t)=cj(t)$ for $j\u2264s$. 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*.

### APPEDNIX C: HAMILTONIAN SIMULATION

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

$z\u0302j$ annihilates $|n\u27e9$ when *n _{j}* = 0. Additionally, the $z\u0302$ operators occur on the right in $M\u0302$, and every term has at least one $z\u0302$ factor,

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

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

Since *N* does not appear in Eq. (C4), we consider $M\u2032$ 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\u2032$.

Evolution by $M\u2032$ 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

is inverted instead.^{14} While *L* is sparse along its columns, $L\u0303$ is not sparse. However, if Hamiltonian simulation of $L\u0303$ can be performed efficiently, then that can be used to implement an efficient QLSA.^{14,15}

The Hamiltonian simulation algorithm by Low and Chuang^{34} can be efficient for some nonsparse Hamiltonians. This includes when the Hamiltonian *H* can be expressed as

where the states $|\chi j\u27e9$ and $|\phi k\u27e9$ 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\u0303$.

Let *W* denote the size of *L* and $d\u2264q(m+1)g+2$ the maximum number of nonzero entries in any column of *L*. We introduce a unitary operation $O\u0302$ such that

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

for some constant Λ, where $Ljki$ is the *l _{i}*th nonzero entry along column

*k*of

_{i}*L*, and $|\varphi j\u27e9$ is any superposition of $|l\u27e9|k\u27e9|x\u27e9$ components with

*k*<

*W*. The sets ${li}$ and ${ki}$ depend implicitly on

*j*; in particular, ${ki}={k|Ljk\u22600}$. Then

where the $|\varphi \u2032j\u27e9$ states satisfy $O\u0302|1\u27e9|\varphi j\u27e9=|1\u27e9|\varphi \u2032j\u27e9$.

Now we introduce states

Also, we define $|\phi j\u27e9$ as the state obtained by applying the operation $+W\u2009(mod\u20092W)$ to the third register of $|\chi j\u27e9$. Then $|\chi j\u27e9$ and $|\phi k\u27e9$ satisfy Eq. (C6) with $H=L\u0303$ and $\Gamma =d\Lambda $. Moreover, they are straightforward to prepare using an implementation of the $O\u0302P\u0302$ operation.

The $P\u0302$ operation requires preparing a superposition containing amplitudes that are proportional to the entries of the classical system tensors *A ^{r}*. Consequently, the cost to implement the $P\u0302$ operation will depend on the details of the classical system. However, there is one result that holds generally: the unitarity of $P\u0302$ 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.