Quantum Algorithm for the Vlasov Equation

TL;DR

This paper introduces a quantum algorithm for simulating the Vlasov-Maxwell system in plasma physics, promising exponential speedups over classical methods for high-dimensional problems, with discussions on extensions and current limitations.

Contribution

The paper develops a quantum algorithm for the Vlasov-Maxwell system, specifically for electrostatic Landau damping, demonstrating potential exponential speedup in high-dimensional plasma simulations.

Findings

Quantum algorithm scales as polylogarithm of velocity grid points

Classical costs scale linearly with grid points and time

Quantum approach faces measurement error challenges

Abstract

The Vlasov-Maxwell system of equations, which describes classical plasma physics, is extremely challenging to solve, even by numerical simulation on powerful computers. By linearizing and assuming a Maxwellian background distribution function, we convert the Vlasov-Maxwell system into a Hamiltonian simulation problem. Then for the limiting case of electrostatic Landau damping, we design and verify a quantum algorithm, appropriate for a future error-corrected universal quantum computer. While the classical simulation has costs that scale as $\mathcal{O}(N_v t)$ for a velocity grid with $N_v$ grid points and simulation time $t$, our quantum algorithm scales as $\mathcal{O}(\text{polylog}(N_v) t/\delta)$ where $\delta$ is the measurement error, and weaker scalings have been dropped. Extensions, including electromagnetics and higher dimensions, are discussed. A quantum computer could efficiently handle a high-resolution, six-dimensional phase-space grid, but the $1/\delta$ cost factor to extract an accurate result remains a difficulty. This paper provides insight into the possibility of someday achieving efficient plasma simulation on a quantum computer.

Equations

$$\frac{\partial f_s}{\partial t}+\mathbf{v}\cdot \nabla f_s+\frac{q}{m}\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right)\cdot \frac{\partial f_s}{\partial \mathbf{v}}=0,$$

$$\frac{\partial f_1}{\partial t}+\mathbf{v}\cdot \nabla f_1=\frac{e}{m}\left(-{\mathbf{E}}_1\cdot \mathbf{v}{f}_M+(\mathbf{v}\times {\mathbf{B}}_0)\cdot \frac{\partial f_1}{\partial \mathbf{v}}\right).$$

$$\hat{\mathbf{k}}$$

$$\hat{\mathbf{v}}$$

$$\hat{\mathbf{E}}$$

$$\frac{\partial \tilde{f}}{\partial t}=-i\mathbf{k}\cdot \mathbf{v}\tilde{f}-\tilde{\mathbf{E}}\cdot \mathbf{v}{f}_M+\frac{1}{c_n}(\mathbf{v}\times {\mathbf{B}}_0)\cdot \frac{\partial \tilde{f}}{\partial \mathbf{v}},$$

$$\frac{\partial \tilde{\mathbf{E}}}{\partial t}=i{c}_n\mathbf{k}\times \tilde{\mathbf{B}}+\int \mathbf{v}\tilde{f}{d}^{3}v,$$

$$\frac{\partial \tilde{\mathbf{B}}}{\partial t}=-i{c}_n\mathbf{k}\times \tilde{\mathbf{E}}.$$

$${\tilde{f}}^{\prime }(\mathbf{v}):=i\sqrt{\frac{\xi }{f_M(\mathbf{v})}}\tilde{f}(\mathbf{v}),$$

$$\frac{\partial {\tilde{f}}^{\prime }}{\partial t}=-i\mathbf{k}\cdot \mathbf{v}{\tilde{f}}^{\prime }-i\mu \tilde{\mathbf{E}}\cdot \mathbf{v}+\frac{1}{c_n}(\mathbf{v}\times {\mathbf{B}}_0)\cdot \frac{\partial {\tilde{f}}^{\prime }}{\partial \mathbf{v}},$$

$$\frac{\partial \tilde{\mathbf{E}}}{\partial t}=i{c}_n\mathbf{k}\times \tilde{\mathbf{B}}-i\int \mathbf{v}\mu {\tilde{f}}^{\prime }\frac{d^{3}v}{\xi },$$

$${\hat{H}}_v=\mathbf{k}\cdot \mathbf{v}+\frac{i}{c_n}(\mathbf{v}\times {\mathbf{B}}_0)\cdot \frac{\partial }{\partial \mathbf{v}}.$$

$$\int g(\mathbf{v}){\hat{H}}_{v}h(\mathbf{v})d^{3}v=\int h(\mathbf{v}){\hat{H}}_v^{\ast }g(\mathbf{v})d^{3}v$$

$$\frac{\partial |EM⟩}{\partial t}=-i{\hat{H}}_l|EM⟩=-i{c}_{n}k\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& -1& 0\\ 0& -1& 0& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{E}_x\\ E_y\\ B_x\\ B_y\end{array}\right].$$

$$\frac{\partial {\tilde{f}}^{\prime }(\mathbf{v})}{\partial t}$$

$${\tilde{F}}_{j_z}$$

$${\tilde{F}}_{j_z}^{\prime }:=i\sqrt{\frac{\Delta v}{G_{j_z}}}{\tilde{F}}_{j_z},$$

$$\frac{\partial {\tilde{F}}_j^{\prime }}{\partial t}=-ik{v}_j{\tilde{F}}_j^{\prime }-i{\alpha }_j\tilde{E}{v}_j,$$

$$\frac{\partial \tilde{E}}{\partial t}=-i\sum _j{\alpha }_j{\tilde{F}}_j^{\prime }{v}_j,$$

$$|x⟩=\eta \left(\sum _{j=0}^{N_v-1}{\tilde{F}}_j^{\prime }|j⟩+\tilde{E}|N_v⟩\right),$$

