Koopman-von Neumann Approach to Quantum Simulation of Nonlinear Classical Dynamics

Quantum computers can be used to simulate nonlinear non-Hamiltonian classical dynamics on phase space by using the generalized Koopman-von Neumann formulation of classical mechanics. The Koopman-von Neumann formulation implies that the conservation of the probability distribution function on phase space, as expressed by the Liouville equation, can be recast as an equivalent Schr\"odinger equation on Hilbert space with a Hermitian Hamiltonian operator and a unitary propagator. This Schr\"odinger equation is linear in the momenta because it derives from a constrained Hamiltonian system with twice the classical phase space dimension. A quantum computer with finite resources can be used to simulate a finite-dimensional approximation of this unitary evolution operator. Quantum simulation of classical dynamics is exponentially more efficient than an Eulerian discretization of the Liouville equation if the Koopman-von Neumann Hamiltonian is sparse. It is possible that utilizing quantum walks and associated techniques could lead to a quadratic improvement over classical time-dependent Monte Carlo methods.


A. Motivation
In principle, future error-corrected quantum computers have the power to simulate quantum mechanical systems exponentially more efficiently [1,2] than computers that are bound to satisfy the laws of classical physics. The recent growth in the capabilities of today's quantum computing devices has spurred great interest in quantum simulation and they have been used to perform key demonstrations of quantum calculations [3][4][5]. Yet, for many fields of science and engineering, including biology, chemistry, and physics, a large share of today's computational resources are used for the simulation of classical dynamics. Hence, it is important to understand whether quantum computers can provide similar gains in efficiency for the simulation of classical dynamics.
Classical dynamical systems are typically nonlinear and many important examples are not Hamiltonian. In fact, they are often dissipative. Since quantum computers can only perform linear unitary operations, it is not clear how nonlinear nonunitary simulations can be performed efficiently. While efficient quantum algorithms for linear ordinary differential equations (ODEs) are known [6,7], an attempt to simulate nonlinear dynamics by measuring the state at each time step and feeding this information into the next time step is known to be rather inefficient [8]. * joseph5@llnl.gov

B. Comparison of Classical and Quantum Resource Requirements
Often one simulates a finite number of trajectories of a classical dynamical system in order to infer statistics of the evolution of the probability distribution function (PDF), f . A probabilistic description of a classical system is actually similar in complexity to that of a quantum system. If there are n classical bits, then describing the classical PDF, over all N = 2 n possible states |j requires the specification of N − 1 real numbers: the f j subject to the constraint j f j = 1. The wavefunction for a pure quantum state holds twice as much information due to the fact that the probability amplitudes are complex. Specifying the wavefunction, requires the specification of N − 1 complex numbers: the ψ j subject to the constraint j |ψ j | 2 = 1 and, since the overall phase does not matter, a constraint on the phase, such as Im (ψ 0 ) = 0; i.e. 2(N − 1) real numbers. Hence, storing the memory and performing an operation on each component of a probabilistic classical system requires half of the resources of that of a pure quantum state.
Specifying a mixed quantum state via the density matrix 2 requires specifying N 2 − 1 real numbers: the ρ jk subject to the constraints of Hermiticity, ρ = ρ † , and unit trace, tr ρ = 1 (and ρ must also be positive semi-definite). The diagonal entries of ρ are sufficient to describe f , which implies that a quantum simulation that experiences decoherence by the end of the calculation could still potentially generate a useful simulation of the PDF. On the other hand, a classical Hamiltonian system is defined on twice the phase space dimension of the quantized version of the classical system, because it includes both configuration space coordinates and conjugate momentum (or velocity) coordinates, in which case, representing the classical system also requires N 2 − 1 real numbers. Thus, the resource requirements for probabilistic classical systems and quantum systems are quite similar.
A clear example of this fact is that simulating the Schrödinger equation requires twice as many resources to store the complex amplitude as simulating the diffusion equation for a real density. Hence, a "Wick rotation" is often employed to convert between the Schrödinger equation and the diffusion equation. For example, this Wick rotation is used to convert between the temporal propagator of quantum mechanics and the thermal partition function of statistical mechanics [9], which is closely related to the propagator of the diffusion equation. It is also used by a number of quantum Monte Carlo algorithms [10] that seek to find eigenstates of the Schrödinger equation by simulating the diffusion equation.
As another example, consider simulation of the Liouville equation, which describes the conservation of the classical probability distribution function on phase space and provides the foundation for nonequilibrium statistical mechanics. Such simulations are commonly performed in diverse fields of science such as population biology, gravitation, molecular dynamics, and plasma physics. If there are M particles traveling in d dimensions, then there are M d classical degrees of freedom, and, for a Hamiltonian system, the Liouville equation is an D = 2M d dimensional partial differential equation (PDE). If one employs an Eulerian discretization of phase space, by using L spatial grid points in each dimension, then this requires an exponential amount of memory and computational work, ∼ L D , to process the data.
Notably, Ref. [11] devised a quantum algorithm for solving the linearized Vlasov-Maxwell system that claims to allow the simulation of Landau damping with exponentially reduced computational work.

