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 , 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 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 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 , 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 to be normalized. The index j runs over the indices of b, and the state corresponds to a sequence of qubits with values matching the bits of the binary representation of j. Whether the 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 , where denotes the complex vector magnitude.16 The QLSA then outputs a state,
where . It is also assumed that the desired output quantity can be efficiently extracted from . A QLSA can then obtain an exponential speedup over a classical algorithm solving 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 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 ωpe is the inverse time unit, and is the unit in which is expressed, where ne is the electron number density. To solve Eq. (1), 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 for , 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 is the vector of variables and 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 is a vector of operators and 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
The reason for denoting these operators by is revealed by the following evaluation:
and therefore
Note that, since the operators have eigenvalues of every variable value taken at any time, continuous evolution of implies that the space spanned by the states is infinite-dimensional. Also, Eq. (6) implies that the operators commute within this space: for every 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 are analytic functions. Then is defined by replacement of z with in the power series representation of . Equation (7) shows that the states evolve linearly according to
where
The finite-time evolution can be expressed as
where depends on the initial state of the classical system, .
Viewing the linear evolution as occurring in a Hilbert space containing the states, the inner product with some state 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 , 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
where
It is straightforward to check that these operators satisfy both parts of Eq. (4) with being the basis state. Additionally, the operators satisfy
which can be used to simplify some expressions.
To concisely express the states, we first write as
For instance, with N = 2, we would have . Now
Additional insight can be obtained after defining a number operator by
Using n to denote the eigenvalues, states can be broken up into components of different n. For example, the n = 1 component of is
and the n = 2 component is
Another observation is that decreases n by one, except the n = 0 component, , 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,
Similarly, a linear combination of variables
is equivalent to with
Polynomials of the variables can also be obtained by adding n >1 components to . For instance, with . More generally, for any polynomial of the variables, there is a specific state for which evaluates that polynomial.
B. Coherent states embedding
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
which amounts to the standard commutation relations and the statement that is the ground state. Additionally, we can use to denote the occupation basis and express the number operator as
Now the eigenvalues n of 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 factor in Eq. (26) can have an effect. For example, yields rather than .
While the states 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 [Eqs. (16) or (26)] changes in normalization when 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 and , where and are dimensionless versions of canonical position and momentum operators, respectively. Then
and both parts of Eq. (4) are met with , the 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 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 must be a real function as well. The states are position eigenstates, and they can be expressed as
Since the translation operator is unitary, the states have the same normalization as .
The evolution operator is
To study this embedding in the occupation basis , we introduce
It is easily checked that these obey
and thus they are standard bosonic lowering operators. Using the expressions for and in the occupation basis, the evolution operator
can also be expressed in an occupation basis.
To express [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 , the required result is well known: the one-dimensional QHO eigenstates in the position basis are
where are the physicists' Hermite polynomials. The form of follows directly:
where we have chosen to drop the unimportant factors of .
Output quantities work somewhat differently in this version of linear embedding. The components of in the subspace, where n denotes the eigenvalues of the total number operator , are polynomials of the variables of degree times a factor of . Therefore, if we write
for a state in the subspace, then is a polynomial of the variables with degree . For example, yields . A potentially simpler way to express an output quantity is as the expectation value of an observable: for some analytic function . However, only can be evaluated within the 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 values, not any value. For this reason, we exclusively study the way of expressing outputs.
The factor in Eq. (36) can complicate the representation of desired outputs, which are frequently polynomials of the variables. However, this issue is avoided if for all x. This ensures that , so the factor in Eq. (36) is constant. Additionally, this has a consequence for the form of the evolution operator: for to hold for all x, it must vanish identically, meaning that all terms are canceled. Therefore,
holds for any vector of commuting operators. Now, [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 [Eq. (33)] have at least one raising and lowering operator. Also, supposing that is a polynomial in x of degree g, then each term has a maximum of g + 1 operators. Therefore, couples between occupation basis components for which n, the eigenvalue, differs by at most g − 1.
The assumption that for all x is nontrivial, but it holds in some important scenarios. For instance, it holds when
for any , 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 that are not in the subset are taken to be zero. Second, can be any linear combination of terms of the described form, and each term can have a different . 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 , 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 condition. Therefore, the factor of is a constant, and it is convenient to divide this out of the state. Then becomes
and the output quantities [Eq. (36)] are updated to
where, as before, is a polynomial of the variables with a form determined by the state .
Evolution in the position-space embedding can also be described in terms of a “Hamiltonian” , and can be written as
It is well known and easily verified that, for operators and satisfying ,
for any analytic function of . Applying this, we find that
where is defined by replacing x with in . For many systems, 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 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 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
which must be applied to obtain an output of the form, yet Eq. (45) cannot be evaluated in general without the full solution to the system. We assume that the described extension is applied when necessary to get a system for which always vanishes. Then,
In the occupation basis, Eq. (46) becomes
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.
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 . If we introduce operators and and a state 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 with a finite difference matrix in the position basis. In this case, the dimensionality of the truncated space scales as , where ϵ is the chosen precision. Therefore, the space complexity, i.e., the required number of qubits, scales as . 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 is a polynomial of degree g in the variables for some small integer g and that , i.e., 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 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 belonging to the subspace of the occupation basis. As usual, n is used to denote the eigenvalues of the total number operator . Second, the evolution operator 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 is purely nonlinear in z. In particular, is a degree g polynomial of the variables. Now the evolution operator can be decomposed as
By itself, generates the linearized evolution of the original dynamical system, i.e., the evolution in the limit of . Since g = 1 for linear evolution, the terms in do not couple between different n. The need to consider components comes from , 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 . Let denote the approximation to c(t) that is obtained by performing the embedded evolution in the truncated space and evaluating 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 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 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,
Noting that
the number of qubits required to represent the 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 and denote and , respectively, after restriction to the truncated subspace. Then the computation can be expressed as
where and are unitary operations that prepare states proportional to and , 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 state. The technique of amplitude estimation can then be applied to estimate this component using 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 . However, the full complex value can also be estimated using a simple algorithm extension.12
The state is a sum of components with particle counts up to m. To represent this using qubits, we can use m registers, each with 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 state is a sum of components (some of which may vanish) with particle counts up to [Eq. (59)]. Therefore, can be represented in the same way as .
Overall efficiency requires efficient preparation of states proportional to and . In particular, we want these operations to cost . 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 qubits with the values of this register representing particle counts. Then, for each particle count , the corresponding component of [e.g., Eq. (19) with t = 0] or [e.g., Eq. (22)] is prepared. Whether that can be done with complexity will depend on the details of the initial state and output quantity, but it is possible in many cases, including when the maximum is a factor of larger than the average , 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 and . A factor , where and , 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, 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 . Equation (68) is minimized by for all j, and the value at the minimum is . Consequently, an m = 2 truncation has . More generally, the scaling of with N is . 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 .
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 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 . 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 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 , to bound the initial state normalization as
Thus, is achieved.
Applying rescaling [Eq. (71)] to Carleman embedding or coherent states embedding results in a factor of γn being applied to components of for each particle count . 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 is rescaled by for each . The general form of the output quantity is
where Cn is a rank-n tensor. Suppose that, prior to any rescaling of the variables, the variable values are bounded as , while the entries of Cn scale as . 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 , and a rescaling with will result in and . Now, the components of in the occupation basis are proportional to the entries of the C tensors (in a manner independent of N). Therefore, the scaling of with N can be obtained by summing all squared absolute entries of Cn for each n, which yields 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 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 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 [Eq. (58)] is formally anti-Hermitian,
Additionally, if we order the basis states based on their particle count n, from low to high, then takes on a block-upper-triangular structure, with one block for each n. This occurs because always decreases n, while leaves n constant. All the blocks along the diagonal are anti-Hermitian since they are components of . Then a block-diagonal, unitary transformation can be applied to diagonalize 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 obtained by truncating to the subspace, which implies that the eigenvalues of are purely imaginary. As a consequence, the condition number κ of the evolution will remain bounded as if is diagonalizable. If happens to be nondiagonalizable, the block structure of limits the size of any Jordan chain to m, which leads to a bound of . 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 evolution scale with N. With the same assumptions as in Sec. III, the general form of the classical system is
where Ar is a rank- 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 , the number of j0 indices such that 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 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 evolution efficiently; this is explored further in Appendix C. Still, the costs grow with the size of the entries of , and thus with the entries of Ar. In particular, there is a cost factor (explained in Appendix C) of
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 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 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 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 costs for doing the simulation classically. We conjecture that there are certain classes of nonlinear systems for which an 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 condition discussed in Sec. II C. Let the variable represent the value of the distribution function at and for some choices of , and . The components of and are integers between zero and some maximum, which can be different for each component. Also, for each expression indexed by and , there is a corresponding vector representation obtained by collapsing the and indices into a single index in some fixed manner. Then we seek to show that , which we can do by breaking up into a number of terms, , and showing that for each i.
Using centered differences, the term takes the form
with ei defined as in Eq. (13). We assume periodic boundary conditions so that the indices are handled cyclically, i.e., and wraparound when they fall outside the range of valid indices. The evaluation of for this term gives
where the last equality holds because each 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 does not depend on the position index ji, which is the only index that differs between and .
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 . We can express for the Eq. (A3) term as
where all terms cancel because the electric field does not depend on the velocity index ki. Therefore, 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 and denote the operators obtained by restricting and , respectively, to the subspace. Since and are finite-dimensional operators with finite entries, they have finite spectral norms, which we denote as and , respectively. The output approximation is
where is the initial state restricted to the subspace. is the same as from Sec. III, now just with the η dependence made explicit. Alternatively, we can write 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 . Meanwhile, the exact output quantity [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 and 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 is infinitely differentiable with respect to η, if there is a unique η = 0 solution for with T finite, then for any integer , there exist functions such that
holds for ; 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 . These conditions are met: is infinitely differentiable in η, the η = 0 solution is just , and we only consider finite simulation times.
Next, since is assumed to be a specified polynomial of , we can plug Eq. (B7) into this polynomial to obtain
where the functions must be the same as in Eq. (B5) due to the uniqueness of asymptotic expansions. Finally, we use that for with the value of s from Eq. (59). This is true because it takes s + 1 applications of to couple from 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 for . 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 for classical systems of the form given in Eq. (81) with each tensor being q-sparse in its first index. Let denote a matrix representing in the occupation basis. We consider only Carleman embedding and coherent states embedding. In both cases,
annihilates when nj = 0. Additionally, the operators occur on the right in , and every term has at least one factor,
This allows us to bound the number of different occupation basis components in , where lies within the truncated subspace, meaning that . The number of ways that a number of particles between 1 and g can be removed from is bounded by
For each removal of particles counted in Eq. (C3), there is a corresponding product of operators which can appear in the terms of with up to q different factors. Therefore, the number of occupation basis components in , which equals the number of nonzero entries in a column of , is bounded by
Since N does not appear in Eq. (C4), we consider to be sparse along its columns. This is not generally true along its rows, i.e., there can be nonzero entries in a row of .
Evolution by 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, is not sparse. However, if Hamiltonian simulation of 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
where the states and 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 .
Let W denote the size of L and the maximum number of nonzero entries in any column of L. We introduce a unitary operation such that
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 , 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 that acts as
for some constant Λ, where is the lith nonzero entry along column ki of L, and is any superposition of components with k < W. The sets and depend implicitly on j; in particular, . Then
where the states satisfy .
Now we introduce states
Also, we define as the state obtained by applying the operation to the third register of . Then and satisfy Eq. (C6) with and . Moreover, they are straightforward to prepare using an implementation of the operation.
The 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 operation will depend on the details of the classical system. However, there is one result that holds generally: the unitarity of 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.