$$\hat{H}=\sum _{j=0}^{N_v-1}{v}_j\left[k|j⟩⟨j|+{\alpha }_j\left(|j⟩⟨N_v|+|N_v⟩⟨j|\right)\right].$$

$$|x{⟩}_s=\eta \left(\sum _{j=0}^{N_v-1}{\tilde{F}}_j^{\prime }|0{⟩}_r|j{⟩}_v+\tilde{E}|1{⟩}_r|0{⟩}_v\right),$$

$$\hat{C}:=H_{bq}\left(\stackrel{\curvearrowleft }{\prod _{z=1}^{L/2}}{B}_b({\theta }_z)V^{†}{B}_b({\varphi }_z)V{B}_b({\vartheta }_z)\right)H_{bq},$$

$$⟨0{|}_a⟨j{|}_{s}U|0{⟩}_a|k{⟩}_s=\beta {\hat{H}}_{jk}:={\hat{H}}_{jk}^{\prime }.$$

$$\begin{array}{c}U:=U_{\text{row}}^{†}{U}_{\text{col}},\\ U_{\text{row}}:=\sum _{i,j}|{\chi }_{i,j}{⟩}_{as}⟨i{|}_a⟨j{|}_s,\\ U_{\text{col}}:=\sum _{i,k}|{\psi }_{i,k}{⟩}_{as}⟨i{|}_a⟨k{|}_s.\end{array}$$

$${\hat{H}}_{jk}^{\prime }=⟨{\chi }_{0,j}|{\psi }_{0,k}⟩.$$

$$R(\varrho )|0⟩\to \varrho |0⟩+\sqrt{1-|\varrho {|}^2}|1⟩$$

$${\hat{H}}^{\prime }=\sum _{j=0}^{N_v-1}\left[{\Omega }_j|j⟩⟨j|+{{\rm Y}}_j\left(|j⟩⟨N_v|+|N_v⟩⟨j|\right)\right]+\hat{D}.$$

$${\Omega }_j=c^2{d}_j^2=k\beta v_j,$$

$${{\rm Y}}_j=\sqrt{1-|c{|}^2}\frac{d_j{b}_j}{\sqrt{N_v}}=\beta v_j{\alpha }_j,$$

$$d_j$$

Full text

I Introduction

Quantum computers show enormous promise for solving classes of problems for which a quantum algorithm can obtain an exponential speedup over the classical counterpart. Naturally, this includes the simulation of quantum systems Feynman (1982); Babbush et al. (2018), but it also includes problems that are not inherently quantum, such as integer factorization Shor (1997) and solving linear systems Harrow et al. (2009); Childs et al. (2017). However, it remains unclear how many classical computational tasks requiring the largest supercomputersKothe (2007) could be sped up using a future error-corrected quantum computer.

One such computationally extreme example is the kinetic plasma problemLapenta et al. (2015); Howard et al. (2015); Ku et al. (2018). High-temperature plasmas occur in a wide range of applications, including fusion, space, solar, and astrophysical contexts. Such plasmas are more complex than conventional fluids because near-range forces are Debye shieldedHückel and Debye (1923); Nicholson (1983). This leads to dynamics far from thermal equilibrium, requiring the time evolution of a smooth, six-dimensional distribution function $f({\bf x},{\bf v},t)$, describing the phase-space density of particles.

High-temperature plasmas are fundamentally described by the Vlasov-Maxwell system of equations which arise from Newton’s second law and Maxwell’s equations, neglecting short-range Coulomb interactions. The Vlasov equationVlasov (1938); Nicholson (1983) is an excellent testbed for applying quantum algorithms to classical computational mathematics because it is a hyperbolic partial differential equation describing a conservation law (conservation of particles) that has features similar to equations in computational fluid dynamics and other applicationsLeVeque (1992).

We study the Vlasov-Maxwell system and find that, when linearized about a Maxwellian equilibrium, the time evolution can be expressed in a unitary fashion. This is then mathematically equivalent to the evolution of a quantum system. On a universal quantum computer, the time evolution of arbitrary quantum systems is achieved with Hamiltonian simulation algorithms. These algorithms have undergone significant development in recent yearsBerry et al. (2015); Berry and Novo (2016); Novo and Berry (2017), culminating in the achievement of optimal query complexity for sparse Hamiltonians Low and Chuang (2017). That may be enough for the simulation of physical quantum systems, but for more general computations the efficient handling of non-sparse Hamiltonians can be important. On this front there has also been progress: a number of algorithms which do not rely on sparsity have been developedChilds (2010); Low and Chuang (2019); Wang and Wossnig (2018). However, efficient simulation is not possible for all non-sparse HamiltoniansChilds and Kothari (2010).

The efficiency of a quantum algorithm that time-evolves the linearized Vlasov-Maxwell system is unclear since the associated Hamiltonian is not sparse. Additionally, encoding initial data into a quantum state and extracting final data out can be expensive. To explore these issues, we fully detail an algorithm for a simple limiting form of the full system. This will correspond to the electrostatic Landau damping problem, which is an important kinetic problem in plasma physics. Although the simulation of Landau damping is not computationally demanding, we show that the quantum algorithm scales differently than its classical counterpart. That could translate into a large speedup in generalizations of the algorithm, which are discussed. This paper is intended as a first step in exploring the potential of quantum computation applied to plasma physics problems.