C. Quantization of Classical Dynamics
There are different strategies that one might potentially employ to represent classical dynamics via quantum mechanics. Perhaps the most obvious method is to use any valid quantization that reduces to the classical system in the limit of vanishing Planck's constant. The quantum dynamics will track the semi-classical dynamics for arbitrarily long times in the limit that → 0. It is known that simulating the quantized version of a classical system can be used to compute dynamical quantities, such as the diffusion coefficient and Lyapunov exponents, more efficiently [12,13]. This approach has the benefit that the Schrödinger PDE is defined on 1/2 of the classical phase space dimension, and, hence, is technically cheaper to simulate than the Liouville PDE. It also shares similarities with the manner in which quantum walks can be used to accelerate classical Monte Carlo methods in order to generically achieve a quadratic speedup [14,15].
However, for finite , it is known that there are important differences between the classical and quantum behavior [16][17][18]. The wavefunction is not as localized as the classical PDF and, due to tunnelling, can spread into regions that are classically forbidden. Thus, the quantized version provides a natural coarse-graining of classical phase space into units of Planck's constant, h, often included as a normalization constant for entropy in classical statistical mechanics. Moreover, different components of the wavefunction carry different complex phase factors which leads to interference effects that are not present in the classical setting. Yet another important qualitative difference is the emergence of dynamical Anderson localization [19,20] and many-body localization [21][22][23], which prevents the wavefunction from sampling the entire classically chaotic region accessible to the classical trajectories. Finally, while this approach is straightforward for Hamiltonian systems, it is not trivial to determine a natural quantization procedure for general non-Hamiltonian dynamical systems.

D. Koopman-von Neumann Approach
Nonlinear classical phase space dynamics can be faithfully embedded within a quantum mechanical systemeven for equations of motion that are not Hamiltonian. Just after the birth of quantum mechanics, Koopman [24] and von Neumann [25,26] realized that classical mechanics can be formulated on Hilbert space in a manner that is exactly analogous to quantum mechanics. The classical Liouville equation, which expresses the conservation of probability on phase space, and its space-time adjoint, which expresses the evolution of a conserved quantity, can be recast as an equivalent Koopman-von Neumann (KvN) Schrödinger equation [20,27]. The KvN Hamiltonian (Eq. 20) associated with the KvN equation (Eq. 17) is linear in the momentum and results from the quantization of an associated constrained Hamiltonian system [28,29] on twice the classical phase space dimension, where the canonical momenta represent Lagrange multipliers that enforce the classical equations of motions as constraints [30]. Heisenberg's uncertainty principle applies to each of the original variables and its conjugate momentum, i.e. the Lagrange multiplier, but does not affect any pair of the original variables. Thus, there is complete fidelity to the classical phase space evolution.
Following the KvN approach in Sec. II leads to the simple derivation of the generalized KvN equation (Eq. 17), which applies to arbitrary classical dynamical systems, starting from the assumption that the probability distribution, f , is the inner product of a complex probability amplitude ψ with its adjoint ψ † , (Eq. 12).
To this author's knowledge, Chirikov, Izrailev, and Shepelyanski [20] were the first to publish the generalized form of the KvN equation, which applies to arbitrary classical dynamical systems, which they attribute to an earlier preprint by Gorban and Okhonin [27]. These authors clearly appreciated the meaning of the KvN equation as a first order Hermitian PDE that exactly preserves the original phase space dynamics of the underlying classical system, and, as one which would display the features of classical chaos rather than of quantum chaos.
Chruściński [28,29] referred to the KvN Hamiltonian as the "Quantum mechanics of damped systems," because a damped system can be embedded within the constrained Hamiltonian discussed in Sec. II E. While this is a valid quantization of the constrained Hamiltonian, it does not have the physical meaning of being a "damped" quantum mechanical system, which would typically be treated using the Master Equation. Rather, it is the KvN Hamiltonian for an arbitrary classical system of ODEs whether damped or not.
Understanding the properties of the unitary KvN evolution operator, i.e. propagator or transfer operator, corresponding to the Hermitian KvN Hamiltonian should be quite interesting from the point of view of of dynamical systems theory. The closely related Perron-Frobenius [31] operator and Koopman operator [32,33] have been studied and used extensively to characterize invariant measures of and to develop reduced order models for nonlinear dynamical systems. For divergence-free, i.e. measurepreserving, flows, these three operators are equivalent. The recognition that the KvN Hamiltonian for dissipative systems still leads to unitary evolution has already found use in the development of high-order conservative numerical discretizations of the advection operator for fluid dynamics [34,35] and plasma physics [36].
The KvN approach can be applied to any system of differential equations that are explicitly presented as first order in time. For example, the method of lines can be used to generate a numerical discretization of a PDE as a finite set of ODEs that can be treated in an equivalent fashion. Hence, the KvN approach can also be applied to important PDEs in classical physics such as the Maxwell's equations in conducting media, the Navier-Stokes equation, the collisional evolution of the kinetic PDF, and general N -body problems in gravitation, plasma physics, and molecular dynamics. Although the probabilistic description of a PDE might appear to be expensive, this is the necessary setting for understanding the evolution of the PDF over regions of phase space and provides the essential framework for uncertainty quantification. It is also necessary for describing interactions with random processes and stochastic forcing terms which are often used to model effects due "noise" and turbulence. In fact, these examples are simply the probabilistic classical field theory and statistical field theory analogs of quantum field theory.