An outline of the paper follows. In Sec. II we analyze the linearized Vlasov-Maxwell system. We then focus on the electrostatic Landau damping problem. In Sec. III we detail the Hamiltonian simulation algorithm used to solve the Landau damping problem. In Sec. IV the results of numerical tests (using a classical computer) are given, and in Sec. V we discuss initialization and measurement. Section VI explores possible extensions of the current algorithm, and finally we summarize our findings in Sec. VII.

II The Vlasov-Maxwell System

We begin with the nonrelativistic Vlasov equation:

[TABLE]

where $f_{s}=f_{s}({\bf x},{\bf v},t)$ is the distribution function for species $s$, i.e. electrons or a type of ions. We solve for the behavior of small-amplitude waves by linearizing Eq. (1) and assuming a Maxwellian background distribution $f_{M}({\bf v})$, at rest. The zeroth-order equation is satisfied only if the background electric field ${\bf E}_{0}$ vanishes. The first-order equation for electrons is

[TABLE]

Assuming stationary ions, the time evolution is governed by Eq. (2) and Maxwell’s Equations. To adapt these to a quantum algorithm we apply a few transformations. We Fourier-transform in space, introducing the wavevector ${\bf k}$, and switch to these dimensionless variables:

[TABLE]

where $T_{e}$ is the electron temperature, $\omega_{pe}$ is the electron plasma frequency, $n_{e}$ is the electron number density, and $\lambda_{De}$ is the Debye length with ions neglected. Applying these changes, Eq. (2) becomes

[TABLE]

where all quantities are now the dimensionless versions, $\sim$ identifies variables that are Fourier components for wavevector ${\bf k}$, the $1$ subscripts have been dropped, and $c_{n}:=\frac{c}{\lambda_{De}\omega_{pe}}$ is the speed of light expressed in units of $\lambda_{De}\omega_{pe}$. Applying the same transformations to Maxwell’s Equations gives

[TABLE]

[TABLE]

Equations (3), (4), and (5) are sufficient to time-evolve the variables $\tilde{\bf E}$, $\tilde{\bf B}$, and $\tilde{f}({\bf v})$. To employ Hamiltonian simulation algorithms, we also need the system’s time evolution to be unitary. Consider the following change of variables:111An arbitrary $e^{i\phi({\bf v})}$ factor can be included in Eq. (6) without breaking the unitarity of the time evolution.

[TABLE]

and the real constant $\xi$ will be chosen later. This turns Eqs. (3) and (4) into

[TABLE]

[TABLE]

where $\mu({\bf v}):=\sqrt{\xi f_{M}({\bf v})}$. The time evolution is unitary if it can be written as $d|x\rangle/dt=-i\hat{H}_{\text{eff}}|x\rangle$ where $\hat{H}_{\text{eff}}=\hat{H}_{\text{eff}}^{\dagger}$ is the effective Hamiltonian of the system, and $|x\rangle$ is a quantum state with the variables encoded as amplitudes. Now, we consider the terms in Eqs. (5), (7), and (8). First, the evolution within the velocity space is generated by $\hat{H}_{v}$ with

[TABLE]

This is Hermitian since

[TABLE]

for normalizable functions $g({\bf v})$ and $h({\bf v})$, using integration by parts. The behavior of the last term in Eq. (9) is further explained by Eq. (53). Next, defining ${\bf k}:=(0,0,k)$, the electromagnetic (EM) part of the evolution can be expressed as

[TABLE]

Clearly, $\hat{H}_{l}$ is a Hermitian matrix. Lastly, there is the coupling between the velocity distribution and the electric field. Here the Hermitian restriction can be expressed in terms of coupling constants $C_{p}({\bf v})$:

[TABLE]

where $E_{p}$ is the $p^{th}$ component of ${\bf E}$. From Eq. (7) we find $C_{p}({\bf v})=\mu v_{p}$ while from Eq. (8) we get $C_{p}({\bf v})=\mu v_{p}d^{3}v/\xi$. Consequently we need to choose $\xi=d^{3}v$. In practice we will have a uniform, finite grid in velocity space in which case $\xi$ will be the volume of a single velocity grid cell. Now, we can formally consider the coupling constants to be components of an operator $\hat{H}_{c}$, which acts on a space that is the direct sum of the spaces in which $\tilde{f}^{\prime}({\bf v})$ and ${\bf E}$ are specified. Then Eq. (12) says that $\hat{H}_{c}$ is Hermitian, so $\hat{H}_{\text{eff}}=\hat{H}_{v}+\hat{H}_{l}+\hat{H}_{c}$ is Hermitian, and thus the evolution generated by $\hat{H}_{\text{eff}}$ is unitary.

II.1 The Landau Damping Problem

Now that we have unitary time evolution, in principle Hamiltonian simulation can be applied to time-evolve Eqs. (5), (7), and (8), but designing an efficient quantum algorithm is still a challenging task. We will demonstrate an algorithm for the electrostatic case with ${\bf B}_{0}={\bf 0}$. This still has the coupling term [Eq. (12)], which is tricky to handle efficiently, and it captures the physics of Landau dampingLandau (1946). In developing the quantum algorithm we target and achieve costs scaling only logarithmically with the grid size, in contrast to the usual linear scaling of classical algorithms. Although a large grid will not be necessary in our simulation of Landau damping, extensions of the algorithm to handle problems that involve very large grids are conceivable. These are discussed in Sec. VI.

In the electrostatic case there is only an electric field along ${\bf k}$, i.e. only $E_{z}$ is nonzero. Moreover, the electric field evolution can be determined using only a 1D velocity space: the velocity space, represented by a uniform 3D grid indexed by ${\bf j}$, can be reduced to 1D by introducing

[TABLE]

[TABLE]

in terms of which Eqs. (3) and (4), now with ${\bf B}_{0}={\bf 0}$, become