E. Semiclassical Evolution
As far as the classical system is concerned, the dynamics of the complex phase factor (Eq. 13) associated with the wavefunction is simply a choice of gauge. In fact, there are alternate action principles that can be used to generate the Liouville equation [37,38] and these formulations can also be used to define the dynamics of the phase factor of an associated KvN system.
Recently, Klein [39] suggested a form for the dynamics of the complex phase factor (Eq. 16) that approximates semiclassical dynamics of a Hamiltonian system and Bondar, Gay-Balmaz, and Tronci, [40], proposed that it is, in fact, the semiclassical dynamics that is physically relevant. These authors showed that the semiclassical KvN approach could be used to couple a semiclassical system to a quantum system in a self-consistent manner. Note, however, that both Refs. [39,40] neglected to discuss the importance of the Maslov index [41][42][43][44] for obtaining the correct semiclassical dynamics.
The focus in this work is on quantum simulation of an arbitrary classical system of ODEs, and, developing a self-consistent framework for the semiclassical dynamics of Hamiltonian systems of ODEs will be left for future investigations.

F. Quantum Simulation
A quantum computer with finite resources can be used to simulate a finite-dimensional approximation of the KvN Hamiltonian operator. One must discretize the equations in a manner that can be represented with a finite number of qubits, which leads to a quantum mechanical coarse-graining and regularization of phase space that is rather different than that due to dissipative effects such as particle collisions. In fact, one should squeeze the initial and final measurement states to achieve the uncertainty limit set by numerical discretization.
If the KvN Hamiltonian is sparse, e.g. local or banded, then the quantum representation of the classical system leads to exponential savings in the memory and computational work [2] required for an Eulerian discretization of the Liouville equation. Classical time-dependent Monte Carlo algorithms [10] can also provide a similar savings. Since quantum algorithms based on quantum walks can often provide a quadratic speedup over classical Monte Carlo algorithms [14,15], it is of interest to determine whether similar gains would generally be expected for classical dynamics.

G. Outline
The next section describes the derivation of the generalized Koopman-von Neumann representation of classical mechanics. Section III discusses the important special case of Hamiltonian dynamics, including canonical Hamiltonian systems, generalized Hamiltonian systems, and general variational systems. Section IV relates the KvN Hamiltonian to an action principle for the wavefunction and uses Noether's theorem to derive a number of conservation laws for the system. Section V discusses the steps necessary to obtain high-fidelity quantum simulation of a classical system. Finally, estimates of the complexity of using the KvN approach to classical dynamics are discussed in Sec. V C. A summary of the conclusions is presented in the final section.
The Einstein summation convention is used throughout.

A. Classical Dynamics on Phase Space
Consider the solution of a system of d classical ODEs with coordinates, where v(x, t) is an arbitrary vector field. There are two different ways to interpret the dynamics of observables or functions on phase space, as illustrated in Fig. 1. From the Lagrangian point of view, one can start with a number of initial conditions and follow the dynamics forward in time. In this picture, observables evolve via the "total" time derivative d/dt. From the Eulerian point of view, one can think of phase space as a manifold with coordinates given by x. In this picture, observables evolve in time through partial differential equations (PDEs) on phase space. These two pictures are related through the equivalence between the total time derivative and the advection operator Any constant of the motion, ϕ(x, t), is simply advected along the flow by the advection operatoṙ The evolution operator, also known as the propagator or transfer operator, corresponding to this PDE is known as the Koopman operator [31]. The expectation value of an observable, O(x, t), is given by integration over the probability distribution function (PDF), f (x, t), is advected by the flow v(x, t), as described by the Liouville equation (Eq. 8).
Thus, as explained in App. A, the PDF is a volume form or phase space density. Conservation of phase space density is described by the Liouville equatioṅ The propagator corresponding to this PDE is known as the Perron-Frobenius operator [31]. The advection of the PDF by the flow is illustrated schematically in Fig.  1. The meaning of the Liouville equation is that the PDF, f (x, t), is an invariant measure of the velocity field, v(x, t), over space-time. Using integration by parts over space-time demonstrates that the the Liouville operator that appears in (Eq. 8) is the anti-Hermitian adjoint of the advection operator Iff the velocity is divergence free, then the advection operator (Eq. 5) and the Liouville operator (Eq. 10) are anti-Hermitian operators, i.e. A = −Â † , and hence, identical, In the divergence-free case, the fact that these operators are anti-Hermitian implies that the corresponding evolution operator is unitary. Discussion of the general form of the Liouville equation on space-time and the corresponding conventions for the divergence operator are discussed in Appendices A-C.