[TABLE]

[TABLE]

where $\alpha_{j}:=\sqrt{\Delta vG_{j}}$, and all quantities that were previously vectors are now their $z$ components, e.g. $j:={\bf j}_{z}$. To apply Hamiltonian simulation, we encode the data as the amplitudes of a quantum state, written in bra-ket notation as

[TABLE]

where $N_{v}$ is the total number of velocity grid cells, and $\eta:=\left(|\tilde{E}|^{2}+\sum_{j=0}^{N_{v}-1}|\tilde{F}^{\prime}_{j}|^{2}\right)^{-1/2}$ normalizes the state. This state evolves via $d|x\rangle/dt=-i\hat{H}|x\rangle$ with a Hamiltonian given by

[TABLE]

III Hamiltonian Simulation

Efficient Hamiltonian simulation of Eq. (18) is somewhat difficult due to the matrix being non-sparse in the sense often assumed for quantum algorithms. Specifically, many Hamiltonian simulation techniques (e.g. Berry et al. (2007); Berry and Novo (2016); Novo and Berry (2017)) have costs that scale polynomially with the maximum number of nonzero entries in any row. Our Hamiltonian has $N_{v}$ nonzero entries in one row due to how the electric field couples with all velocities.

To avoid a $\text{poly}(N_{v})$ cost scaling, we apply the Hamiltonian simulation technique developed by Low and Chuang (2017) (detailed further in Low and Chuang (2019, 2017); Low et al. (2016)). This technique allows the Hamiltonian to be encoded in a general way, which creates possibilities for efficiently handling non-sparse cases. It also has optimal query complexity, scaling linearly with time plus a term logarithmic in the error. Still, we must show how the individual queries can be performed efficiently. In what follows, we detail our Hamiltonian simulation implementation and analyze its efficiency.

Our algorithm operates on a quantum state with $n_{v}+7$ qubits where $n_{v}:=\log_{2}(N_{v})$. We divide this into five registers. Three of these, labeled by $b$, $q$, and $r$, are single-qubit registers, $a$ has four qubits, and $v$ has $n_{v}$ qubits. The main data are stored in the $r$ and $v$ registers as

[TABLE]

where $s$ is used to denote the combined $r$ and $v$ registers. It is on this $s$ register that the Hamiltonian [Eq. (18)] acts, and the $s$ state components with indices in $(N_{v},2N_{v})$ are unused. The other registers have ancilla qubits. They are initialized as $|0\rangle$, and at the end of the algorithm if any of them are measured to be nonzero then the algorithm has failed and must be rerun.

We describe the algorithm from the top down. Only the query implementation $U$ at the bottom is of our design. The higher operations are not specific to our Hamiltonian and are described by Low and Chuang (2019), but we include them for completeness. Hamiltonian simulation is performed by the circuit

[TABLE]

where $H_{bq}$ denotes a Hadamard gate, $\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\ 1&-1\end{pmatrix}$, applied to both the $b$ and $q$ qubits; $B_{b}(\Phi)$ is the phase shift gate, $\begin{pmatrix}1&0\\ 0&e^{i\Phi}\end{pmatrix}$, applied to the $b$ qubit; and $V$ is the quantum circuit depicted in Fig. 1. The ↶ over the product means that it is expanded in right-to-left fashion, e.g. the $z=1$ term is the first circuit operation. The angles used in the phase shift gates are determined222$\theta_{j+1}:=\pi+\varphi_{2j+1}$, $\vartheta_{j+1}:=-\varphi_{2j}$, $\phi_{j}:=-\theta_{j}-\vartheta_{j}$ by the phase vector $\vec{\varphi}$ which is computed as explained by Low et al. (2016). $\vec{\varphi}$ depends only on the simulation time and error threshold, and its length $L$ is proportional to the query complexity. The determination of $L$ is described later, leading to the query complexity estimate of Eq. (39).

We define a query to be a single call to the operator $U$ or $U^{\dagger}$. The Hamiltonian is encoded in the top-left block of $U$:

[TABLE]

The constant factor $\beta$ occurs because we are limited to encoding matrix rows and columns with normalizations of at most one. $\beta$ just amounts to a rescaling of the simulation time, which is corrected for by simulating for a time $t^{\prime}:=t/\beta$. To implement $U$ we use a strategy from Low and Chuang (2017) of breaking it into components $U_{\text{row}}$ and $U_{\text{col}}$, each of which is a unitary state preparation operator:

[TABLE]

The Hamiltonian that will be simulated is determined by the prepared states:

[TABLE]

There are still many possibilities for implementing the state preparation operators. Our choices are shown in Figs. 3 and 4. We introduce a rotation operator $R(\cdot)$, taken to be in SU(2), defined by its action:

[TABLE]

for some $|\varrho|\leq 1$. When $\varrho$ is an efficiently-computable function of the qubits of a different register we call this a variable rotation. We only use values of $\varrho$ that are purely real or imaginary. This corresponds to applying $e^{-i\hat{\sigma}_{y}\arccos{\varrho}}$ or $e^{-i\hat{\sigma}_{x}\arccos{\text{Im}(\varrho)}}e^{i\hat{\sigma}_{z}\pi/2}$, respectively. By computing the rotation angles in temporary registers and applying rotations controlled on the angle qubits, these can be implemented efficiently. For our variable rotations, the input register has $n_{v}$ qubits, and the angles can be computed at a cost poly($n_{v}$) [assuming for simplicity that they are computed to poly($n_{v}$) bits of precision], so the variable rotations also cost $\text{poly}(n_{v})=\text{polylog}(N_{v})$.