B. Koopman-von Neumann Hamiltonian
The Koopman-von Neumann (KvN) approach to classical mechanics parallels the development of the postu-lates of quantum mechanics. The probability distribution function, f , is defined as the inner product of the complex probability amplitude, ψ(x, t), with its adjoint, ψ † (x, t), Thus, the probability amplitude or wave function is closely related to notion of the "square root" of f . The expectation value of any phase-space observable, O(x, t), is given by integration over the probability distribution function (PDF), f (x, t), This requires that the PDF is normalized to yield unit probability after integration over all of phase space. Assume that the phase, ϕ, satisfies the general equation of motion [39] Then, taking the time derivative of Eq. 13 and using Eqs. 8 and 16 proves that ψ satisfies the equatioṅ Simply multiplying this equation by i yields the generalized Koopman-von Neumann equation The KvN equation has the form of a Schrödinger equation As can be proven using integration by parts, the KvN Hamiltonian operator,Ĥ, is Hermitian over the Hilbert space inner product defined on any two functions of phase space via Hence, define the usual momentum operator,P = −i ∇, and position operator,x = x, and promote any function of the coordinates to an operator, e.g.v = v(x, t) and W = W (x, t), via the formal Taylor series expansion. Thus, one arrives at the generalized Koopman-von Neumann Hamiltonian operator Because the KvN Hamiltonian is Hermitian, the corresponding KvN evolution operator is unitary. Koopman and von Neumann only considered the case of Hamiltonian dynamics, but the generalized form leads to a unitary evolution operator for any set of ODEs, even if they are not Hamiltonian. In contrast to the usual Schrödinger equation, the KvN Hamiltonian is linear in momentum,P, because the KvN equation is a first order PDE rather than a second order PDE. Hence, there are a number of important differences between the behavior of the KvN equation and the usual Schrödinger equation as well as in the mathematical analysis of the two types of equations.
The linear dependence on momentum implies that Heisenberg's equations of motion for the original phase space variables are exactly the same as the original classical equations of motion, Thus, thex operators have the same solution as the original classical equations of motion. In contrast, Heisenberg's equations of motion for the conjugate momenta generically receive a quantum correction when compared to the classical limit (Eq. 37).
If there is a natural volume form on phase space, then the PDF and the wavefunction can also be treated as scalar fields. The Koopman-von Neumann approach for scalar fields is discussed in App. B. The Koopman-von Neumann equation on space-time is discussed in App. C.

C. Koopman-von Neumann Evolution Operator
The phase space evolution described by Eqs. 8, 16, and 17 can be solved using the method of characteristics. The solution of the ODE's in Eq. 4 has the form x = ξ(x 0 , t) for the initial conditions x(t 0 ) = x 0 . The inverse relation x 0 = ξ −1 (x, t) allows one to express the solution of these equations in terms of the determinant The PDF is determined by and, if one defines then the phase is determined by Thus, the evolution of the amplitude, ψ, can be determined from the definition in Eq. 13 The choice of the complex phase of J 1/2 0 is irrelevant for the classical system; for example, one could formally use |J 0 | 1/2 instead. However, the phase shift due to the square root of the Jacobian, ∆ϕ = −νπ/2, is given by the integer Maslov index, ν, which counts the number of zeros of the Jacobian along the trajectory since the starting point [31,[41][42][43][44]. This phase shift must be correctly accounted for in order to obtain the correct semiclassical phase factor.
The KvN evolution operator,Û, is defined to satisfy the operator equation Using the general solution given above, the KvN evolution operator can be written as Again, because the generalized KvN Hamiltonian (Eq. 20) is Hermitian, the KvN evolution operator is unitary.

D. Semiclassical Dynamics and the Phase Factor
The dynamics of the phase factor,φ = −W/ , has no effect on the classical dynamics. Within the confines of classical dynamics, the phase, ϕ, is not measurable, and, hence, the choice of W is equivalent to a choice of gauge.
On the other hand, the semiclassical approximation to quantum dynamics represents the propagator as a sum over classical paths [9,31,44], and, in this case, the phase factor is important for describing interference between paths. Similarly, when the classical system is coupled to a quantum system, then the dynamics is also sensitive to the semiclassical phase factor introduced by the semiclassical system [40]. In other words, Schrödinger's cat has rather different semi-classical phase factors depending on whether it is dead or alive, and, because the cat is entangled with the quantum part of the system, this can lead to measurable interference effects due to the semiclassical phase of the cat itself. If the classical system is itself a Hamiltonian system, with the Hamiltonian, H(q j , p j , t), specified as a function of generalized coordinates, q j , and conjugate momenta, p j , then there is a natural choice of phase factor [39,40]. If the phase velocity is determined by the classical Lagrangian, L, then the phase factor is equal to the classical action and corresponds to the semiclassical phase. The semiclassical phase factor is related to the momentum through the Eikonal approximation, and satisfies the Hamilton-Jacobi equation, Here, the partial derivatives are taken as constant with respect to the initial value of the coordinates, q j 0 . Taking the total time-derivative of ϕ yields Eq. 31.
For short time intervals, only a single classical path will generically contribute to the propagator. However, for longer time intervals, multiple classical paths can contribute, each with a separate phase. In order to obtain the semi-classical phase factor for each branch, one must take care to ensure that the Maslov index [41][42][43] is correctly accounted for, by choosing the correct sign of the square root of the Jacobian, J 1/2 0 , in Eq. 28. This generates an additional −π/2 phase jump whenever J 0 passes through zero and leads to an additional overall phase shift, ∆ϕ = −νπ/2, where ν counts the number of times that J 0 passes through zero along the trajectory [31,44]. (This important fact appears to have been neglected in the discussion of Refs. [39,40].) It is also possible to give the amplitude additional index structure, ψ ijk... , so that the various components transform as a representation of the Poincaré group under changes of reference frame. Clearly, this extension is important for describing the semi-classical evolution of higher spin fields as well as for describing interactions with gauge fields. While this does not affect the classical dynamics, the semiclassical dynamics could certainly be an interesting avenue to pursue in future work.
In the more general non-Hamiltonian case considered here, a natural choice ofφ is not immediately obvious. The trivial choice, W = 0, corresponds to the constrained classical action (Eq. 36) considered in the next subsection.
The general case is also related to the quantum mechanics of a charged particle in an N -dimensional vector potential, A(x, t), and scalar potential, Φ(x, t), in the limit that the canonical momentum is much smaller than the vector potential, |P | |qA|, where q is the charge of the particle. In this limit, one can identify the relations, where m is the mass of particle. Thus, in this limit, the choice of φ = −W is equivalent to a choice of scalar potential, Φ.