The Hamiltonian implemented by $U$ takes the form

[TABLE]

Here $\hat{D}:=\sum_{j=1}^{N_{v}-1}\Omega_{j}|j+N_{v}\rangle\langle j+N_{v}|$ is unimportant since it acts within the unused subspace of the state. We obtain our desired Hamiltonian [Eq. (18)] provided

[TABLE]

[TABLE]

where $c$, $b_{j}$, and $d_{j}$ specify the rotations as depicted in Figs. 3 and 4. Since we also have the restriction $|b_{j}|,|d_{j}|<1$, we choose

[TABLE]

where $v_{\text{max}}:=\max_{j}|v_{j}|=(N_{v}-1)\Delta v/2$ and $g_{\text{max}}:=\max_{j}|v_{j}G_{j}|$. Substituting Eq. (31) into Eqs. (29) and (30) gives an equation for $c$:

[TABLE]

where $\Gamma:=(k^{2}v_{\text{max}})/(\Delta vN_{v}g_{\text{max}})$. The last unknown, $\beta$, is given by $\beta=c^{2}/(kv_{\text{max}})$.

Since we must simulate $\hat{H}^{\prime}$ for a time $t^{\prime}:=t/\beta$, the simulation costs scale with $1/\beta$. Consequently, the dependence of $\beta$ on problem parameters, including $k$, $N_{v}$, and $\Delta v$, is important. One can show333using $1\leq(1+\sqrt{1/\Gamma})(\sqrt{1+4/\Gamma}-1)\Gamma/2\leq 5/4\;\forall{\Gamma\geq 0}$ that

[TABLE]

where

[TABLE]

Crucially, since $\Delta vN_{v}=2v_{\text{max}}+\Delta v$, $\Lambda$ does not increase with decreasing $\Delta v$. Therefore using an extremely high-resolution grid does not increase the simulation query complexity. Meanwhile, the space and gate complexity scale only logarithmically with the grid resolution through the number of qubits $n_{v}$.

Since $\|\hat{H}^{\prime}\|\leq 1$ is required by the Hamiltonian simulation technique, the ideal $1/\beta$ would be $\|\hat{H}\|$. An upper bound of $\|\hat{H}\|$, obtained by adding the norm of the diagonal part to the norm of the off-diagonal part, is

[TABLE]

Note that this expression turns into the $\Lambda$ expression if we replace $v_{j}^{2}G_{j}$ with $v_{\text{max}}g_{\text{max}}=v_{\text{max}}\max_{j}|v_{j}G_{j}|$. So $\Lambda>\Lambda^{\prime}$ but in most scenarios the difference is only a factor of order one. The variable rotations used to implement $U$ move some of the input state out of the good subspace, which leads to a smaller $\beta$. More general state preparation techniques (e.g. Soklakov and Schack (2006); Grover and Rudolph (2002)) could improve (increase) $\beta$, but for our purposes the simpler implementation choices are sufficient.

IV Numerical Validation

For testing, a code that implements the Hamiltonian simulation gates as linear algebra operations performed on a classical computer has been developed. This is compared to the result obtained by directly computing $e^{-i\hat{H}t}|x\rangle_{s}$ and to theoretical calculations. We test with the following parameter choices:

[TABLE]

where $t_{0}=0$ denotes the initial time, and the initial electric field is from Poisson’s equation in 1D. The evolution of the electric field in this case is illustrated in Fig. 5.

The quantum algorithm simulates the entire time period as one operation, accepting as parameters the total time and a maximally allowed error, with the error $\varepsilon$ for our circuit $\hat{C}$ defined to be

[TABLE]

where $\hat{P}_{g}$ is the projector onto the good subspace, specifically the space with all ancillas being zero, and $|x_{g}\rangle$ is any input within the good subspace, i.e. $\hat{P}_{g}|x_{g}\rangle$ = $|x_{g}\rangle$. For a given error tolerance $\epsilon$ the circuit can be designed to guarantee $\varepsilon\leq\epsilon$, based on the error bounds of the Hamiltonian simulation technique. We validate our implementation by running with a range of error tolerances. The results are displayed in Fig. 6.

The query complexity of our algorithm is $2L$ where $L=2(n-1)$ and the positive integer $n$ is chosen such that the error boundLow and Chuang (2017),

[TABLE]

is below the chosen tolerance $\epsilon$. In our implementation we invert $b(n)\leq\epsilon$ numerically, but in the limit of $t^{\prime}\gg 1$ a simple approximation is also effective. Note that $b(n)\leq e^{-\Delta n}b(n_{0})$ with $n_{0}:=et^{\prime}/2+C$, $\Delta n:=n-n_{0}\geq 0$, and $C\geq 0$. This suggests taking $n=n_{0}+\ln(1/\epsilon)$. Plugging in $\epsilon=1$ gives the bound $C\geq n-e(n!/32)^{1/n}$. The simple choice of $C=1.5$ satisfies this for all $n$. Then $b(n)<\epsilon$ for any $\epsilon\leq 1$ provided that

[TABLE]

which leads to a query complexity $Q$ with

[TABLE]

A somewhat stronger asymptotic bound, obtained by Low and Chuang (2019), is

[TABLE]

but for $t^{\prime}\gg 1$ and reasonable values of $\epsilon$ Eq. (39) is fairly accurate, and the $t^{\prime}$ term dominates. Our test case has $t^{\prime}\approx 105.7$. The main features, shown in Fig. 6, of total query counts around 600 and error decreasing by approximately a factor of $e$ for every four additional queries, are consistent with Eq. (39). All this does not depend on our specific Hamiltonian. However, as shown in Sec. III, the queries for our Hamiltonian can be implemented efficiently on a quantum computer, i.e. with costs scaling only as polylog($1/\Delta v$), while classically they cost $\mathcal{O}(1/\Delta v)$.

Having verified that the Hamiltonian simulation is accurate, now we interpret the results. Physically, the system has Langmuir oscillations that decay due to Landau damping. As Fig. 5 shows, after a brief initial stage (i.e. the first period or so) the electric field is well described by a damped sinusoidal. Specifically, we can fit the function $Ae^{-\gamma t}\cos(\omega t-\rho)$ to obtain parameters of interest, namely the frequency $\omega$ and damping rate $\gamma$:

[TABLE]

Simple theoretical estimatesNicholson (1983), translated into our dimensionless variables and with $k=0.4$, give

[TABLE]

These estimates rely on $k<<1$, so they can only give rough agreement. Precise values of $\gamma$ and $\omega$ can be found by numerical integration of the dispersion relation; in our units: $1+\frac{1}{\sqrt{2\pi}k^{2}}\int_{-\infty}^{\infty}\frac{v\exp(-v^{2}/2)dv}{v-(\omega-i\gamma)/k}=0$. This gives $\omega=1.28506$ and $\gamma=0.06613$, in agreement with the fit results. These comparisons show that our algorithm accurately reproduces Landau damping.

Note that this test case was chosen to be simple and easy to compute. The number of operations is not large, both for classical and quantum algorithms. Even the classical code that simulates the quantum computation completes in around a second. This is possible since we used a low-resolution velocity space, i.e. $N_{v}=32$. The classical costs would also be much higher for the multidimensional generalizations discussed in Sec. VI.

V Initialization and Measurement

On a quantum computer, we would perform three steps: state preparation, Hamiltonian simulation, and measurement. State preparation costs depend on the prepared state, and in the worst case these costs scale with the state space size, but many classical probability distributions can be prepared efficientlySoklakov and Schack (2006); Grover and Rudolph (2002). For simulating Landau damping we use an initially Maxwellian $\tilde{F}$, and the preparation can be done efficiently with Fig. 7. That circuit succeeds with probability $\langle|h_{j}|^{2}\rangle$ where

[TABLE]

Inputting $\tilde{F}^{\prime}_{j}\propto\tilde{F}_{j}/\sqrt{G_{j}}\propto e^{v_{j}2/4}$ one finds $1/\langle|h_{j}|^{2}\rangle=\mathcal{O}(v_{\max{}})$ so Fig. 7, which has a gate complexity $\text{poly}(n_{v})=\text{polylog}(N_{v})$, needs to be repeated $\mathcal{O}(v_{\max{}})$ times444or $\mathcal{O}(\sqrt{v_{\max{}}})$ rounds of amplitude amplificationBrassard et al. (2002). before running the rest of the algorithm. Note that this cost just adds to the overall algorithm cost and does not scale with the simulation time. Consequently, it does not change the algorithm’s asymptotic complexity.

Obtaining output from the algorithm is more difficult. In general quantum state tomography can be applied to learn the full state, but we do not want a scaling with the state space size. If we just want the electric field, that is a single amplitude, so the technique of amplitude estimationBrassard et al. (2002) is appropriate. This can produce an estimate $\tilde{p}$ of the probability $p$ of measuring the state to be within a specified subspace, with error bounded by

[TABLE]

where $M$ is the number of iterations applied. Each iteration involves running the original algorithm forwards and backwards once, thus the full costs are multiplied by $2M$. From Eq. (44), $M=\mathcal{O}(1/\delta)$, a quadratic speedup over direct sampling which would instead require $\mathcal{O}(1/\delta^{2})$ repetitions to estimate the outcome probabilities with a standard deviation of $\sigma=\Theta(\delta)$.

If we want to efficiently estimate some probability to a specified relative accuracy, we also need to ensure that it is not too small, since the costs scale with the absolute accuracy $\delta$. For our initially Maxwellian $\tilde{F}$, and with sums over velocity space approximated by integrals, one finds

[TABLE]

Since this is not small, if we want to measure $\tilde{E}$ to a fixed accuracy $\delta^{\prime}$ relative to the sinusoidal envelope, the measurement adds an $\mathcal{O}(1/\delta^{\prime})$ cost factor provided the simulation time $t$ is chosen such that $e^{\gamma t}=\mathcal{O}(1)$, i.e. the electric field is only moderately damped.

Both direct sampling and amplitude estimation will only give information about magnitudes such as $|\tilde{E}|$. In general one may want information about the phase of a complex amplitude. This can be obtained by extending the original algorithm as shown in Fig. 8. We introduce another qubit, labeled by $c$, and a phase $\zeta$. If the original algorithm produced an amplitude $\nu$ for the $|0\rangle_{m}$ basis state, the extended algorithm modifies this via

[TABLE]

Then the usual measurement process, e.g. amplitude estimation, can be applied to estimate $\left|\nu+e^{i\zeta}\right|$. By repeating this for $\zeta=\{0,\pm 2\pi/3\}$, we can determine the location of $\nu$ within the complex plane, rather than just its magnitude. Explicitly, if we write the new measured magnitudes as

[TABLE]

we can solve to obtain

[TABLE]

If each of the squared magnitudes is measured to an accuracy $\delta$, the same $\delta$ as in Eq. (44), then we obtain an estimate $\tilde{\nu}$ of $\nu$ which obeys the bound $|\tilde{\nu}-\nu|\leq 8\delta/3$. Thus we can learn the value of $\nu$ within the complex plane to within a distance $\delta$ at a cost $\mathcal{O}(1/\delta)$. This can be applied to any chosen basis state, not just $|0\rangle_{m}$, by specifying the target state when performing amplitude estimation. For instance, we could use this technique to obtain the sign of Im($\tilde{E}$).