E. Constrained Hamiltonian
In the limit → 0, the KvN Hamiltonian of Eq. 20 becomes the classical Hamiltonian Because it is linear in the momentum, P, this corresponds to a constrained Hamiltonian system. Any set of N classical ODE's can be generated by using the action principle corresponding to this Hamiltonian [30] S[x, P; t] : which can be interpreted as as a sum over constraints. Variation with respect to the Lagrange multipliers, P j , enforces the classical equation of motion for the coordinates, x j . Variation with respect to the coordinates, generates the classical equation of motion for the Lagrange multipliersṖ = −(∇v) · P − ∇W.
The equations of motion ensure that the coordinate transformation between different points in time always remains symplectic, and, in fact, canonical. The constrained Hamiltonian can be quantized via any acceptable quantization procedure that reduces to Eq. 35 in the limit → 0. The symmetric Weyl quantization rule leads toĤ in Eq. 20. The constrained Hamiltonian in Eq. 35 implies that the quantized variablesx i commute with one another. In contrast, other types of constraints, e.g. holonomic constraints, would require modifying the Poisson bracket to the Dirac bracket [45]. This would lead to nontrivial commutators and nontrivial uncertainty relations which could potentially cause deviations from the exact classical dynamics.

F. General Examples
Consider a single autonomous ODE in a single variable, The KvN Hamiltonian operator is simplŷ The equation of motion can be derived from the constrained Hamiltonian, H = P v. In addition to Eq. 38, the canonical momentum satisfies the equatioṅ The classical solution to this equation is simply Since these equations are symplectic, the phase space area is conserved, as can be checked directly from the equations of motion. Heisenberg's equation of motion is A specific instance of this case is given by v = γx, where the sign of γ determines whether the motion exhbits exponential growth or damping. The constrained Hamiltonian is simply H = P γx and the classical solutions on extended phase space are The corresponding KvN Hamiltonian iŝ Consider the general linear equation seṫ The constrained Hamiltonian leads to the additional equations of motion, In this case, the equations of motion have a simple anti-Hermitian block off-diagonal form. The corresponding KvN Hamiltonian is given bŷ Since the last two examples had linear equations of motion resulting from Hamiltonians that are quadratic forms, Heisenberg's equations of motion are the same as the classical equations.

A. Canonical Hamiltonian Systems
Hamiltonian systems of differential equations are of fundamental importance to physics and mathematics. Canonical coordinates are defined by pairs of configuration space coordinates, q j , and the corresponding conjugate momenta, p j . In canonical coordinates, Hamilton's equations take the canonical forṁ Hamilton's equations are divergence-free in canonical coordinates, which implies that the natural measure is the canonical measure d d qd d p. The expectation value of any phase-space observable, O(q j , p k , t), is given by integration over the probability distribution function (PDF), Similarly, the inner product of any two functions is defined via Since Hamilton's equations are divergence-free, the Liouville equationḟ is equivalent to it's space-time adjoint, advection by the flow. Here, the canonical Poisson bracket is defined as Let us now introduce a new convention for the Lagrange multipliers P j , Q k so that they are conjugate to the q j , p k in the "canonical" manner The constrained Hamiltonian takes the form The equations of motion of the Lagrange multipliers arė The usual Dirac quantization procedure for Hamiltonian systems promotes the classical Poisson bracket to commutation relations. Instead, the KvN approach only promotes the Poisson brackets of the Lagrange multipliers (Eq. 55) to commutation relationŝ Since the canonical equations of motion are divergencefree, the KvN Hamiltonian is simplŷ In this formulation, the Heisenberg uncertainty principle only applies to the pairs q j , P j and p j , Q j , but not to the pairs of q j , p j .

B. Canonical Hamiltonian Examples
Consider the classical harmonic oscillator, with Hamiltonian The dynamics can be exactly represented by the constrained Hamiltonian The additional nontrivial equations of motion for the Lagrange multipliers arė Hence, the dynamics on the extended phase space consists of two harmonic oscillators which rotate in the same direction. The quantum Hamiltonian iŝ where tan(θ) = p/q. Consider the dynamics of an integrable system of nonlinear oscillators defined by the Hamiltonian H 0 (J j ), where J j are the conserved action variables corresponding to each oscillator. The evolution of the conjugate phases, θ j , is determined by the frequencies ω j 0 := ∂ Jj H 0 . If one defines the number operators aŝ then, this leads to the KvN Hamiltonian While the number operators are conserved like the action coordinates, the angle operators defined bŷ satisfy the Heisenberg equations of motion

C. Generalized Hamiltonian Systems
In general, a Hamiltonian system of equations is defined both by the Hamiltonian H(x, t) and by the Poisson bracket [46,47] {a, b} = J jk ∂ j a∂ k b.
The Poisson 2-vector, J jk (x, t), must be an antisymmetric tensor and must satisfy the Jacobi identity 0 = {{a, b} , c} + {{b, c} , a} + {{c, a} , b} (70) which yields the relation The classical equations of motion generated by the Hamiltonian H(x, t) arė Following the general procedure outlined previously, the corresponding constrained Hamiltonian is This leads to the following dynamics of the Lagrange multipliers:Ṗ The KvN Hamiltonian for generalized Hamiltonian systems isĤ Once again, this form follows the convention for general systems of differential equations, but does not correspond to the convention used for canonical Hamiltonian systems used in Sec. III A. When the generalized Hamilton's equations of motion are non-singular, they have a canonical definition of the volume form. Thus, as explained in App. E, one can also derive a corresponding KvN equation where the PDF and the wavefunction are treated as scalar fields.