Now that all components of the algorithm have been discussed, we can determine its computational complexity, assuming the output is the final electric field. In principle the scaling with respect to any algorithm parameter can be considered, but some parameters are more relevant than others. When $k\gtrsim 1$ the electric field is rapidly damped away, so we assume $k=\mathcal{O}(1)$. Next, since the velocity components beyond $v_{\max{}}$ decay as $e^{-v_{\max{}}^{2}/2}$, a large $v_{\max{}}$ is generally unnecessary. In contrast, a large $N_{v}$ is needed to resolve the fine velocity-space oscillations that develop due to the $-ikv_{j}\tilde{F}^{\prime}_{j}$ term in Eq. (15). Therefore, we treat $v_{\max{}}$ as a constant and give the scaling with respect to $N_{v}$, the simulation time $t$, and the error tolerances $\epsilon$ and $\delta$. Taking $\epsilon=\delta$, the asymptotic gate complexity is

[TABLE]

If we assume that $t>\log(1/\delta)$, which is true for large enough $t$ and reasonable values of $\delta$, this simplifies to $\mathcal{O}(\text{polylog}(N_{v})t/\delta)$.

The final electric field is just a simple example of an output for the algorithm. Given that $\tilde{E}$ oscillates in time, a more appropriate output may be

[TABLE]

where the times $t_{i}$ are distributed over some time range of interest. In fact, this can be obtained by a simple extension of the original algorithm. Suppose each $t_{i}$ can be efficiently computed from its index $i$, and the times can be represented with integers $\tilde{t}_{i}:=2^{m}t_{i}/t_{\max}-1$ where $t_{\text{max}}:=\max_{i}t_{i}$. We can add an index register $i$ in a uniform superposition of the indices and a time register $t$, with $m$ qubits, in which the $t_{i}$ are computed. Then, controlled on each time qubit, the original algorithm is run for the time interval specified by that bit. The query bound of Eq. (39) changes to $Q\leq 2et_{\text{max}}^{\prime}+m(4\ln(1/(m\epsilon))+6)$. After the Hamiltonian simulation, the state normalization within the subspace spanned by $\{|i\rangle_{i}|t_{i}\rangle_{t}|N_{v}\rangle_{s}\}\forall{i}$ and all ancilla qubits $|0\rangle$ will be $\eta\tilde{E}_{\text{rms}}$ (to within error $\epsilon$). Amplitude estimation can be applied to determine this normalization with a better scaling than direct sampling. Thus we can get $\tilde{E}_{\text{rms}}$ with the same dominant cost scaling as the original algorithm, i.e. $\mathcal{O}(\text{polylog}(N_{v})t_{\text{max}}/\delta)$.

VI Generalizations

In Sec. II we initially worked in 3D velocity space, with a uniform ${\bf B}_{0}$ and general field perturbations ${\bf E}_{1}$ and ${\bf B}_{1}$. In principle that problem can still be solved efficiently on a quantum computer using Hamiltonian simulation. A difficulty is that the dropped terms are large: $c_{n}\gg 1$ for nonrelativistic plasmas and $\|\frac{\partial}{\partial{\bf v}}\|=\mathcal{O}(1/\Delta v)$. Consequently, the spectral norm of the Hamiltonian would be dramatically increased by including these terms, and Hamiltonian simulation costs scale with that norm. However, there is a potential solution: we can break the Hamiltonian into the new large piece $A$ and the old small piece $B$ and apply the Hamiltonian simulation technique of Low and Wiebe (2018). In that paper, the interaction picture is applied along with time-dependent Hamiltonian simulation to allow separate handling of the two components. We would take $A$ to be the Hamiltonian in $d|x\rangle/dt=-iA|x\rangle$ where the variables are encoded in $|x\rangle$ and are evolving via

[TABLE]

[TABLE]

while $B$ would be the Hamiltonian for the other terms of Eqs. (7) and (8).

As is, our algorithm handles evolution by $B$ in the 1D case. The extension of the velocity space to 3D is simple: the 1D velocity index is replaced by a flattened 3D velocity index. The matrix would have $\mathcal{O}(1/\Delta^{3}v)$ nonzero entries in one row and column, but the coupling constants of Eq. (12) scale as $\sqrt{\Delta^{3}v}$ so we still have $\|B\|=\mathcal{O}(1)$. Then we can obtain similar efficiency, i.e. costs linear in time and logarithmic in error and grid resolution, provided that we can efficiently implement $e^{iA\tau}$ with $\tau=\mathcal{O}(1)$, corresponding to the time evolution by the added terms alone.

The EM evolution of a single Fourier mode is given by Eq. (11). This is trivial to solve; the pairs $\{E_{y},B_{x}\}$ and $\{E_{x},B_{y}\}$ each undergo rotations about the x-axis in their Bloch spaces. Therefore rotation gates can be applied to perform this evolution, circumventing the usual scaling with spectral norm. We still need a way to efficiently perform the evolution of Eq. (9) (which commutes with the EM evolution). We can reexpress this using

[TABLE]

where ${\bf\hat{L}_{v}}:=\frac{{\bf v}}{i}\times\frac{\partial}{\partial{\bf v}}$ is the velocity space analog of the angular momentum operator. Thus the evolution caused by this term is simply the rotation of velocity space around ${\bf B}_{0}$, i.e. it generates cyclotron motion. Once again we can hope to implement this efficiently in a direct way, avoiding spectral norm scaling. For instance, if we work on a cylindrical velocity grid, symmetric about the ${\bf B}_{0}$ axis and with velocity cells spaced uniformly in $\phi$, then certain discrete time steps just correspond to addition of a binary constant to the qubits associated with the $\phi$ index. We leave the details for potential future work.