D. General Variational Systems
A general variational system of equations is defined by extremizing the action principle [47] S[ The action is defined in terms of the Poincaré 1-form, α = α j dx j . (For an introduction differential forms and exterior calculus, see Refs. [48,49].) This leads to the equations of motion where the symplectic 2-form, Ω := dα, which defines the Lagrange bracket, is the exterior derivative of the Poincaré 1-form; i.e. in components These equations can be written in the form of Hamilton's equations of motioṅ with the identification of the Poisson tensor as the inverse of the Lagrange form, J = Ω −1 . The Jacobi identity follows from the fact that the symplectic 2-form is closed, since it is exact dΩ = d 2 α = 0. The canonical equations of motion result from the Poincare 1-form, α = p j dq j . This yields the canonical symplectic 2-form Ω = dp j ∧ dq j (82) and the canonical Poisson 2-vector where, now, the X j are the Lagrange multipliers that enforce the variational equations of motion. The action principle is now The Hamiltonian form in Eq. 73 and the variational form in Eq. 84 are equivalent with the definition Due the fact that the Poisson tensor is conserved by the equations of motion, dΩ/dt = 0 (proven in App. D), the Lagrange multipliers, X j , satisfy the equations of motioṅ The analogous quantized operators,X j , should be given a Hermitian form Thus, the KvN Hamiltonian for general variational systems isĤ = 1 2 P jv j +v jP j +Ŵ (89) The final form corresponds to the convention that was used for canonical Hamiltonian systems in Sec. III A. For general variational systems, there is a canonical volume form induced by the symplectic 2-form. Thus, as explained in App. E, one can also derive a corresponding KvN equation where the PDF and the wavefunction are treated as scalar fields.

A. Action Principle for General Systems
The Lagrangian that generates the KvN equation for the wave function, Eq. 16, as well as the Liouville equation and its adjoint, Eqs. 6 and 8, is given by where integration by parts is used to derive the final expression. The Hamiltonian that generates these partial differential equations is then given by where integration by parts is used to derive the final expression. The definition of the Hamiltonian also allows one to derive the expression There are a number of additional action principles that have been discussed in the literature [37,38,40]. After conversion to Koopman-von Neumann form, these action principles can be interpreted as providing an alternate dynamics of the phase factor. Hence, it may be interesting to explore the KvN dynamics that is implied using these alternate formulations as well.

B. Action Principle for Canonical Hamiltonian Systems
For the case of a canonical Hamiltonian system, with Hamiltonian, H(q j , p j , t), integration by parts allows one to express the Hamiltonian in the form With the gauge choice φ = −W = 0, the probability density, f H , that determines the KvN Hamiltonian density, is This will only be equal to f , iff the phase ϕ is canonically conjugate to ln (f ) Due to the Jacobi identity, this condition is preserved by the dynamics, so that if Eq. 100 is true as an initial condition, then it is true for all time.
In the canonical case, the Poisson bracket can also be expressed as a divergence, so that and this form of f H as a divergence is also correct for generalized Hamiltonian dynamics. Ref. [40] claimed that Eq. 100 is not compatible with the fact that f H can be written as a divergence, Eq. 101, and, thus, as a surface integral, which they assumed can be made to vanish at the boundary. However, in order to satisfy Eq. 100, the phase function, ϕ, must be multiple-valued, so that the surface integral does not vanish. Alternatively, assume that the semiclassical phase factor is chosen, so that it is equal to the classical action [39,40] In this case, using integration by parts, the difference between f H and f is determined to be [40] f Because the semi-classical phase factor also satisfies the Hamilton-Jacobi equation the difference, f H − f , vanishes identically and the final expressions for the Hamiltonian in Eq. 95 and Eq. 96, precisely yield the classical energy. This condition can be used to simplify a number of the results of Ref. [40].

C. Symmetries and Conservation Laws
Noether's theorem implies that a symmetry of the action leads to a conservation law for a corresponding density, Q. In this case, there are a number of conservation laws in the form of a space-time divergence The conservation of probability density, f , also implies that, if the system results from a classical Hamiltonian, H(q j , p j ), that is independent of time, ∂ t H = 0, then the classical energy density, Hf , is conserved. Similarly, if the classical Hamiltonian is independent of one of the coordinates, q j , so that ∂ q j H = 0, then the classical momentum density, p j f is conserved. If W is chosen to be a function of conserved quantities, e.g. H and/or p j , etc., then the classical Hamiltonian conservation laws imply that the KvN counterparts also hold true.

A. Heisenberg Uncertainty
The only nontrivial commutators in the extended phase space are Hence, the only nontrivial Heisenberg uncertainty relations are where the uncertainty in quantity A is defined via Thus, it is possible to make a simultaneous measurement of all of the classical phase space variables, x j , or all of the Lagrange multipliers, P j .

B. Numerical Discretization
In order to represent the classical phase space dynamics with a finite number of qubits, one must construct a finite-dimensional numerical approximation of the KvN Hamiltonian. Assume that each coordinate, x j , is periodic on the length, X j max , and is represented with an integer number of levels, L j , so that the level spacing is ∆x j = X j max /L j . Then, the Fourier representation of the conjugate momenta, implies that these coordinates are also periodic and have L j levels with level spacing ∆P j = h/X j max and the range P j,max = hL j /X j max . Thus, as illustrated in Fig. 3, the phase space uncertainty due to the discreteness of the representation, 3: A numerical discretization of phase space implies finite numerical widths ∆x and ∆P . However, the phase space area ∆x∆P = h/L due to the discretization is much smaller than than the Heisenberg limit allows. Squeezed states can be used to reduce the uncertainty to the numerical discretization limit, so that σx = ∆x. ∆x j ∆P j = h/L j , is much smaller than the Heisenberg limit allows.
Often, one considers coherent states to be the analog of classical states. However, a coherent state, which saturates the Heisenberg bound, will have a width that is large compared to the classical level spacing, i.e. σ x j /∆x j = (L j /4π) 1/2 and σ Pj /∆P j = (L j /4π) 1/2 . Since one is only interested in the dynamics of the original phase space, one can use squeezed states to reduce the uncertainty in the x j coordinates of interest and increase the uncertainty in the P j coordinates. Squeezing the uncertainty in x j by the factor of L −1/2 j reduces the quantum uncertainty to the limit set by numerical discretization, σ x j = ∆x j . This increases the uncertainty in P j by the factor L 1/2 j so that σ Pj = L j /X j max = (L j /4π)∆P j . This implies that the relative uncertainty satisfies σ x j /X j max = 1/L j while σ Pj /P j,max = 1/4π. Thus, one can use squeezed states as initial conditions and as final states in order to perform measurements that saturate the uncertainty limit set by numerical discretization.

C. Complexity Estimates
It is already known that there are quantum algorithms that can speed up the solution of a linear system of ordinary differential equations [6] and partial differential equations, such as wave equations [7]. Can this be extended to nonlinear systems?
The KvN representation of classical dynamics implies that Lloyd's proof [2] that a quantum computer can be used to accelerate quantum simulation can also be applied to classical dynamics, even for non-Hamiltonian systems. If the quantum Hamiltonian in Eq. 20 is sufficiently sparse, e.g. because it is local, then one can break the operator into a small number, m, of non-commuting parts,Ĥ = m j=1Ĥ j , where m is independent of N . Application of the Trotter-Suzuki product formula can be used to generate an efficient simulation of the unitary evolution with error proportional to j =k Ĥ j ,Ĥ k .
Recent results for general s-sparse Hamiltonian simulation [50] using n qubits imply that the quantum simulation complexity is proportional to sn Ĥ t max , neglecting polylogarithmic factors. This implies large gains over an Eulerian discretization of the Liouville equation.
For example, consider the case of M particles that experience local interactions between them, traveling in 2d phase space dimensions, including both configuration space and momentum (or velocity) space, so that the total phase space dimension is D = 2dM . If an Eulerian discretization is employed on a grid of L = 2 points in each direction, then this requires L D = 2 D degrees of freedom, and D qubits to represent the PDF. If the kinetic energy is quadratic in velocity, then the maximum kinetic energy scales as dL 2 , so that Thus, as compared to a classical Eulerian discretization of the Liouville equation, this implies an exponential speedup in D and a polynomial speedup in L. If the number of dimensions, D, is large, then the degree of the polynomial speedup, D/2, is also large. Similar estimates can also be obtained for the max norm of the commutator. Note, however, that some important calculations can require a large number of time steps [13], potentially scaling as a power of D, and this would reduce the expected savings to polynomial at best. However, the more interesting comparison is between the quantum simulation and the best probabilistic classical algorithm. For high dimensional PDEs, the best probabilistic classical algorithms that are known are typically some form of generalized time-dependent Monte Carlo (MC) algorithm, broadly including particle-based techniques such as particle-in-cell (PIC) and molecular dynamics within the framework of Markov Chain Monte Carlo (MCMC). With appropriate assumptions, classical MC algorithms can also provide a large speedup over Eulerian discretization [10] because the accuracy of the results will scale as the square root of the number of samples, i.e. the computational work scales as ∼ ε −2 , where ε is the required accuracy. However, for quantum algorithms based on quantum walks, the accuracy generically decreases as the number of samples itself [14,15], which leads to a quadratic improvement in computational work, ∼ ε −1 , for a given accuracy. Thus, it is possible that quantum simulation using the KvN approach to classical dynamics could lead to a quadratic improvement in com-putational cost relative to classical time-dependent MC algorithms.
In the typical case of quantum simulation and in the typical case of linear evolution, one can often argue that if the state of interest is smooth, e.g. because it is a ground state or a low-lying excited state, then the max norm may be a severe overestimate of the errors involved in computing the measured expectation values. If the solution is well-resolved, then Fourier harmonics of the probability distribution should decay exponentially with increasing L. This can potentially reduce the actual expectation value of the error to become exponentially small and can potentially imply an exponential speedup [51]. However, for nonlinear classical dynamics, the PDF generically develops finely mixed phase space structure that approaches a fractal distribution in the infinite time limit [52]. This implies that this kind of argument is not likely to apply in the setting of nonlinear classical dynamics.
There are also additional restrictions on the ability of quantum algorithms to speed up the calculation. A general quantum simulation program must have initialization, simulation, and output (measurement) subprograms. Care must be taken to ensure that the complexity of the initialization and output stages are less than or equal to the simulation stage. This implies that it is not desirable to set each initial condition individually or to measure the entire PDF over all states. The input states should be relatively easy to construct so that the initialization step is "sparse." Similarly, the output should consist of a relatively small number of measurements, e.g. low-order moments of PDF, so that the measurement step is also "sparse." Finally, note that, if the underlying system is Hamiltonian and simulating the quantized Hamiltonian system is sufficient for the intended calculation, then, because the KvN approach leads to a system with twice the phase space dimension, simulating the quantized Hamiltonian is the more efficient computational approach. If the computational gains versus classical Monte Carlo are only quadratic, then doubling the phase space dimension would effectively eliminate the advantage. Moreover, while quantizing the Hamiltonian clearly yields selfconsistent semiclassical dynamics, the standard semiclassical approximation to the evolution operator [31] is not guaranteed to be unitary. Thus, a consistent approach to semiclassical evolution using the KvN framework has yet to be developed.

VI. CONCLUSION
In conclusion, quantum computers can be used to simulate nonlinear non-Hamiltonian classical dynamics on phase space by using the generalized Koopman-von Neumann formulation of classical mechanics. The Koopmanvon Neumann formulation implies that the classical phase space dynamics expressed by the Liouville equa-tion can be recast as an equivalent Schrödinger equation for the wavefunction on Hilbert space. The wavefunction completely specifies the probability distribution function and its dynamics is generated by a Hermitian Hamiltonian operator and a unitary evolution operator. Thus, a quantum computer with finite resources can be used to simulate a finite-dimensional approximation of this unitary evolution operator.
The conservation of probability on phase space can be expressed as an equivalent Schrödinger equation that is linear in momenta. The equivalent Schrödinger equation corresponds to the quantization of a constrained Hamiltonian system with twice the dimension of the original phase space, where the conjugate momenta act as Lagrange multipliers that enforce the equations of motion on the original phase space. Heisenberg's uncertainty principle applies to each variable and the corresponding Lagrange multiplier, but not between any of the variables of the original phase space. Squeezing the uncertainty of the original phase space variables and enhancing the uncertainty of the conjugate Lagrange multipliers allows the quantum uncertainty to be reduced to the limit set by numerical discretization. Hence, there is complete fidelity to the classical phase space dynamics.
Quantum simulation of classical dynamics is exponentially more efficient than an Eulerian discretization of the Liouville equation if the Koopman-von Neumann Hamiltonian is sparse. It is possible that utilizing quantum walks and associated techniques could also lead to a quadratic improvement over classical time-dependent Monte Carlo methods. Exploring the advantages and disadvantages of the Koopman-von Neumann representation for specific examples of classical dynamical systems and developing quantum simulation algorithms for an accurate approximation to the semiclassical evolution operator are important directions for future work.
The final term vanishes due the fact that the symplectic form is closed, dΩ = d 2 α = 0. The first two terms cancel due to the equations of motion.
The time derivative of the Poincaré form is Due to the equations of motion, the time derivative is the differential of the Lagrangian Since d 2 L = 0, this yields an alternate derivation of the conservation of the symplectic form. Conservation of the symplectic 2-form also implies conservation of the Poisson 2-vector, due to the definition, J = Ω −1 .