In general the distribution function is 6D. Our algorithm applies to a single spatial Fourier mode ${\bf k}$, but it can easily be extended to many modes. For applying the algorithm to a superposition of modes, the only adjustment is that all quantities which depend on ${\bf k}$ need to be computed in superposition using operations controlled on the new ${\bf k}$ register(s). In principle this can still be done efficiently, but we will not work this out here. Additionally, if the initial variables are specified as functions of position, we could use the quantum Fourier transformNielsen and Chuang (2010) to efficiently convert this to a Fourier space representation before running the algorithm. Another consideration is the magnitude of ${\bf k}$. If we do not assume $|{\bf k}|\lesssim 1$ then the simulation costs go up since the Hamiltonian spectral norm scales linearly with $|{\bf k}|$. Physically, $|{\bf k}|\lesssim 1$ is justified in that this is enough to resolve $\lambda_{De}$. Moreover, when $|{\bf k}|\gtrsim 1$ the electrostatic waves are rapidly damped away, as suggested by Eq. (42). Then even with many modes, having $|{\bf k}|\lesssim 1\;\forall{\bf k}$ is reasonable.

Another significant restriction made in Sec. II was the assumption of a Maxwellian background distribution. This makes the system unconditionally stable. The alternative, where some growing wave modes exist, violates unitarity since then a state that holds the perturbed quantities must grow exponentially if it contains a growing mode component. Such cases cannot be handled directly with Hamiltonian simulation. A completely different algorithm, possibly utilizing a quantum linear systems algorithm, would be required, and it is unclear whether it would scale well. This is a topic of future work. Still, our algorithm could be applied to other stable background distributions. The background distribution enters Eqs. (5), (7), and (8) only through its velocity derivative. In our dimensionless variables, $\frac{\partial f_{M}}{\partial{\bf v}}=-{\bf v}f_{M}$. For a more general background distribution $f_{0}$ we can make the replacement $f_{M}\rightarrow g_{0}$ where

[TABLE]

By redoing the steps up through Eq. (12) one finds that the new Hamiltonian is Hermitian provided $g_{0}({\bf v})\geq 0\;\forall{\bf v}$. For background distributions that satisfy this restriction, Hamiltonian simulation can still be applied. Furthermore, in the electrostatic case and with ${\bf B}_{0}={\bf 0}$ the background distribution enters only through $\frac{\partial f_{0}}{\partial v_{z}}$. Here the appropriate replacement is

[TABLE]

and the unitarity restriction is $G_{j_{z}}\geq 0\;\forall{j_{z}}$. Our algorithm can be easily applied to such distributions; only the specific values of $b_{j}$ and $\beta$ would be altered.

Thus far we have solved for the motion of electrons, with ions assumed stationary. If we remove this restriction then each species obeys its own version of Eq. (1), and they are coupled through the fields. Keeping the same dimensionless variables, Eq. (3) for each species is

[TABLE]

where $r_{s}:=\frac{m_{e}q_{s}}{m_{s}e}$. If we absorb the $r_{s}$ factor into the background distribution and magnetic field, we recover the original mathematical form, and the prior manipulations are still applicable, including the discussed extensions. The rescaling factors $\mu(\bf v)$ become species dependent through the $r_{s}$ factor. One may be concerned about the sign change for positive species, given the importance of signs for unitarity, but since there is a corresponding sign change in Eq. (4), there are no new issues here. The generators of the evolution of each $\tilde{f}_{s}$, including the couplings with $\tilde{E}$, are anti-Hermitian, therefore the combined evolution of all species and fields is unitary. Additionally, the terms in the evolution, and thus the matrix entries, for other species are relatively small since $r_{s}\ll 1$ for ions. Consequently the Hamiltonian spectral norm, along with the space and time complexity of the algorithm, would be nearly unchanged by the inclusion of multiple species.

VII Summary

In high-temperature plasma physics, evolution of the particle distribution function is governed by the Vlasov-Maxwell system of equations, but simulating this classically with the full 6D phase space is extremely expensive. For instance, if the phase-space grid requires 1000 cells in each dimension, just storing $f({\bf x},{\bf v},t)$ is beyond exascale capabilityService (2018). However, on a quantum computer, the distribution function can be encoded as the amplitudes of a quantum state, requiring exponentially fewer qubits than classical bits. Moreover, we have shown that the linearized Vlasov-Maxwell system of equations (with $f_{0}=f_{M}$) produces a unitary evolution of the field variables and a rescaled version of the distribution function. Thus Hamiltonian simulation algorithms can be applied to perform this classical physics simulation.

For one limiting form of the full system we demonstrated in detail how the Hamiltonian simulation can be performed efficiently, such that the quantum algorithm gets an exponential speedup in space and time complexity over the classical simulation with respect to the velocity grid size. To obtain this efficiency we used a recently developed Hamiltonian simulation technique, and the efficient handling of the more general problem appears to require even more sophisticated Hamiltonian simulation techniques. An additional challenge is the extraction of an output. For measuring a value encoded in a single quantum amplitude to within absolute error $\delta$, $\mathcal{O}(1/\delta)$ iterations of the original algorithm are required. This can significantly diminish the quantum speedup, depending on the desired accuracy.

Acknowledgements.

Research supported in part by the United States Department of Energy under grant DE-SC0020393.