Quantum Computing For Fusion Energy Science Applications

Quantum computing for fusion energy
science applications
Cite as: Phys. Plasmas 30, 010501 (2023); https://doi.org/10.1063/5.0123765
Submitted: 01 September 2022 • Accepted: 15 December 2022 • Published Online: 27 January 2023
I. Joseph, Y. Shi, M. D. Porter, et al.
COLLECTIONS
Paper published as part of the special topic on Papers from the 2021-2022 Sherwood Fusion Theory Conferences
This paper was selected as Featured
ARTICLES YOU MAY BE INTERESTED IN
Progress toward fusion energy breakeven and gain as measured against the Lawson criterion
Physics of Plasmas 29, 062103 (2022); https://doi.org/10.1063/5.0083990
The importance of laser wavelength for driving inertial confinement fusion targets. I. Basic
physics
Editorial: Preface to the 30th volume of Physics of Plasmas

Phys. Plasmas 30, 010501 (2023); https://doi.org/10.1063/5.0123765 30, 010501
© 2023 Author(s).
Physics of Plasmas REVIEW scitation.org/journal/php
Quantum computing for fusion energy science

applications
Cite as: Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765
Submitted: 1 September 2022 . Accepted: 15 December 2022 .
Published Online: 27 January 2023
I. Joseph,a) Y. Shi, M. D. Porter, A. R. Castelli, V. I. Geyko, F. R. Graziani, S. B. Libby,

and J. L. DuBois
AFFILIATIONS
Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, California 94551, USA
Note: This paper is part of the Special Topic: Papers from the 2022 Sherwood Fusion Theory Conference.
a)
Author to whom correspondence should be addressed: joseph5@llnl.gov
ABSTRACT
This is a review of recent research exploring and extending present-day quantum computing capabilities for fusion energy science applications.
We begin with a brief tutorial on both ideal and open quantum dynamics, universal quantum computation, and quantum algorithms. Then, we
explore the topic of using quantum computers to simulate both linear and nonlinear dynamics in greater detail. Because quantum computers
can only efficiently perform linear operations on the quantum state, it is challenging to perform nonlinear operations that are generically
required to describe the nonlinear differential equations of interest. In this work, we extend previous results on embedding nonlinear systems
within linear systems by explicitly deriving the connection between the Koopman evolution operator, the Perron–Frobenius evolution operator,
and the Koopman–von Neumann evolution (KvN) operator. We also explicitly derive the connection between the Koopman and Carleman
approaches to embedding. Extension of the KvN framework to the complex-analytic setting relevant to Carleman embedding, and the proof
that different choices of complex analytic reproducing kernel Hilbert spaces depend on the choice of Hilbert space metric are covered in the
appendixes. Finally, we conclude with a review of recent quantum hardware implementations of algorithms on present-day quantum hardware
platforms that may one day be accelerated through Hamiltonian simulation. We discuss the simulation of toy models of wave–particle interac-
tions through the simulation of quantum maps and of wave–wave interactions important in nonlinear plasma dynamics.
C 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://
V
creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0123765
I. INTRODUCTION squeezed states, as sensitive probes can reduce the detection threshold
A. Motivation to the Heisenberg limit 1=S. Such techniques have been used to
improve the sensitivity of the Laser Interferometer Gravitational-wave
The three pillars of quantum information science (QIS)—quantum
sensing, quantum communications, and quantum computing (QC)— Observatory (LIGO) gravitational wave detector by a factor of 2 for
promise to have transformative impact on science, engineering, and nearly a decade.2 New detectors based on nitrogen vacancy centers
technology as we know it.1 This article presents a pedagogical introduc- (NV centers) in diamond have provided unprecedented sensitivity for
tion to quantum computing and reviews recent research to develop and measurements of magnetic and electric fields as well as temperature
apply quantum algorithms and utilize quantum computing hardware and pressure. Advances in atom interferometry3 have led to a revolu-
platforms for fusion energy science (FES) applications. The interesting tion in gravitational and inertial (gyroscopic) sensing, now transition-
results obtained so far make it hopeful that QIS may one day lead to ing to real world applications, such as navigation without the Global
game-changing capabilities for FES. Positioning System (GPS), rapid passive sensing of mass distributions,
Before diving into quantum computing (QC), let us briefly men- and underground structure discovery, as well as basic science applica-
tion the other two pillars of QIS. First, quantum sensing techniques are tions, such as gravitational wave detection in the frequency regime
already being used today to improve measurement sensitivity by between LIGO and the Laser Interferometer Space Antenna (LISA).
orders of magnitude. Instead of being limited
pffiffiffi by the central limit theo- Second, quantum communications offer the possibility of secure
rem to yield a noise-to-signal ratio of 1= S, where S is the number of and intrinsically parallel information transfer.4 The goal is to trans-
samples, using intrinsically quantum entangled states, such as form and transport quantum information over long distances and to
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-1

C Author(s) 2023
V
transduce information between different quantum hardware plat- distribution functions (PDFs), which are simply advected along the tra-
forms. This field generated some of the earliest rigorous proofs that jectories. While this is certainly true in the traditionally classical realm
combining entangled states with quantum communication protocols of plasma physics, it is also quite true in the traditional quantum realm
could surpass their classical counterparts in their ability to carry and of chemistry and materials science. The computational workhorses of
manipulate information. The application to network security is con- density functional theory (DFT) and molecular dynamics crucially
sidered to be so important that many researchers are already working rely on the ability to efficiently evaluate nonlinear functions and to
on the quantum internet. evolve nonlinear differential equations. The same is true for evolving
Finally, quantum computing, which is the focus of this review classical kinetic equations that describe the evolution of nuclear and
article, promises to efficiently perform quantum algorithms that have chemical species with reaction rate theory.
the power to achieve polynomial to exponential improvements in An essential limitation of the standard quantum computing para-
complexity while manipulating exponentially large quantum memory digm is that, aside from measurements, only linear unitary operations
resources.5–7 The recent growth in the capabilities of today’s quantum can be performed efficiently in the quantum Hilbert space. If one
computing devices has spurred great interest in quantum simulation, believes that quantum mechanics (QM) is the correct physical theory
and they have been used to perform key demonstrations of quantum of nature, then linear dynamics is all that is needed to describe the evo-
calculations.8–10 The potential impact of this so-called quantum lution of the universe. (A notable exception to this rule is that some
advantage is so exciting that many government and private industry physicists believe that the measurement process might require nonlin-
laboratories around the world are working to perfect the technology. ear dynamics to be described completely.) However, according to the
In fact, in recent years, there have been claims that some demonstra- exact factorization theorem14,15 the process of generating reduced
tions have already passed the point of quantum supremacy,10–12 the order models typically generates nonlinear interactions between sub-
point at which quantum computers can surpass even the world’s best systems. This is the reason that essentially all useful physical models
supercomputers at performing certain computational tasks. that achieve significant complexity reduction, e.g., fluid dynamics,
While there has certainly been enormous progress over the last density function theory (DFT), and kinetic theory, are nonlinear.
two decades, we shall see that the practical use of today’s hardware
platforms is still quite limited due to the lack of fault-tolerant error C. Quantum computing capabilities
correction. Because present-day and near-term quantum devices offer
Many scientific computations of interest to FES can be solved
large quantum memory registers but lack error correction, the present
more efficiently using quantum algorithms (see the excellent overviews
era has been dubbed the noisy intermediate-scale quantum (NISQ)
given in Refs. 16 and 17). Quantum algorithms promise to accelerate
era.13 The search for the ultimate physical platform for realizing a
many of the algorithms that are central to scientific computing, includ-
quantum computer spurs research onward to invent new quantum
ing the Fourier transform, sparse linear solvers, and sparse
materials, new hardware platforms, and new algorithms that can sur-
Hamiltonian simulation. Quantum algorithms also exist for linear
pass today’s limitations.
ODEs18 and linear PDEs19,20 that can be cast in the form of initial
value problems (IVPs), boundary value problems (BVPs), and eigen-
B. Fusion energy science needs value problems (EVPs).
Fusion energy science (FES) is defined by the mission of For example, for plasma and materials science applications, some
Achieving a safe sustainable fusion energy source for the foreseeable IVPs are naturally posed as being Hamiltonian and, hence, can natu-
future. Hence, it encompasses all fields of science required to achieve rally be presented as a unitary evolution resulting from a Hermitian
this goal, notably including plasma physics, nuclear physics, materials Hamiltonian. If the Hamiltonians of interest have a sparse structure,
science, and chemistry. Thus, the FES mission requires developing then they can be solved using efficient quantum algorithms for
accurate predictive models within all of these disparate fields, to form Hamiltonian simulation,21,22 which we will discuss in detail in Sec. IV.
combined whole-device models that ultimately span a wide range of Notably, Ref. 23 developed a quantum algorithm to solve the linear-
physical regimes from classical to quantum. ized Vlasov–Poisson system and Ref. 24 developed a quantum algo-
To stress the difficulty of this endeavor, FES requires calculations rithm to solve the problem of cold plasma wave propagation and
of strongly correlated quantum materials such as superconductors at linear mode conversion. Yet, the majority of IVPs for reduced models
temperatures below 1 meV to the collective motions of plasma at include dissipation. Since the dynamics is not Hamiltonian, they need
1–10 keV as well as nuclear reactions at 10s of MeV. Whereas mag- to be treated with specialized methods, for example, the general pur-
netic confinement fusion uses low-density plasmas on the order of pose quantum linear differential equation solver,18 discussed in Sec.
1020 m3, inertial confinement fusion seeks to obtain a compression of III E 7.
particle density on the order of 1000 times that of solid matter, for par- For many fields of science and engineering, including biology,
ticle densities beyond 1030 m3 . Thus, this represents a range of 9 chemistry, and physics, a large share of today’s computational resour-
orders of magnitude in energy and 10 orders of magnitude in density ces are used for the simulation of classical nonlinear dynamics. Hence,
(without considering the density of nuclear matter). it is important to understand whether quantum computers can pro-
Note, however, that, today, the vast majority of scientific com- vide similar gains in efficiency for the simulation of nonlinear prob-
putations in the fusion energy sciences are classical in nature. The lems. Nonlinear operations can be performed by forming the tensor
meaning of this statement is that such calculations seek to predict the product of identical states. However, due to the no-cloning theo-
time-dependent evolution of nonlinear sets of both ordinary differen- rem,25–27 it is not possible to make copies of a quantum state; rather,
tial equations (ODEs) and partial differential equations (PDEs). one must prepare multiple identical states from scratch. Thus, if a non-
Ensembles of systems are typically evolved using classical probability linear operation is to be iterated, e.g., in time, then this will require
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-2

C Author(s) 2023
V
preparing a number of identical states that is exponentially large in the dynamics, including generalized approaches to Koopman–von
number of iteration steps.28 Neumann (KvN) and Carleman embedding, as well as explain their
Efficient quantum algorithms for simulating nonlinear physical differences. The penultimate section, Sec. VI will discuss recent find-
processes have only recently begun to be invented and improved17,29–34 ings in applying several leading quantum computing hardware plat-
and some of these efforts were initiated for FES applications.16 The new forms to quantum simulation problems of interest to FES. Appendix B
algorithms all use the strategy of embedding the nonlinear system extends the KvN framework to the complex-analytic setting relevant
within an infinite-dimensional linear system that is then truncated to to Carleman embedding. Appendix C proves that different choices of
finite dimension. For systems of nonlinear ODEs, the general method of the relevant Hilbert space bases, i.e., reproducing kernel Hilbert spaces
using a quantum computer to solve for the evolution of the wavefunc- (RKHS), depend on the choice of Hilbert space metric. The concluding
tion corresponding to a classical PDF was first clearly described in Refs. section summarizes our findings and offers perspectives on the future
17 and 29, using what is known as the Koopman–von Neumann outlook for quantum computing for FES applications.
embedding. Using an alternative approach, Ref. 31 developed concrete
quantum algorithms for dissipative differential equations with quadratic II. QUANTUM INFORMATION
nonlinearity based on an approach known as Carleman embedding,
Benioff,48 Manin,49 and Feynman,50,51 were some of the first sci-
which is applied to simulate the Burgers equation. Connections between
entists to think seriously about the potential of using machines that
these embedding techniques are discussed in Ref. 32, and comparisons
satisfy the laws of quantum mechanics for performing scientific calcu-
of these approaches are elucidated using numerical examples in Ref. 33.
lations. Feynman also left us with one of the great quotes in the field:50
Building upon ODE solvers, quantum algorithms for solving PDEs,
Nature is not classical, dammit, and if you want to make a simulation
including the Navier–Stokes equations, have also been explored.35,36
of nature, you’d better make it quantum mechanical, and by golly it’s a
Apart from deterministic methods, algorithms based on random
wonderful problem, because it does not look so easy. Understanding
processes have also been developed. For example, quantum versions of
this statement is the key to understanding why scientists believe that
Monte Carlo algorithms can generally be performed with a quadratic
quantum computers can, in principle, be so powerful.
speedup over their classical counterparts.37 Similar concepts can then
be used to speed up calculations of stochastic differential equations
(SDEs).38 General multi-level Monte Carlo methods for stochastic dif- A. The postulates of mechanics
ferential equations can also been accelerated.39 In fact, quantum algo- To see why quantum computers are fundamentally different from
rithms based on these concepts have been proposed for computing classical computers, let us begin by recalling the postulates of quantum
turbulent mixing and the reaction rate within turbulent flows.40,41 mechanics6,7,52 and compare them with postulates of classical mechanics
The algorithms mentioned above requires fault-tolerant quantum (CM). A primer on Hilbert spaces is given in Appendix A. In the follow-
computers to implement various oracles that are required by quantum ing, we assume that the reader is familiar with Dirac’s bra-ket notation.
subroutines. Alternatively, on NISQ devices, which have been shown
to be capable of performing classical-quantum hybrid optimization 1. Postulates of quantum mechanics (QM)
problems,42,43 advancing nonlinear equations may be treated as an
optimization problem using variational algorithms.44 In this case, eval- In the following formulas, the Hilbert space will be treated
uating the nonlinear functional only requires a fixed number of copies as if it were discrete (which is typically true for digital models of quan-
per time step. Hence, there have been recent proposals for simulating tum computation), but if it is continuous, then the sums should be
the Navier–Stokes equations using the variational technique.45–47 replaced by integrals.
Even though simulating nonlinear problems is perhaps one of
1. Physical states.
the most useful applications of scientific computation, relevant quan-
(a) Pure physical states, denoted by the wavefunction w, are
tum algorithms are still in their infancy and there is much work to be
rays in Hilbert space.
done in terms of improving their speed and accuracy and verifying • Pure states are vectors in Hilbert space, generically, a
their performance on quantum computing hardware platforms. Once
linear superposition of basis elements that span the
fully vetted and understood, they can be used as subroutines to provide
Hilbert space.
a foundation for quantum programs that perform more useful and • Pure states are elements of projective Hilbert space; i.e.,
more complex tasks.
the overall complex constant is unimportant, and the
normalization can be set to unity, jjwjj ¼ 1.
D. Overview of contents • The Hilbert space is typically a space of functions that is
Section II will provide a brief introduction to the principles of associated with a configuration manifold, m, but it can
quantum information and explain the reason that quantum memory also be a discrete set.
registers can be considered to be much larger than their classical coun- (b) Mixed physical states are probability distributions over
terparts. Then, Sec. III will discuss the digital circuit model of universal pure states and, hence, are represented by a positive self-
quantum computation and review the key quantum subroutines that adjoint operator on Hilbert space, q† ¼ q, called the prob-
are useful for accelerating classical computations and for building ability density matrix.
quantum programs. Next, Sec. IV will discuss quantum approaches to • Due to Hermiticity, mixed states can be diagonalized
simulating differential equations, the type of application that is used X
most by FES researchers today. Then, Sec. V will unify several recent q¼ pj jwj i hwj j: (1)
proposals for using quantum algorithms to simulate nonlinear j
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-3

C Author(s) 2023
V
Because the eigenvalues, pj 0, must be non-negative • Ideal closed systems evolve via a linear unitary,
this represents a probability distribution over a set of U1 ¼ U† , transformation of Hilbert space,
pure states, jwj i, with probability pj .
• Probability is normalized to unity, wðtÞ ¼ Uwðt0 Þ qðtÞ ¼ Uqðt0 ÞU† : (7)
X
Tr ðqÞ ¼ pj ¼ 1: • Continuous time evolution is generated by a Hermitian
j Hamiltonian operator, H
• Proper mixed states must satisfy Tr ðqa Þ < 1 for a > 1
ih@t w ¼ Hw ih@t q ¼ ½H; q; (8)
and Tr ðqa Þ > 1 for 0 < a < 1.
• Mixed stated are required to represent the outcome of
where the commutator is defined via
measurements.
(c) The Hilbert space of a composite system is the tensor prod- ½A; B :¼ AB BA: (9)
uct of the Hilbert spaces of the subsystems.
• The Hilbert space dimension of the composite system is The usual statement of Newton’s laws53,54 discusses the kind of equa-
the product of the dimensions of the subsystems. tions of motion that are allowed rather than the underlying mathemat-
• Entangled states are states that cannot be represented as ical assumptions. This does not make it easy to compare the
tensor products of the states of the individual subsys- underlying mathematical assumptions of quantum and classical
tems; i.e., states that are a nontrivial superposition of mechanics. Perhaps, this is one of the reasons that quantum mechan-
the states of the subsystems. ics took some time to be developed and understood. If one were to
2. Measurement of physical observables. attempt to write the postulates of classical mechanics in a form that
(a) Physical observables are self-adjoint operators O† ¼ O on parallels that of quantum mechanics, then one might arrive at the
Hilbert space. following.
(b) Measurements yield eigenvalues of the operator corre-
sponding to the observable 2. Postulates of classical mechanics (CM)
Ojwa i ¼ ajwa i: (2) In the following formulas, the classical phase space will be treated as if
it were continuous, but if it is discrete, then the integrals should be
(c) The probability of observing a particular eigenvalue, a, is
replaced by sums.
given by the overlap between the eigenstate and the physi-
cal state 1. Physical states.
X (a) Pure physical states are points in a set, z 2 m, called the
ProbðaÞ ¼ hwa jqjwa i ¼ pj jjhwa jwj ijj2 : (3) phase space.
j
• Phase space is typically a continuous and smooth space
Therefore, the expectation value of measuring an observ- such as a differentiable manifold, but it can also be a dis-
able is given by crete set.
X
hOi ¼ Tr ðqOÞ ¼ pj hwj jOjwj i: (4) (b) Mixed physical states are represented by a probability dis-
j tribution function, pðzÞ, over phase space.
Ð
(d) Immediately after a measurement of eigenvalue a, the • Probability is normalized to unity pdz ¼ 1.
wavefunction collapses to the state jwa i. (c) The phase space of a composite system is the tensor prod-
uct of the phase spaces of the subsystems.
3. Temporal evolution. • The dimension of the composite phase space is the
(a) General open systems product of the dimensions of the subsystems.
• Physical states evolve via linear superoperations that act
2. Measurement of physical observables.
on the density matrix (a) Physical observables are real functions of coordinates in
qðtÞ ¼ eðt; t0 Þqðt0 Þ: (5) phase space, O : m ! R.
(b) In principle, measurements of the coordinates, z, can be
The evolution must preserve the Hermiticity, positivity, made with arbitrary precision, and, so, a measurement of
and trace of any density matrix. an observable at the point z yields its value, OðzÞ.
• Continuous time evolution is generated by the generator (c) The probability of measuring a particular value of the
of a linear superoperator observable, OðzÞ, is given by the probability, pðzÞ, of
observing the pure state at the point z. Therefore, the
dq=dt ¼ lðtÞqðtÞ: (6) expectation value of an observable is given by
The evolution must preserve the Hermiticity, positivity, ð
and trace of any density matrix. hOi ¼ Opdz: (10)
(b) Ideal closed systems
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-4

C Author(s) 2023
V
3. Temporal evolution. information has effects that are similar to coarse-graining. Thus,
open quantum systems can display many different types of non-
(a) General open systems
• Physical states evolve via (generically nonlinear) coordi-
ideal and non-unitary behavior.
In fact, the act of measurement is the primary example of a non-
nate transformations, f : m i ! m,
ideal but linear evolution process. Because the measurement apparatus
zðtÞ ¼ fðzðt0 Þ; tÞ: (11) in QM Postulate 2 is assumed to be classical, the probability of each
eigenvalue is described by classical probability theory. Thus, even if
• Continuous time evolution is generated by (generically one begins with a pure state, after the measurement is performed, the
nonlinear) differential equations resulting knowledge about the system must become a mixed state: a
probability distribution over possible eigenfunctions of the observable
dz=dt ¼ Vðz; tÞ: (12)
operator that was measured.
(b) Ideal closed systems Decoherence refers to the way in which small but persistent inter-
• Ideal closed systems are Hamiltonian systems that actions with the environment cause an almost irretrievable loss of
evolve via Hamilton’s equations of motion, which yield information for open quantum systems. The phrase decoherence was
infinitesimal symplectic transformations of phase space. introduced by Zeh55 and has been explored in great detail by
• In canonical coordinates, z ¼ fq; pg, composed of equal Zurek56,57 and collaborators. Physicists and mathematicians, such as
numbers of configuration space coordinates, q, and Choi, Gorini, Krauss, Kossakowski, Lindblad, Redfield, Stinespring,
momentum space coordinates, p, Hamilton’s equations and Sudarshan, explored possible types of evolution laws, called master
of motion for Hamiltonian HðzÞ are equations, that an open quantum system can undergo.52 The goal is to
develop a master equation that describes the loss of information in a
dq=dt ¼ @p H dp=dt ¼ @q H: (13) manner similar to the way the Fokker–Planck equation describes the
loss of information due to diffusive processes for classical probability
• Evolution over a finite time interval yields a (generically theory.
nonlinear) symplectic transformation of phase space, If one were to use an accurate quantum model of the environ-
f : m i ! m. ment, then the state of the system can interact in many different
ways with the environmental degrees of freedom. For simplicity, it
B. Open quantum dynamics is customary to attempt to determine a reduced order model that
1. Open systems and decoherence has the expected qualitative behavior. For example, it is customary
to model the interactions with the environment as a quantum pro-
Many physicists and mathematicians are familiar with the temporal evo- cess: a linear transformation of the density matrix that has no
lution of ideal closed systems, e.g., in the classical CM Postulate 3b and memory, i.e., a Markovian process. Note, however, that using a
quantum QM Postulate 3b contexts. The more general formulation for more realistic reduced order model of the environment that has
open classical systems in CM Postulate 3a is also quite natural. This form nontrivial dynamics will typically generate a nonlinear model that
applies to non-ideal open classical systems, which typically experience is not Markovian.
various forms of dissipation, such as friction, diffusion, and collisions.
In contrast, the discussion of quantum dynamics in QM 2. CPTP quantum processes
Postulate 3a was far too brief to do justice to this important topical
area. While the mathematical justification is on solid footing today, in Linear operators that act on the density matrix, rather than on
terms of applications to real-world quantum hardware platforms, this pure states, are referred to as quantum processes, quantum channels,
field is still under active development. or superoperators. If the dimension of the Hilbert space is N, then
For closed ideal quantum systems, QM Postulate 3b requires there are N 2 unitary operations that describe ideal evolution.
time evolution to be unitary and, hence, reversible. However, simi- However, there are many more, N 4 N 2 , quantum processes that
lar to the way that classical chaos can scramble classical informa- describe decoherence.
tion, quantum dynamical processes can scramble quantum If the density matrix is to remain a positive Hermitian operator,
information. Classically, the combination of the butterfly effect then the superoperator must itself be Hermitian and positive. In fact, it
(positive Lyapunov exponents) and coarse-graining due to the is common to demand a stronger criterion called completely positive
finite precision of measurements leads to irreversibility. However, (CP)58,59 that results from considering the evolution of the system
quantumly, the unitary dynamics only has vanishing Lyapunov when coupled to a reservoir that represents the environment with a
exponents, and, so, while the evolution may track the semiclassical Hilbert space of arbitrary dimension, d. In this case, the state is
dynamics for long periods of time, any apparent exponential insta- q 1d , and evolution of the subsystem with superoperator, e, while
bility associated with the Lyapunov exponent must always eventu- the reservoir is unchanged results in the state eðqÞ 1d . The map,
ally cease. eðqÞ, must be completely positive in order for the coupled map,
The paradox is resolved by the fact that most quantum sys- eðqÞ 1d , to be positive for the total system reservoir Hilbert
tems are actually open; i.e., the observer does not have access to space. It turns out that complete positivity is a stronger condition than
information to all parts of the system. The parts of the system that positivity alone as there are examples of maps, e.g., the transpose oper-
either are not or cannot be measured are usually referred to as the ation, that are positive but not completely positive. In order to con-
environment. The existence of this unmeasured or hidden serve probability, the superoperator must also be trace-preserving
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-5

C Author(s) 2023
V
(TP). Such superoperators are called completely positive trace-preserv- including the unit matrix, which possesses a trace, was already used in
ing (CPTP) operators. the construction of the damping operator A.) This evolution is CPTP
A Hermiticity preserving superoperator must have the form52,59 as long as the rate matrix, ij, is Hermitian and positive semidefinite.
N2
Information loss occurs whenever the rate matrix has positive
X
eq ¼ vjk Lj qL†k ; (14) eigenvalues.
j;k¼1 Again, because the rate matrix is Hermitian, this form can be
simplified by diagonalizing jk with a unitary transformation.
where fLj g defines a complete basis of operators of size N N and the Applying this transformation to the collapse operators allows one to
matrix vij is Hermitian. In order for the map to be CP, the matrix vij
find a new basis, L^ a , that simplifies the equation to a sum over the
must also be positive semi-definite. (Note that completely positive refers
eigenvalues of the vij matrix, which we denote as a
to the action on the total Hilbert space rather than the definiteness of
the eigenvalues.) In order to be trace preserving (TP), the condition X †
n
†
o
@t q ¼ ½H; q=ih þ a L^ a qL^ a L^ a L^ a ; q =2 : (19)
N2
X a
1¼ vjk L†k Lj (15)
Since the eigenrates a are positive, the square root of the rates can be
j;k¼1
absorbed into the definition of the collapse operators, as in the Krauss
must hold. The trace preservation criterion represents N 2 constraints. form above. While the diagonal form offers simplicity, it hides the
Since there are N 2 unitary operations, this implies that there are N 4 complexity associated with the myriad ways in which the environment
2N 2 non-unitary decoherence operations. can potentially interact with a quantum system to form an arbitrary
Because the vjk matrix is Hermitian, this form can be simplified rate matrix. A pedagogical introduction to the GKLS master equation
by diagonalizing this matrix with a unitary transformation. Applying is given in Ref. 63.
this transformation
to the operator basis allows one to find a new There are a few important environmentally induced decoherence
basis, L^ a , that simplifies the equation to a sum over the eigenvalues processes that affect almost all quantum systems: (1) relaxation: of
of the vjk matrix, which we denote as va higher energy levels to lower energy levels, (2) excitation: of lower
energy levels to higher energy levels, and (3) dephasing: the diffusion
N2
X
eq ¼
†
va L^ a qL^ a : (16) of the complex phases within a linear superposition. For a harmonic
a¼1
oscillator, the relaxation collapse operator is the destruction operator,
a, the excitation collapse operator is the creation operator, a† , and the
Because the eigenvalues va are positive, the so-called Krauss form60 dephasing collapse operator is the number operator, N ¼ a† a. While
can be achieved by absorbing the square root of the eigenvalue into excitation is important in general contexts it is often much slower than
1=2
the definition of the operators L^ a ! va L^ a . the other two processes for typical quantum hardware platforms.
Note that completely positive only requires the eigenvalues to be
positive semi-definite va 0. If the map is positive but not completely
C. The qubit
positive, then it must have at least one negative eigenvalue.
Just as the bit is the smallest unit of classical information, the
3. CPTP evolution qubit is the smallest unit of quantum information. Both are shown
schematically in Fig. 1. The bit can only ever be in one of two classical
The Gorini–Kossakowsi–Lindblad–Sudarshan (GKLS) master states: 0 or 1. In the language of quantum information, we might call
equation61,62 represents the most general form of a continuous CPTP these classical states j0i and j1i. However, it would be clearer to label
evolution of the density matrix when the system is separable from the these states as density matrices j0ih0j and j1ih1j to signify the fact
environment and when their interactions are Markovian. Just as the that, in the classical world, superpositions are not allowed and, hence,
Fokker–Planck equation describes processes such as diffusion due to these two states are mutually exclusive. In other words, in the classical
collisions for the Liouville evolution of the classical PDF, the Lindblad world, only one of the two possible pure states j0ih0j and j1ih1j is
equation adds nonunitary CPTP processes to the von Neumann equa- allowed to exist at any given moment in time.
tion for the density matrix. The GKLS equation is given by To make a physical bit from a physical device, one must devise a
2
physical system that is subject to a potential that has two deep wells, as
NX 1
shown in Fig. 1. Then, one must ensure that the temperature of the
@t q ¼ ½H; q=ih fA; qg þ jk Lj qL†k ; (17)
j;k¼1
system and, more generally, the possible range of energy exchanges
with the environment remains far below the energy barrier separating
where the anticommutator is defined via fA; Bg :¼ AB þ BA. Thus, the two wells.
in addition to unitary evolution determined by the Hamiltonian H, If one constructs a similar potential with a quantum mechanical
one adds a general CP evolution as well as the TP damping term device, then the wavefunction can potentially be in a superposition of
2
all accessible states, including the lowest energy level, called the ground
X
N 1
state or vacuum state, as well as the higher energy levels, called the
A¼ jk L†k Lj =2 (18)
excited states, also shown in Fig. 1. To make a physical qubit from a
j;k¼1
physical device, one must confine the wave function to two states. For
required to preserve the trace. In this case, the collapse operators, fLj g, this two-state system, j0i is the ground state and j1i is the first excited
define a basis over trace-free matrices of size N N. (The freedom of state. Again, one must control the temperature and the interactions
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-6

C Author(s) 2023
V
Given the definition above, when measurements are made along

the ^z axis, the probability of j0ih0j is p0 ¼ cos2 ðh=2Þ ¼ ð1 þ cos hÞ=2
and the probability of j1ih1j is p1 ¼ sin2 ðh=2Þ ¼ ð1 cos hÞ=2. In
fact, the probability of measuring a given outcome is simply given by
the binomial probability distribution with probability p1 of obtaining
the value 1 on each trial. For S trials, the binomial distribution has
mean p1 S and variance p0 p1 S, so the variance of the mean decreases as
the central limit theorem would predict, p0 p1 =S. In agreement with
QM Postulate 2: immediately after a measurement, the state is known
with certainty, so the variance must vanish. As S ! 1 for fixed prob-
ability, the binomial distribution is well-approximated by a normal
Gaussian distribution, but this approximation breaks down for small
sample sizes or if one of the probabilities becomes vanishingly small.
For example, if S ! 1, but the product k ¼ p1 S is held fixed, then
the binomial distribution tends toward the Poisson distribution with
mean and variance k.
There is also the possibility of rotating the measurement basis to
point along the x^ or y^ directions. In this case, the probability is con-
trolled by the angle / rather than h. In practice, it usually easier to
keep the measurement apparatus fixed and, instead, to rotate the state.
For example, to measure the probabilities along the 6^ x axis, simply
rotate the state counterclockwise around the y^ axis by 790 .
FIG. 1. A classical bit (top left) can only ever be found in the states j0ih0j or j1ih1j It is useful to explicitly define the effect of a rotation as a unitary
(bottom left). Physically, this is created with two deep potential wells that trap the operation on the qubit. Define the Pauli vector as the vector of Pauli
physical state in one region or the other. A quantum qubit (top right) can be in any matrices ~ r ¼ ðrx ; ry ; rz ÞT . A rotation around the unit vector n
^ by the
superposition of states j0i and j1i (bottom right). Physically, this is created using
the ground and first excited states of a quantum system. angle h is defined by
Rn^ ðhÞ ¼ eih^n ~r =2 ¼ cos ðh=2Þ1 i sin ðh=2Þ^

n ~
r: (23)
with the environment so that the probable range of energy exchanges Note that a rotation by h ¼ 2p multiplies the qubit by the overall
with the environment is much less than the energy differences between phase eip ¼ 1. This explicitly demonstrates the fact that the special
these two states as well as the energy differences between these two unitary group, SU(2), is actually a double cover of the special orthogo-
states and all other states. Moreover, one must also be careful to con- nal group of rotations SO(3).
trol environmental interactions that cause decoherence and destroy A classical probability distribution over qubit wavefunctions is
superposition and entanglement. equivalent to a positive Hermitian matrix called the density matrix
The wavefunction of a qubit lives in a two-dimensional complex
Hilbert space q00 j0ih0j þ q01 j0ih1j
q¼ : (24)
i/=2 i/=2
þq10 j1ih0j þ q11 j1ih1j
jwi ¼ sin ðh=2Þe j0i þ cos ðh=2Þe j1i: (20)
Due to Hermiticity, the off diagonal elements must be complex conju-
Due to the fact that the wavefunction is a ray in Hilbert space, i.e., it gates, q10 ¼ q01 and the diagonal elements must be real. These diago-
actually lives in projective Hilbert space, the overall complex constant nal elements represent the probability of measuring the two states,
is unimportant and is normalized to unity. Note, however, that for j0ih0j and j1ih1j, and are guaranteed to be positive due to the positiv-
describing a normalized vector in Hilbert space, and for describing ity of the density matrix. In order for the probabilities to sum to one,
unitary transformations of physical states, one must include the overall there is the constraint Tr q ¼ q00 þ q11 ¼ 1. The off diagonal ele-
complex phase factor, eiv jwi. ments represent non-classical superposition states that have no analog
Thus, the qubit is described by two real parameters that form the in the classical world.
surface of a sphere, called the Bloch sphere. By convention, the þ^z Again, due to Hermiticity, the density matrix can be diagonal-
direction that points toward the top of the sphere corresponds to the ized, so that it has the form
j0i state while the ^z direction that points toward the bottom of the
sphere corresponds to the j1i state. Thus, an arbitrary qubit state can q ¼ pþ j^
n ih^
n j þ p j^
n ih^
n j; (25)
also be defined by the direction that it points in; if the unit direction
^ , then the state is denoted j^
vector is n n i. For example, j^z i :¼ j0i and where j6^ n i are orthogonal states and p6 are the classical probabilities
j^z i :¼ j1i, but remember that these states are orthogonal. Similarly, of each of these states occurring. Any 2 2 density matrix can also be
one can define the states written in the form
q ¼ ½1 þ ðpþ p Þ^
n ~
r =2: (26)
x i ¼ ðj0i 6 j1iÞ=21=2 ;
j6^ (21)
y i ¼ ðj0i 6 ij1iÞ=21=2 :
j6^ (22) Thus, the qubit density matrix can be visualized as an arrow of length
jpþ p j pointing in the n
^ direction, called the Bloch vector. A proper
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-7

C Author(s) 2023
V
mixed state must have jpþ p j < 1. A pure state can only have numbers. Hence, the quantum wavefunction holds exactly twice the
eigenvalue 1 for one entry and eigenvalue 0 for the other, so that information of the classical PDF due to the fact it has a complex phase
Tr qa ¼ 1 for any a > 0. In contrast, for a proper mixed state, both associated with each probability amplitude. Yet, the classical PDF
eigenvalues must be less than 1, so that Tr ðqa Þ < 1 for any power should really be compared to the quantum density matrix, which is
a > 1 [and Tr ðqa Þ > 1 for any power 0 < a < 1]. specified by N 2 1 ¼ 22n 1 real numbers and, thus, has quadrati-
Given the discussion above, a classical probability distribution cally more information.
function (PDF) over classical bits is simply a diagonal density matrix In many ways, these simple facts are the reason for all of the
excitement about quantum computing. This capacity of a quantum
p ¼ Diag ðqÞ ¼ q00 j0ih0j þ q11 j1ih1j: (27) memory register for holding an exponentially large amount of infor-
Because it is diagonal, it can also simply be considered to be a function mation and the potential ability of a quantum computer to efficiently
of the classical states alone, p : f0; 1g ! R, with the two values manipulate this exponentially large amount of information is amazing
pð0Þ ¼ p0 ¼ q00 and pð1Þ ¼ p1 ¼ q11 . Again the normalization con- from the standpoint of classical computation.
dition, Tr ðpÞ ¼ 1, ensures that the probabilities sum to unity. This also implies that performing a direct classical simulation of
a multi-qubit quantum system is exponentially hard! Although many
D. Multiple qubits ingenious algorithms have been invented for approximately solving
these kinds of problems, progress with direct solvers is limited due to
In the following, we assume that individual qubit states can be the exponentially large size of the required Hilbert space. Classical
physically prepared and measured in the same basis each time, which computations much use drastically reduced models to even attempt
is referred to as the computational basis. Two qubits live within a four- the simulation of large systems. Now we can understand the first part
dimensional Hilbert space, which is the tensor product of the two indi-
of Feynman’s statement: If you want make a (quantum-mechanical)
vidual qubit Hilbert spaces. Hence, basis states that span the two qubit
simulation, you’d better make (the computational engine) quantum
Hilbert space can be denoted jq1 q2 i ¼ jq1 i jq2 i, where q1;2 refer to
mechanical… so that it has a sufficiently large Hilbert space and so
the state of qubit 1 and 2, respectively. There are four possible mea-
that this Hilbert space can be manipulated efficiently.
surement outcomes, which correspond to four possible classical bit
Note that classical algorithms that can efficiently approximate
strings. However, there are many more quantum states that are not
the classical PDF are also quite powerful. From the simple scaling of
classical, corresponding to three complex degrees of freedom.
the amount of information in the classical PDF vs the quantum den-
In addition to single-qubit superpositions, the key resource to be
sity matrix, one might expect that quantum algorithms can generically
exploited is the ability to create entangled states. Entangled states are
outperform classical probabilistic algorithms by this same quadratic
defined as states that cannot be represented as tensor products of the
factor. As will be discussed in detail later, this is indeed the case.37,64
states of the individual subsystems. Thus, a state such as j^ x xî
¼ ðj00i þ j01i þ j10i þ j11iÞ=2 is not entangled, even though each
qubit is in a superposition state itself. In contrast, a two-qubit state E. Qudits
such as the Bell state ðj00i þ j11iÞ=21=2 is truly entangled. When Just as classical information with d possible values is referred to
quantum information is entangled, but only part of the system can be as a dit, the quantum analog of a wavefunction with d possible states is
measured, this leads to kind of information hiding that causes deco- called a qudit. For both the classical and quantum cases, systems with
herence. Note that the concept of entanglement explicitly refers to multiple bits and dits can be used to simulate each other with polyno-
consideration of the tensor product structure of the Hilbert space and, mial efficiency, so there is no change in computational complexity
hence, the basis that is being used. class. In the classical world, transistors are so cheap that essentially all
In order to find the size of a multi-qubit wavefunction, simply modern computing technology is based on bits. In the quantum world,
apply QM Postulate 1c: multiply the dimensions of the Hilbert spaces experimentalists are still striving to reduce the error budget to the
of the subsystems and apply the ray condition. As illustrated in Fig. 2, point at which error correction can be applied. For this reason, a num-
for n qubits, there are N ¼ 2n computational basis states, which is ber of researchers have proposed using the higher energy levels of a
equivalent to the number of possible measurement outcomes, and, quantum system to encode information, as illustrated in Fig. 3. For
hence, to the number of possible classical bit strings. For n qubits, the certain quantum device platforms, these higher energy states are natu-
wavefunction is specified by N ¼ 2n 1 complex numbers. Although rally present, so it could be a good idea to effectively pack more qubits
classically, the state is only ever in one of the possibilities, a classical into a single device, potentially allowing less error-prone communica-
probability distribution (PDF) of n bits is specified by N ¼ 2n 1 real tion between these states.65 However, since decoherence rates are
FIG. 2. The Hilbert space of multiple

qubits is the tensor product of the Hilbert
spaces of each qubit in the system. The
total composite Hilbert space dimension is
2n where n is the number of qubits.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-8

C Author(s) 2023
V
required to generate the whole group. Another important example is

the Hadamard gate

1 1 1
H :¼ ZRYðp=2Þ ¼ 1=2 : (28)
2 1 1
Thus, starting from the j^z i :¼ j0i state pointing along the ^z axis, this
generates the uniform superposition state j^ x i. Similarly, starting from
the state j^z i ¼ j1i, this generates the superposition state j^ x i.
Adding just one nontrivial two-qubit entangling gate between
nearest neighbor qubits is enough to generate all other operations. The
two-qubit Hilbert space is a four-state system and these operations
have the form of a 4 4 matrix. Common choices are CNOT ¼ CX
and CZ
0 1 0 1
1 0 1
B0 1 C B 1 C
CNOT :¼ B @
C CZ :¼ B C: (29)
FIG. 3. Qudits use more than two states, e.g., j0i; j1i; j2i, ….
0 1A @ 1 A
1 0 1
Another useful operation is the SWAP operation that simply
often faster for higher energy levels, there is clearly a trade-off in per- exchanges information between qubits
formance that must be investigated for each experimental hardware 0 1
1
platform. B 0 1 C
SWAP :¼ B @
C:
A (30)
III. QUANTUM COMPUTATION 1 0
A. Models of universal quantum computation 1
Universal quantum computation5–7,52 refers to models of quan- If U22 is a single-qubit unitary and 122 is the single-qubit unit
tum computing that can efficiently simulate a universal classical com- matrix, then a controlled U operation has the form
puter, i.e., a Turing machine. There are a number of different models
of universal quantum computation that have been proven to be equiv- 122
CU :¼ : (31)
alent to one another in computational capabilities; i.e., they can simu- U22
late each other with polynomial complexity. Popular models include This shows that the operation U is performed if the first qubit is in the
the digital circuit model,5–7 the quantum walk model,66,67 the adiabatic j1i state. To make this clearer, use the notation CUjk to indicate the
model,68 the quantum annealing model,69,70 and the measurement- fact that the state of the kth qubit is controlled on the state of the jth
based model.71 The digital circuit model, also known as the digital qubit.
quantum computational model, has power and simplicity and is per- A sequence of unitary operations can be visualized through a cir-
haps closest to digital classical computation based on classical circuits, cuit diagram. Operations are applied from left to right in a circuit
which is familiar to most scientists and engineers today. diagram—which is opposite to the usual right to left operator ordering
The basic idea of the digital circuit model is that quantum states convention in mathematical expressions. A circuit diagram has a num-
can be transformed efficiently via linear unitary operations: w ¼ Uw0 . ber of lines going from left to right that each indicate individual qubits.
In fact, according to QM Postulate 3b, this is the only operation that Blocks in the middle of the diagram represent single-qubit operations.
can be performed without making a measurement. After making a Blocks that are controlled on another qubit have a line from the block
measurement, the state generically becomes a density matrix q, but in to the control qubit. For example, consider the creation of the
between measurements, linear unitary operations can be performed: entangled Greenberger–Horne–Zeilinger (GHZ) state72 on three
q ¼ Uq0 U† . The fact that a quantum computer can do this efficiently qubits defined by ðj000i þ j111iÞ=21=2 (an analog of the Bell state
is unbelievable from the classical standpoint, because it allows one to for two qubits). This state can be created with the circuit diagram
efficiently perform matrix–vector operations with exponentially large shown in Fig. 4. The instructions are to (i) apply H to the first
matrices and exponentially large vectors. qubit; (ii) perform a CNOT1;2 operation, where the state of the sec-
While there are a large number, N 2 , of independent unitary oper- ond qubit is controlled by the state of the first qubit; (iii) perform a
ations, they are generated by a small, OðnÞ Oðn2 Þ; number of CNOT1;3 operation, where the state of the third qubit is controlled
basic operations called a gate set. Single qubit operations can be per- by the state of the first qubit; and (iv) finally, the quantum state
formed efficiently with just two standard gates that perform rotations can be verified by measuring all qubits. If one rotates the states
of the Bloch sphere. Rotations by angle h around the x^; y^, and ^z axes before measurement, then one can test the consequences of Bell’s
of the Bloch sphere, are denoted RXðhÞ; RYðhÞ, and RZðhÞ, respec- theorem.73
tively. A rotation by angle h ¼ p, simply generates the Pauli opera- The fundamental discovery of early research on quantum
tions, denoted X ¼ rx ¼ NOT; Y ¼ ry , and Z ¼ rz . A standard algorithms is that many useful computations can be performed
result in Lie group theory says that only two such operators are using relatively few Oðpoly ðnÞÞ of these basic gate operations.5–7
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-9

C Author(s) 2023
V
• Directly sampling a quantum state to estimate the probability

amplitude costs the same number of measurements as it does
classically.
As explained previously, the probability distribution for measure-
ments of the state of a qubit is given by the classical binomial dis-
tribution. For probabilities far from 0 or 1, there is no escape
from the central limit theorem: if the measurement outcome has
a finite mean and finite variance, theprelative
ffiffiffi accuracy with which
it can be measured converges as 1= S where S is the number of
samples.
• Directly measuring exponentially small probability amplitudes is
FIG. 4. An example of a circuit diagram for the creation of the entangled GHZ state
ðj000i þ j111iÞ=21=2 . Apply H to the first qubit and CNOT gates between the sec- exponentially hard.6
ond and third qubits and the first qubit. Finally, measure the result. Since these values would only be observed with exponentially
small probability, one must make an exponentially large number
Key resources are the ability to create superpositions and entangle- of measurements to observe the values with any accuracy.
ment and to use interference to one’s advantage. This is often • Nonzero failure probability.7
called quantum parallelism. The second fundamental discovery is Quantum algorithms are probabilistic and, as such, they have
that any reversible classical computation can also be performed.5–7 a certain probability of failure. This failure probability can be
This requires first computing and storing the result and then made quite small by batching many runs together (see the
uncomputing or reversing the order of operations on the original clear explanation of the Powering Lemma in Ref. 37).
memory register.74 However, generically, the failure probability does not vanish
However, emulating a classical computer usually comes with and this feature is quite different from classical deterministic
non-negligible overhead, and so, one might prefer to use a classical algorithms.
computer to perform certain operations such as classical logic. For Thus, it is misleading to say that the universal quantum compu-
this reason, the designers of quantum computer platforms usually tation model offers a method of efficiently simulating any quantum
consider a hybrid architecture where classical processing units system. In fact, in order to remain efficient, this model can only
(CPUs) sit beside as well as drive quantum processing units approximately simulate certain quantum systems and it can never
(QPUs) much like today’s heterogenous CPU and graphical proc- actually exploit the entire Hilbert space. Although we began with the
essing unit (GPU) architectures. Like CPU/GPU platforms, there notion that there is an exponential amount of information that can be
will be a large overhead in input/output (I/O) operations that stored and manipulated, we now see that the actual amount of such
exchange data between CPUs and QPUs. In fact, the I/O will be information that can be manipulated and measured is much smaller
exponentially large if one attempts to exchange the entire quantum than one might have hoped.
random access memory (QRAM). In today’s world, QPUs are also In fact, it is useful to compare the power of the quantum comput-
far more expensive than CPU’s so hybrid architectures are an obvi- ing model to classical probabilistic models of computation, such as
ous choice. While it may be difficult to envision today, perhaps at Monte Carlo methods. For probabilistic models, one never manipulates
some point in the future, quantum computers will become useful all of the information required by a direct calculation; rather, one relies
enough and cheap enough that they will become ubiquitous. on the central limit theorem to provide a slow but steady increase in rel-
ative accuracy as the number of individual samples is increased.
B. Limitations of the quantum computing model The quantum computational model is certainly a probabilistic
model and many runs of the same algorithm are needed to yield a
Now it is very important to discuss some key limitations of the
final answer. After the required measurements are made for each
universal quantum computing paradigm.5–7
run, one can associate the outcome with a quantum trajectory that
• Approximating an arbitrary unitary is exponentially hard. is consistent with the outcome.75 It is only by assembling a large
This implies that only certain classes of unitaries, precisely those number of quantum trajectories that an accurate answer is
that are composed of Oðpoly Þ basic gate operations, can be con- achieved. In a certain sense, the key difference is the fact that quan-
structed “easily.” tum trajectories can exploit superposition states and interference
• Initializing all quantum information is exponentially hard. effects rather than classical trajectories which are limited to com-
Initializing quantum information requires applying an arbitrary putational basis states alone.
unitary operation (usually starting from the ground state). So, Nonetheless, there are many elegant mathematical theorems that
this also implies that initializing all information in quantum prove that the quantum computational model can exceed even the
memory is exponentially hard. classical probabilistic counterpart.5–7 Generically speaking, the ability
• Measuring all quantum information is exponentially hard.6 of the quantum density matrix to store quadratically more information
Measuring all quantum information requires making an expo- than the classical PDF enables up to a quadratic speedup over classical
nential number of measurements. Moreover, each measurement probabilistic methods.76 However, there are certain classes of problems
generically requires a different unitary transformation to be for which superpolynomial and even exponential speedup has been
applied. proven to hold.77
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-10

C Author(s) 2023
V
C. Quantum error correction: Improving necessary to represent a logical qubit is estimated81 to be on the order
the information confinement time of 103 104 . Hence, the number of qubits required to perform a use-
From the inception of quantum computing, it was recognized ful calculation is estimated to be on the order of 105 106 or more.84
that the ability to correct quantum errors is indispensable for building Unfortunately, this large overhead cost has been investigated and
large-scale quantum computers. Quantum error correction (QEC) pro- interpreted91 as implying that the first successful error-corrected quan-
tocols5–7 seek to undo errors that are induced by the inability to enact tum applications should offer speedups that are beyond quadratic, e.g.,
the desired unitary operations with perfect precision as well as the at least quartic.
uncontrollable interactions with the environment that cause decoher- After decades of theoretical development of error correction
ence. The general strategy is to use an encoding that repeats informa- codes and experimental demonstrations of primitive steps, a surface
tion so that it can be checked against various types of error syndromes. code was recently demonstrated in a superconducting circuit,92 and a
This implies that there is a distinction between the number of physical color code was recently demonstrated on a QCCD ion-trap quantum
qubits that comprise the quantum hardware platform and the number computer.93 These recent demonstrations are starting to approach the
of logical qubits that can encode quantum information for a given time threshold where logical qubits outperform physical qubits.
period with negligible error rates. The subject of QEC is vast and we Since each operation takes a certain physical time period, the
cannot do justice to the wide range of ideas necessary to understand limit on the error rate can also be translated to a limit on the norm of
the contemporary approaches to the field. The interested reader can the effective decoherence rate matrix, jjjj. Thus, the basic idea of
begin by consulting standard textbooks,5–7 books on quantum error quantum error correction is to increase the information confinement
correction,78–80 and relatively recent review articles.81–85 time until this limit can be achieved. Perhaps surprisingly, this goal is
In general, quantum error correction protocols seek to undo very similar to the Lawson criterion, which requires a net energy-
errors that are induced by the inability to enact the desired unitary producing fusion reactor to achieve a specific energy confinement
operations with perfect precision as well as the uncontrollable interac- time. Within the field of fusion energy, in addition to inherent physical
tions with the environment that cause decoherence. As discussed in limitations of the reactor design, it is anomalous and often turbulent
Sec. II B, if the system of interest has N states, then, each time one transport processes that typically control the energy confinement time.
attempts to enact a quantum gate, there are N 4 N 2 ways, including Similarly, within the field of quantum computing, in addition to inher-
both decoherence operations and accidental unitary operations, that ent physical limitations of the hardware platform, it is decoherence
the intended gate can go awry. If one is attempting to obtain an expo- due to anomalous interactions with the environment that controls the
nential speedup by utilizing a large number of qubits, so that N ¼ 2n , information confinement time. Thus, it will be interesting to gauge the
then there are an exponential number, N 4 ¼ 24n of gate features relative progress in the two fields by measuring the rate of improve-
that must be controlled. Just obtaining enough measurements to fully ment of the respective confinement times.
characterize the operation of the gate then becomes an intractable
problem. Now we can understand the second part of Feynman’s quote: D. A brief history of quantum algorithms
by golly it’s a wonderful problem because it doesn’t look so easy. The field of quantum algorithms has seen a rapid development in
For example, from the point of view of gate design, it is very effi- interest over the last two decades and it is not possible for a short
cient to have all-to-all connectivity between qubits, as is the case for introduction such as this to do justice to the wide range of brilliant
ion trap platforms.86 This eliminates the need for many extra SWAP ideas that have begun to bear fruit. Nonetheless, when compared to
gates that would otherwise be needed to shuttle quantum information classical algorithms, the number of essentially quantum subroutines
between different parts of the memory register. Because these SWAP are actually relatively few. In part, it is the vast number of ways that
gates must often be performed in serial, eliminating them also these subroutines can be configured as well as the number of applica-
increases the potential for parallelization of multiple gate operations. tions that these quantum subroutines can speed up that leads to the
However, from the point of view of error correction, it is far supe- considerable complexity in the field. Moreover, since the field is still
rior for qubits to be well-separated both in terms of their resonant fre- undergoing rapid development, there may be many new concepts wait-
quencies and physical locations. For example, this can be ing to be discovered.
accomplished by carefully choosing and tuning the resonant frequen- The challenge of performing computations using an exponen-
cies and by physically arranging qubits into a lattice, as is the case for tially large quantum Hilbert space was reconceived as a potential
superconducting platforms.87 Another example is the quantum charge opportunity by Feynman50,51 and Manin49 in the early 1980s.
coupled device (QCCD) ion-trap architecture,88 which manages its Feynman reported that his interest in the subject was piqued by
error budget by moving inactive ions to distant locations from active Charles Bennett, who had worked to understand reversible classical
gate regions. If one can limit the amount of significant crosstalk computation74 and who wanted to understand how the physics would
between qubits to m nearest neighbors alone, then the number of gate work when using the laws of quantum mechanics.51 In fact, it was
features that need to be controlled is reduced to 24m per qubit. The key Charles Bennett and Giles Brassard that made the first inroads into
point is that now m is independent of the size of the memory register. proving rigorous theorems that demonstrated that quantum commu-
Known quantum error correction algorithms can be proven to be nication channels had more information processing capacity than clas-
successful, allowing one to store and manipulate quantum information sical communication channels.6,7
indefinitely, as long as the error rates are below a specific limit depend- In the 1980s, Benioff48 and Deutsch94 defined quantum Turing
ing on the protocol, typically at 103 105 per operation. For the machines and showed that that they were as powerful as classical
best performing approaches known today, e.g., the topological surface Turing machines. In 1992, the Deutsch–Josza algorithm95 was discov-
code89 with magic state distillation,90 the number of physical qubits ered, a concrete example of a toy problem that could be solved faster
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-11

C Author(s) 2023
V
with a quantum computer than a classical one. Soon after, the Eventually, general-purpose Hamiltonian simulation algorithms
Bernstein–Vazirani algorithm96 and Simon’s algorithm,97 were discov- (HSAs) that exploit sparsity in the Hamiltonian itself were devel-
ered as well. It turns out that these simple algorithms rely on the ability oped.115 These algorithms were then improved by using quantum
to perform the Fourier transform more efficiently on a quantum walks,116,117 new representations of Hamiltonians as a linear combina-
computer. tion of unitaries118,119 (LCoU), and as unitary approximations of the
In 1994, Peter Shor’s algorithm98 for factoring integers and taking truncated Taylor series of the propagator.120 The latest development of
discrete logarithms exponentially faster than a classical computer sur- these algorithms eventually achieved nearly optimal dependence on all
prised the world with the realization that quantum computers can parameters.21,22,121,122
break classical security protocols relatively easily. The key step is esti- Hamiltonian simulation powers the Harrow, Hassidim, and
mating the periodicity of a sequence with unknown period and the Lloyd (HHL) sparse quantum linear solver algorithm123 (QLSA) as
key discovery is that there is an exponentially efficient factorization of well as improved versions with superior performance.124–127 Next,
the discrete Fourier transform on quantum memory registers—today sparse quantum linear differential equation solver algorithms (QLDSA)
called the quantum Fourier transform (QFT). After using the QFT, the were developed that use linear solvers in combination with sparse
period can be estimated using phase estimation protocol to estimate Hamiltonian simulation algorithms.18,128
the eigenvalue of a Hamiltonian, given an eigenvector. Shor then went Today, there are new classes of optimal Hamiltonian simulation
on to proving that there were quantum error correction protocols that methods powered by quantum signal processing (QSP)22 and qubitiza-
could be used to protect quantum information,99 which gave the tion.129 The term qubitization refers to the use of an ancillary qubit, as
abstract theory a sense of physical realizability, perhaps for the first in quantum walks, while the term quantum signal processing refers to
time. the way in which this ancillary qubit is manipulated to perform useful
The next breakthrough was Lov Grover’s search algorithm100,101 calculations. The combination of these methods can generate nonlin-
which speeds up unstructured search by the square root of the number ear quantum eigenvalue transformations (QEVT) of a block-encoded
of items. Key elements of the search algorithm were then separated operator that is a normal matrix, ½H; H† ¼ 0. This approach was then
into algorithms for amplitude amplification, quantum counting, and generalized to generate nonlinear quantum singular value transforma-
amplitude estimation.102–104 These subroutines form the basis for tions (QSVT) for an arbitrary matrix.126 The use of QEVT and QSVT
quantum algorithms for computing sums and integrals64,105–107 and was shown to provide a grand unification of many quantum algo-
forms the basis for quadratically speeding up classical Monte Carlo rithms,122,126 including search, amplitude amplification. matrix multi-
algorithms.37 plication, matrix inversion, phase estimation, and QFT. QSP has
In the late 1990s, Seth Lloyd and Daniel Abrams discovered already been put to good use within quantum algorithms designed to
efficient algorithms for simulating klocal Hamiltonian interac-
solve plasma physics problems.23,130
tions.108,109 The basic idea is to split the Hamiltonian into a rela-
tively low number of terms that do not commute with each other
and to then use the Lie–Trotter–Suzuki decomposition, aka E. Key quantum subroutines
Trotterization, to approximate the exact unitary as a product of There are a number of useful textbooks5–7 and review
terms. For example, for Hamiltonians defined on a 1D lattice with articles131–133 on quantum algorithms. Although the quantum compu-
nearest neighbor interactions, one can split the Hamiltonian into tational model might take some time to become familiar with and new
a sum of three pieces, wherein all terms within each piece com- quantum algorithms with improved complexity are under rapid devel-
mute: a diagonal piece, an even index off diagonal piece, and an opment, it is helpful to realize that there are actually only a few essen-
odd index off diagonal piece. As another example, the single par- tially quantum subroutines that power the vast majority of quantum
ticle Schr€odinger equation can be simulated efficiently110 in a algorithms.
manner that mimics the phase space path integral by splitting the
kinetic and potential energy terms into two pieces. Then one can 1. Quantum Fourier transform: Cost of O(log2N) rather
repeatedly use the QFT and inverse QFT to switch back and forth than classical N log N
between the momentum and position basis, so that each term is
diagonal in the respective basis. While operator splitting meth- In the 1990s, it was discovered that the Fourier transform can
ods, such as Strang splitting and higher order Suzuki methods,111 be factorized in an extremely efficient manner on a quantum mem-
are familiar to applied mathematicians, recent work112 has led to ory register. Many Hamiltonians of interest can be simplified and,
new understanding of error convergence for these methods. hence simulated efficiently, if the Fourier transform is exponen-
Quantum walks113 are a vast generalization of the amplitude tially cheap and if one does not need to extract intermediate infor-
amplification paradigm that allow discrete state changes based upon mation. In fact, the key subroutine that powers Shor’s algorithm
the outcome of flipping a quantum coin or rolling quantum dice; i.e., for factoring integers and taking discrete logarithms is the Fourier
on the essential use of ancillary qubit registers. Quantum walks were transform.
originally constructed for the purpose of understanding simple models There are many variations of the Fourier transform that
of quantum physics. Notably, Richard Feynman developed a quantum arise from the character theory of commutative groups. The
walk model for how Dirac fermions traverse a discrete lattice.114 By most familiar of these is the Walsh–Hadamard transform, which
the early 2000s, the essential ideas behind quantum walks were codi- arises from the group ðZ=2ZÞn . The Walsh–Hadamard trans-
fied into a powerful set of algorithms113 and were then proven to offer form on multiple qubits, Hn , is simply the tensor product of
a universal computing model of their own.66,67 Hadamard gates
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-12

C Author(s) 2023
V
X
Hn ¼ nk¼1 H: (32) QFT ¼ nk¼1 Rðbxck hÞH j0ihjj (44)
j
The subscript n is not necessary, because the number of qubits to be X
acted upon should be clear from the context. Evidently, this transform ¼ Rðbxc1 hÞH Rðbxc2 hÞH … ½Hj0ihjj: (45)
factorizes so well that it can be performed with single-qubit rotations j
alone. It is simple to show by induction that this generates a uniform
Now, it is easy to see that each rotation should itself factorize as the
superposition of states from vacuum
product of n k 1 controlled rotations, CR‘;k , between the kth
X
N qubit and all higher ‘ th qubits. The operator
jsi :¼ Hj0i ¼ jxi=N 1=2 : (33)
x¼1
^ n
Rðbxck hÞ ¼ ‘¼kþ1 CR‘;k ðx‘ hÞ (46)
Then, a quantum subroutine can work on the entire superposition at
once: w ¼ Ujsi. This simple initialization step is used by many quan-
tum algorithms. is almost correct, but it applies the phase to the kth qubit rather than
The discrete Fourier transform arises from the modular group the n–kth qubit and the desired rotations were supposed to be imple-
Z=NZ and is defined via mented in the opposite order within the tensor product. This is easy to
X correct by first computing the result for both k and n–k and then
hkjQFTjwi ¼ N 1=2 wj e2pijk=N : (34) swapping the results. In other words, simply add a final stage of
j SWAP gates at the end defined by
For the binary version of the algorithm, N ¼ 2n , and, one can express
bn=2c
the index in binary notation via SWAPbn=2c :¼ ‘¼1 SWAP‘;n‘ : (47)
X
n
The final factorized expression is
x :¼ j=N ¼ xj 2j ¼ 0:x1 x2 x3 …xn : (35)
j¼1
X
n
The quantum Fourier transform (QFT), which is nothing but the QFT ¼ SWAPbn=2c nk¼1 Rðbxc
^ k hÞH j0ihjj: (48)
usual Fourier transform, factorizes into the following form: j¼1
X The main QFT stage has n single-qubit Hadamard gates and
QFT ¼ N 1=2 e2pijk=N jkihjj; (36)
nðn 1Þ=2 controlled two-qubit rotations, for a total of nðn þ 1Þ=2
jk
operations. The final SWAP stage consists of bn=2c SWAP operations.
Xh ‘
i
¼N 1=2
n‘¼1 e2pijk‘ 2 jk‘ i hjj; (37) The swap stage can often be eliminated with a logical reordering of
jk qubits in the next subroutine or by reordering the measurements that
2 3 read out the ancillary register. As shown in Fig. 5, this circuit has a
X X
1
‘ simple graphical structure.
¼ 4‘¼1 2
n 1=2
e 2pijk ‘ 2
jk‘ i5h jj; (38)
j k‘ ¼0
2. Phase estimation
Xh
‘
i
¼ n‘¼1 21=2 j0i þ e2pij2 j1i h jj: (39) Phase estimation allows one to determine the eigenvalue of a uni-
j
tary operator, given the ability to prepare an eigenvector. If the appli-
Now the state repeated within the tensor product is simply cation of the unitary, e.g., through Hamiltonian simulation, is
Rð2pj2‘ ÞHj0i ¼ Rð2n‘ xhÞHj0i, where h ¼ 2p=N ¼ 2p=2n and inexpensive, then phase estimation is also inexpensive. Assume that
the single qubit phase rotation operator is defined as jwa i is an eigenvector of the unitary operator U, so that Ujwa i
¼ eia jwa i. Then repeated applications of the unitary yields informa-
k k
Rð/Þ ¼
1 0
: (40) tion about the less significant digits of a, i.e., U2 ¼ e2 ia jwa i. This
0 ei/ information can be stored in an ancillary quantum memory register
by performing the operations
Thus, one can also write
X X
m X
m
QFT ¼ nk¼1 Rð2nk xhÞH j0ihjj (41) jwi ¼
j1
jji CU2 jwa i ¼
j1
ei2 a
jji jwa i: (49)
j j¼1 j¼1
X
¼ Rð2n1 xhÞH Rð2n2 xhÞH … ½Hj0ihjj: (42) Now performing the QFT on the ancillary register yields the binary
j
The effect of the factor 2nk is to truncate the phase to the last n–k digits of the eigenphase
binary digits of x, i.e.,
X
m
X
nk QFTjwi ¼ ak jki jwa i: (50)
bxck :¼ b2nk xc ¼ xkþj 2j1 ¼ 0:xk …xn (43) k¼1
j¼0 The circuit diagram for quantum phase estimation (QPE) is shown in
so that Fig. 6.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-13

C Author(s) 2023
V
FIG. 5. Circuit diagram for the quantum

Fourier transform QFT. Here, Rk
:¼ Rð2nk hÞ ¼ Rð2p=2k Þ, where the
phase gate Rð/Þ is defined in Eq. (40).
Thus, if the eigenvector can be prepared efficiently and the uni- where f ðxÞ is a binary function that returns the values 0 or 1. In partic-
tary can be simulated efficiently, then the QFT can be used to extract ular, the oracle that marks a single particular state jyi will be denoted
the eigenphase efficiently. Note that there is nothing intrinsically Ojyi if f ðxÞ ¼ dx;y . (Note the minus sign.)
quantum about this algorithm other than these assumptions. In fact, The idea behind amplitude amplification is to generate a rotation
Kitaev et al.5 presented a version of this algorithm that simply uses that amplifies the amplitude of the subspace of interest, called the good
direct measurements of the phase cos ð2k aÞ and sin ð2k aÞ dispenses subspace. Now, the product of two reflections is a rotation around the
with the ancillary qubit register altogether. point of intersection of the two lines of reflection. The rotation angle is
Phase estimation for certain classes of generalized eigenvalue twice the size of the angle between the two lines of reflection. Clearly,
problems, of the form commonly found in St€ urm–Liouville problems, the oracle is a reflection operator, so, if one can find another suitable
magnetohydrodynamics (MHD), and kinetic theory was developed in reflection, then it might be possible to generate a useful rotation. In
Ref. 134. fact, there are an infinite number of suitable choices.
Let S denote a unitary transformation that creates the state jsi
pffiffiffiffi ¼ Sj0i from the vacuum state j0i. This state should be a superposition
3. Amplitude amplification: Cost of N rather than of all states to be searched over; for example, the choice S ¼ H, which
classical N generates the uniform superposition state jsi ¼ Hj0i, is often a useful
Understanding the key components of Grover’s search algo- choice. The amplitudes and phases of the superposition may be chosen
rithm led to the discovery of amplitude amplification and ampli- almost arbitrarily, but only amplitudes that are not exponentially small
tude estimation. While the utility of Grover’s search is often will have a reasonable probability of being searched. The search oracle
debated because the cost of data entry (for an exponentially large operator based on this state is
phone book) must be amortized over an exponentially large
Ojsi ¼ 2jsihsj 1 ¼ SOj0i S† ; (52)
number of searches, the utility of amplitude amplification and
estimation is not in question. Amplitude amplification can be where Oj0i ¼ 2j0ih0j 1 marks the vacuum state j0i. The generalized
used to amplify the overlap between two wavefunctions quadrati- Grover walk operator, G, alternates between the search oracle and the
cally faster than this can be performed classically. This is the key function oracle
part of Grover’s search that allows it to obtain the correct answer
with high probability with quadratically fewer queries that the G ¼ Ojsi Of : (53)
classical counterpart. Now assume one is given an initial state jw0 i. Clearly, this can be
Assume that one has a function oracle operator that marks chosen decomposed into the sum of two parts: a part in the marked subspace,
states by flipping their sign, i.e., denoted jgoodi, and a part in the unmarked subspace, denoted jbadi, i.e.,
Of jxi ¼ ð1Þf ðxÞ jxi; (51) jw0 i ¼ cos ðh=2Þjbadi þ sin ðh=2Þjgoodi: (54)
FIG. 6. Circuit diagram for quantum phase

estimation (QPE), using QFT to readout
the phase to an ancillary memory register.
If the number of ancillary states is a power
of 2, as in the case of an ancillary qubit
register, then the initial QFT can be
replaced by H.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-14

C Author(s) 2023
V
The action of the generalized Grover walk operator is to rotate the eigenvector of G would allow one to directly readout the value of the
amplitude of the good and bad states via amplitude, cos ðhÞ. In this case, we cannot construct an exact eigenvec-
tor, but we can construct the initial state jw0 i ¼ jsi, which is a sum of
jwk i ¼ Gk jw0 i (55) two eigenvectors, as seen from Eqs. (54) and (55). The result of apply-

1 1 ing phase estimation to jsi, as shown in Fig. 8, yields an equal superpo-
¼ cos k þ h jbadi þ sin k þ h jgoodi: (56) sition of the states jhi and jp hi in the ancillary register.104 Clearly,
2 2
both states yield the same value for the amplitude.
The amplitude amplification process is illustrated in Fig. 7.
For example, for Grover’s search, choose the walk operator to be 5. Quantum walks: Quadratic speedup
defined by the uniform superposition state generated by S ¼ H acting
on the vacuum. Assume that one is guaranteed that the function oracle Quantum walks on graphs are a vast generalization of the ampli-
Of marks M out of N states. If one begins with an initial uniform tude amplification paradigm. In its discrete-time version, quantum
superposition, then sin ðh=2Þ ¼ M=N. By choosing the appropriate walks allow state changes based upon the outcome of flipping a quan-
value of k, i.e., tum coin113 (or by rolling quantum dice). Quantum walks accelerate
the solution of various problems, e.g., the search hitting time, typically
k ¼ round ðp=2hÞ 1=2; (57) quadratically. A heuristic explanation is that the ability to explore
pffiffiffiffiffiffiffiffiffiffiffi
lim k ¼ ðp=4Þ N=M ; (58) superposition states increases the effective speed of the walker.
M=N!0 If one includes Grover’s search, amplitude amplification, and
estimation within this group, then the majority of quantum algorithms
one will generate a state that has the probability of measuring good
that achieve a quadratic speedup today use quantum walks in an
states amplified to a high level. Then, simply performing measure-
essential manner.131 Quantum walks have been used to speed up
ments will yield these marked states with high probability. The final
search on graphs,135,136 accelerate simulation of Markov chains,137 and
limiting form occurs in the limit that the dimension of the good sub-
Chebeyshev cooling schedules for quantum annealing,138–140 as well as
space is much smaller than the overall dimension M=N ! 0.
to solve the element distinctness problem.141 Notably, they are useful
for Hamiltonian simulation116,142 as well as for quantum signal
4. Amplitude estimation: Cost of 1/accuracy instead
processing.22
of classical 1/accuracy2
Although quantum walks are deterministic, they share similari-
Once amplified, the less significant digits of the amplitude can be ties with classical random walks which are used as the basis for classi-
measured efficiently. Hence, the amplitude can be measured with qua- cal probabilistic Monte Carlo algorithms that simulate the evolution of
dratically fewer steps than it can by directly sampling the probabil- a classical PDF in time. Due to the Feynman–Kac formula, this evolu-
ity.104 In fact, whenever the answer that one seeks in encoded in a tion can often be ascribed to Langevin dynamics defined by a stochas-
probability amplitude, then the only way to extract this answer more tic forcing of a deterministic system, i.e., a stochastic differential
efficiently than by direct sampling is by using amplitude estimation. equation (SDE). This temporal evolution of the PDF can be used as
Algorithms for quantum counting,104 which can be used to compute the basis of search algorithms that seek to find marked states within a
sums and integrals up to quadratically faster than classical probabilistic certain hitting time. Another important goal is to evolve the PDF
algorithms37,64 fall into this class. Clearly, these algorithms are essential toward a desired statistical distribution which may physically represent
for many scientific computing applications. a canonical ensemble, such as the Boltzmann distribution, the ground
Let the unitary S generate the state of interest jsi ¼ Sj0i whose state of a wavefunction, or the solution to an optimization problem.
overlap with the good subspace, marked by Of , is to be measured. After a few mixing times, the desired PDF is achieved with high accu-
Assuming that both S and S† can be performed efficiently, the search racy and one can generate samples from the PDF and/or calculate the
oracle Ojsi and the Grover walk operator, G, can be constructed partition function. As the temperature is reduced, the PDF corre-
efficiently. sponds more and more closely to the ground state, which leads to the
The essential idea of amplitude estimation is to determine the simulated annealing approach to optimization and ground state
least significant bits of h through the action of powers of the Grover preparation.
walk operator Gk . For example, using phase estimation on an
FIG. 8. Circuit diagram for quantum amplitude estimation (QAE), using QPE to
readout the amplitude. The Grover walk operator operator, G in Eq. (53), is defined
by the choice of oracles and the choice of initial superposition operator S. In this
diagram, the multiply controlled Gx block is defined as the block of controlled-U
FIG. 7. Circuit diagram for quantum amplitude amplification. To amplify to order operations sandwiched between QFT’s within the QPE diagram, Fig. 6, where G is
unity, the amplification block is to be repeated a number of times, k, of order 1/ used in place of U. If the number of ancillary states is a power of 2, as in the case
jamplitudej. of an ancillary qubit memory register, then the initial QFT can be replaced by H.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-15

C Author(s) 2023
V
Szegedy introduced137 a quantum walk analog for any Markov 7. Linear equation solvers, linear differential equation
process defined by a stochastic transition matrix. Like Grover’s walk, it solvers, and matrix operations
is composed of two reflections, but now on the tensor product of two
Harrow, Hassidim, and Lloyd (HHL)123 invented the first general
Hilbert spaces. He proved that the hitting and mixing times of
quantum algorithm for solving sparse linear equations, Ax ¼ b. The
the quantum walk are generically improved quadratically because the
HHL algorithm estimates the eigenvalues of A and replaces them by
spectral gap for the quantum walk is the square root of that for the
their inverse using an ancillary qubit as well as a technique called post-
random walk. This quadratically improves the complexity of using
selection, i.e., measurement of the state of the ancillary qubit. The
simulated annealing to find the solution to combinatorial optimization
method was subsequently improved using the methods of variable
problems in addition to the quadratic speedup in accuracy obtained
time amplitude amplification124 and linear combinations of uni-
through amplitude estimation.138–140 The improvement in conver-
taries.125 General matrix operations can be performed using the QSP
gence also improves the Chebeyshev cooling schedule.37 The complex-
and qubitization techniques.122
ity of general multi-level Monte Carlo methods for SDEs can also After the introduction of linear solvers and advanced
generically be improved quadratically.39 Quantum algorithms based Hamiltonian simulation techniques, methods for solving linear differ-
on these concepts have recently been explored for simulating turbulent ential equations were developed.18,128 For hyperbolic partial differen-
mixing and for computing the reaction rate within turbulent tial equations, such as general wave equations, Ref. 19 found similar
flows.40,41 results. More recent linear differential equation solvers have much
improved dependence on precision by using a linear combination of
6. Hamiltonian simulation: Polynomial to operations18 and by using spectral methods.20,148
superpolynomial speedup The quantum linear solver algorithm (QLSA) and quantum lin-
Simulation of many-body quantum systems was the inspiration ear differential equation algorithm (QLDSA) have an exponential
for the development of quantum algorithms for quantum computers. improvement in the complexity of solving for a quantum encoded ver-
sion of the solution. Note, however, if the amplitudes of the entire
The first efficient algorithms proposed were based on the realization
solution are required, then the speedup is reduced to quadratic. For
that many Hamiltonians of interest have a sparse decomposition into
example, Ref. 149 found that, after properly accounting for precision,
a sum of non-commuting parts, where each part is composed of com-
the improvement for solving elliptic PDEs using finite elements is only
muting terms, and, hence, each part can be computed individually as a
quadratic.
product of terms.108 Then, one can use operator splitting techniques,
such as the Trotter–Suzuki decomposition and Lie group decomposi-
tion formulas,143 to approximate the unitary corresponding to the full 8. Nonlinear solvers
Hamiltonian as a product of terms with sufficiently small time steps. a. Nonlinear operations. The QM Postulates require linear evolu-
In fact, while these techniques are well known to the operator splitting tion of the state, and, if the evolution is ideal, the linear operation must
community, a rigorous formulation of the approximation theory was be unitary. Hence, nonlinear operations of the entire state cannot be
only recently developed.112 performed with a speedup over classical algorithms. As discussed pre-
For example, for the Schr€odinger equation for particles interact- viously, nonlinear quantum eigenvalue transformations (QEVT) of a
ing with few-body potentials, the Hamiltonian breaks into a sum of normal matrix129 and nonlinear quantum singular value transforma-
the kinetic energy of all particles and the potential energy of all par- tions (QSVT)126 can be performed efficiently for a block-encoded
ticles. One can then efficiently switch back and forth between the matrix. This approach has begun to yield efficient methods for per-
momentum and position basis in order to evaluate each as a simple forming nonlinear operations on blocks of the wavefunction.150,151 An
phase shift.110,144 There is also a similar decomposition for lattice alternate approach, based on considering a nonlinear operation on a
Hamiltonians with nearest neighbor coupling. These techniques can given set to be a linear operation on the Hilbert space of functions on
also be applied to both bosonic and fermionic systems as well as to sys- the set was explored by Ref. 152. This method is the analog of the
tems of distinguishable particles.109 Lattice simulations may be directly Koopman–von Neumann approach, discussed below, for finite time
applied to model relativistic quantum plasmas.145,146 intervals.
Quantum algorithms for general sparse Hamiltonians were then There is good potential for these methods to develop into the
developed by a number of authors.115,147 Today, optimal strategies are ability to efficiently perform arbitrary nonlinear operations. If this is
based on Childs’ generalization of Szegedy’s approach to define a the case, then the ability to perform a nonlinear operation, combined
quantum walk for an arbitrary Hamiltonian.116,142 Using strategies with the ability solve linear equations, allows one to develop methods
such as the linear combination of unitaries118 and Taylor series for solving nonlinear equations using either fixed-point Picard itera-
approximations,120 the algorithms were improved to have nearly opti- tion or Newton’s method.
mal dependence on all key parameters, including condition number,
precision, and failure probability, in Refs. 117, 119, and 120 and, even- b. Nonlinear differential equation solvers. The first concrete quan-
tually, to optimal dependence in Refs. 21 and 121. Optimal depen- tum algorithm to simulate nonlinear differential equations was pro-
dence on all parameters was actually first achieved for a large class of posed in Ref. 28, but had an exponential complexity in time. This is
Hamiltonians using the methods of qubitization and quantum signal because the entire wavefunction needed to be stored for each time step
processing (QSP).22,129 This method works by using QSP to correct in order to compute the required nonlinear functions. After this result,
the eigenvalue distribution of the Childs walk so that it generates a the quantum algorithms community began to focus on solvers for lin-
more accurate approximation of the underlying Hamiltonian. ear differential equations, as described in Sec. III E 7.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-16

C Author(s) 2023
V
The next conceptual breakthrough was provided by Ref. 29, associated with the circuits that are to be optimized to minimize the
which proposed an algorithmic approach to solve nonlinear differen- cost of the objective function. Quantum resources are useful because
tial equations by mapping the nonlinear problem to the linear they allow the use of Hamiltonians with many degrees of freedom and
Koopman–von Neumann representation of classical dynamics. After of ansatzes that are classically expensive. Thus, they have a few key
this step, Hamiltonian simulation could be directly applied to the sys- independent parts: (i) constructing the cost function via Hamiltonian
tem. Then, Liu et al.31 derived an explicit algorithm for differential simulation, (ii) constructing a suitable ansatz for the initial wavefunc-
equations with quadratic nonlinearities based on the closely related tion, (iii) developing efficient measurement protocols to determine the
technique of Carleman linearization. The connections between these cost function and, potentially, its derivatives with respect to optimiza-
techniques will be elucidated in detail in Sec. V. tion parameters, and (iv) the classical optimization algorithm.
Another algorithmic approach for solving the dissipative nonlin- While it is often hoped that these algorithms will be useful for
ear differential equations that arise from many-body physics was also achieving at least quadratic quantum advantage on near-term NISQ
recently advocated by Ref. 153. The basic idea, which goes back to ini- hardware platforms, there are a number of outstanding challenges that
tial considerations of quantum computing,154,155 is that such a descrip- must be considered and these challenges can make it difficult to obtain
tion is an approximation of many-body quantum physics. Hence, a rigorous proofs of their complexity. Notable issues are associated with
simulation of the many-body quantum physics Hamiltonian should developing efficient ansatzes for the initial wavefunction, developing
lead to an approximate solution of the classical system in the appropri- efficient measurement techniques, avoiding barren plateaus with little
ate limit. information on the optimization problem, and the presence of many
Yet another approach for solving the Navier–Stokes equations of local minima, which make the problem intrinsically difficult. In fact,
fluid dynamics was explored by Gaitan156 based on a method for inte- the classical-quantum hybrid scheme may not have any advantage
grating differential equations developed by Kacewicz.38 In turn, the when the classical optimization problem is NP-hard.162
time-integration method developed by Kacewicz relies on quantum IV. QUANTUM SIMULATION
algorithms developed by Heinrich and Novak for computing sums
A. General framework
and integrals.37,64 As explained previously, these methods can obtain
up to a quadratic quantum speedup over Monte Carlo methods for The quantum simulation application space is closest to the typical
high-dimensional systems with non-smooth right hand sides. applications used within FES. Several authors17,29,31,32,163 have outlined
However, a key point is that this speedup occurs when only partial a framework for the quantum simulation of nonlinear dynamics. As
information is given to specify the right hand side,157 i.e., only a certain discussed in Sec. III, there are a number of subroutines to be employed
number of derivatives exist and the position of possible discontinuities that have steadily improved in computational complexity and accu-
in the last derivative are not known in advance. Thus, such gains can racy. A summary of the key steps is discussed in this section.
only achieved if there are non-smooth initial conditions whose evolu- The key conceptual steps are as follows:
tion is too difficult to track in time or if there is a stochastic forcing of (1) Determine an embedding of the dynamics into a linear differen-
the differential equations. Stochastic forcing is in fact useful for studies tial system of equations.
of fluid and plasma turbulence as well as for stochastic differential (2) Numerically discretize the linear system.
equations in general. However, it is important to keep in mind the fact
that this method would not accelerate the solution of ODEs and PDEs If the system is intrinsically linear, like a quantum Hamiltonian,
with analytic coefficients and initial conditions. then the first step is relatively straightforward, but the choice of repre-
sentation can still have an impact on complexity, which is especially
9. Variational quantum algorithms important when resources are limited. If the system is nonlinear, then
a linearization must be determined. In either case, it is important to
Variational quantum algorithms133,158–160 are typically formu- choose the basis so that the resulting equations have a sparse structure
lated as hybrid quantum-classical algorithms for solving optimization or an efficient Trotterization that can be utilized by the quantum simu-
problems. All of the quantum subroutines discussed above require lation algorithms.
large-scale fault-tolerant quantum computers which will not be avail- The key computational steps are as follows:
able in the near term. As a middle ground for the NISQ-era, hybrid
(1) Prepare the initial state.
quantum-classical algorithms have been proposed where quantum
(2) Evolve the system in time using an efficient linear differential
devices are envisioned as co-processors that accelerate the otherwise
equation solver, or if the system is Hamiltonian, an efficient
classical calculations. It is hoped that these will serve as NISQ algo-
Hamiltonian simulation algorithm.
rithms for difficult many-body problems in quantum chemistry and
(3) Estimate the observable to be measured efficiently, typically by
quantum materials science, typically involving fermions. They may
using an ancillary quantum memory register to encode the
even be useful as methods for solving nonlinear problems.44 For exam-
observable and performing amplitude estimation.
ple, variational simulation of the Navier–Stokes equations has recently
been nvestigated.45–47 Variational simulation of nonlinear stochastic If the linearized system is sparse, then using the sparse linear dif-
differential equations has also been explored.161 ferential equation solver ensures that the key steps are efficient.31 If the
The idea is that there is a cost function of high complexity that is linearized system is both sparse and Hamiltonian, then using the
evaluated by a quantum computer which depends on a number of sparse Hamiltonian simulation algorithm ensures that the simulation
parameters that is minimized by using a classical optimization algo- step is efficient.29 State preparation is required for the initialization
rithm. For example, the parameters may be in the form of phase angles step and is often required for the final measurement stage. For an
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-17

C Author(s) 2023
V
efficient algorithm, one must choose initial conditions that can be pre- finite collections of orthogonal polynomials, and complex analytic
pared using an efficient method such as Hamiltonian simulation or a function spaces, such as Segal–Bargmann coherent states, Bergman
quantum walk.29 spaces, and Hardy spaces. For example, as proven in Appendix C, the
Many works on quantum computing algorithms consider the Carleman representation is based on the standard Hardy space, a com-
output of a subroutine to be a quantum encoding of the solution, i.e., plex analytic RKHS that uses a Sz€ego kernel. This RKHS is based on
wavefunction or density matrix that can then be used as input for the collection of monomials, z k , where evaluation is given by the Sz€ego
another subroutine. While this is quite reasonable for intermediate kernel for the boundary of the unit disk in the complex plane.
stages of a quantum program, for almost all scientific calculations of Determination of which representations are more or less efficient for
interest, the output usually consists a set of real floating point numbers representing different applications is an active research area.
that quantify the system of interest. Whenever these data are encoded The Stone–von Neumann theorem implies that all representa-
in the amplitude of the wavefunction, then one requires a final stage tions that satisfy the (Weyl exponential form of the) canonical commu-
where amplitude amplification and estimation is used to extract the tation relations (CCR) are unitarily equivalent.166 While there are
quantities of interest. important exceptions, the vast majority of the representations used in
practice are unitarily equivalent to the position and momentum opera-
B. Sparse representation and numerical discretization tors; i.e., there is a unitary transformation that relates the two different
representations. For example, the number basis defined by normalized
Choosing a finite representation of the wavefunction is clearly
Hermite polynomials, Hk ðxÞ are related to the position representation
required for numerical simulation. Choosing a basis that keeps the
operations sparse is clearly an important consideration. In general, a by a unitary transformation. Similarly, the coherent states are related
complete basis allows one to represent wðxÞ, as a to the number states via a unitary transformation. The combination of
the wavefunction,
these two transformations generates another unitary transformation
linear combination of basic elements, /k ðxÞ so that
from position space to coherent state space known as the
X
N Segal–Bargmann transform.
wðx; tÞ ¼ wk ðtÞ/k ðxÞ: (59) There have already been some interesting uses of RKHS’ within
k¼1 quantum computing for nonlinear dynamical systems. Reference 163
The variables x are indices representing the computational basis. If the used a real RKHS that used smoothed basis functions to develop a
computational basis represents position, then it might be convenient well-defined Hamiltonian simulation algorithm for integrable systems
to use the position basis /k ðxÞ ¼ dkx or the momentum basis, that displays smooth pointwise convergence in the limit N ! 1.
/k ðxÞ ¼ exp ð2pikxÞ=N 1=2 . More generally, the discretized Hilbert Reference 31 used the Carleman embedding, which is based on the
space could span any finite collection of linearly independent func- standard Hardy space, in developing a quantum algorithm for dissipa-
tions, such as orthogonal polynomials, Bessel functions, etc. One could tive systems, such as the logistic equation and the Burgers equation.
also use finite difference or finite element basis functions, i.e., orthogo- Reference 32 considered the possibility of exploring other representa-
nal polynomials (or, more generally, orthogonal functions) over sub- tions of the CCR. As proven in Appendix C, all of the concrete exam-
domains with finite support. If the differential equation to be solved ples discussed in Ref. 32 are well-known examples of complex analytic
has an efficiently computable sparse representation in terms of these RKHS’.
basis functions, i.e., s terms, then the resulting system of equations will
be s-sparse. C. State preparation: Initialization and measurement
One can also use an overcomplete basis, such as the coherent In the typical models of quantum computing, the initial wave-
states. Coherent states are an important example of a reproducing function must start in a standard state, such as the ground state j0i.
kernel Hilbert space (RKHS), which is defined by the fact that the This is useful from the point of view of a physical device, because, if
pointwise evaluation of functions is continuous.164,165 In this case, the the temperature is much lower than the energy difference between the
Riesz representation theorem implies that there is a continuous linear first excited state and the ground state, then there is an exponentially
kernel function, Ky ðzÞ, that allows one to evaluate a function at the small likelihood of finding the initial wavefunction in an excited state.
point y via Due to dissipative decoherence processes such as relaxation, simply
ð waiting long enough allows one to re-initialize the quantum program
wðyÞ ¼ hKy jwi ¼ Ky ðzÞwðzÞdz: (60) from j0i. Active restart techniques intentionally use decoherence to
speed up the process of resetting to vacuum.
This is called the reproducing property of the kernel. The definition of Then, to actively begin the program, one must typically initialize
the bilinear kernel function in terms of the Hilbert space inner prod- the wavefunction in a specific state. For simple programs, such as
uct, Kðy; zÞ ¼ hKy jKz i, demonstrates that the kernel function must be search, one might use the Walsh–Hadamard transform to prepare a
symmetric and positive definite. The Moore–Aronszajn theorem states uniform superposition of a quantum memory register. For a simula-
that every symmetric positive definite kernel defines an RKHS. tion, one would typically initialize the memory register with a state
There are actually many types of RKHS’ that are used within the corresponding to the initial condition of interest. As discussed previ-
machine learning community for efficient representations of data in ously, this initial condition should be efficiently computable, perhaps a
higher-dimensional spaces. One may expect that RKHS’ should play a smooth function with some random noise which might be useful for
similar role in developing efficient representations for quantum com- triggering the dynamics of interest, e.g., as in simulations of turbu-
puting and should certainly be of interest for quantum machine learn- lence. As discussed previously, many useful functions can be com-
ing (QML). Examples of RKHS’ are bandwidth-limited functions, puted efficiently by using quantum walks or Hamiltonian simulation.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-18

C Author(s) 2023
V
In order to obtain output from a given memory register, a projec- time steps, of the order of dt h=jjHjj. Hence, the total number of
tion step is usually required to obtain the data of interest. In practical time steps, Nt ¼ Dt=dt jjHjjDt=h, will typically be large. For exam-
models of quantum computing, there is typically only a possibility of ple, to observe temporal variations on the longest relevant timescale
performing measurements in the computational basis. Hence, in order would require Dt h=minðeigenvaluesðHÞÞ, so the number of steps
to extract other properties of the wavefunction, e.g., measurements will typically scale as the condition number Nt jðHÞ. Hence, while
along another axis of the Bloch sphere, one must perform judiciously the QLDSA is more costly than Hamiltonian simulation, particularly
chosen unitary transformations before the measurement is taken. This with regard to condition number, the QLDSA is similar to that of
final projection step can either be thought of as a transformation of Hamiltonian simulation in terms of complexity in time steps. Thus,
the data by some unitary operation U or as a transformation of the which choice is optimal may depend on other desired features of the
state to be measured against with the inverse operation U† . solution.
D. Simulation E. Estimation of physical observables

After defining an efficient representation, the simulation step has In fact, there are some relatively simple methods for extract-
now reduced the problem to solving for the evolution of a linear set of ing observables from a set of measurements.29 To simplify the dis-
ODEs. In components, these ODEs can be written as cussion, here we consider the case where all observables commute
X with each other (are in involution). The expectation value of a local
ihdwj =dt ¼ Hjk wk ¼ h/j jHjwi: (61) observable
k X
hOi ¼ OðxÞpðxÞ (64)
If the evolution is unitary, then one can directly use the best avail- x
able sparse Hamiltonian simulation algorithm (HSA).108,109,117 This
constructs an approximation of the form can be found by performing a reversible classical computation of the
Ð functions
wðtÞ ¼ Uapprox w0 ’ tei Hdt=h w0 ; (62)
/ðxÞ ¼ O1=2 ðxÞwðxÞ; (65)
where t is the time ordering operator. The result will typically be qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
approximate due to the need to simulate the Hamiltonian efficiently. If /0 ðxÞ ¼ 1 j/ðxÞj2 : (66)
the Hamiltonian can be broken into a small number of pieces where First, form the states
each piece is easy to construct, e.g., because it is sum of terms that
commute with each other, then one can use a Trotter–Suzuki decom- X rX
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
position to generate an accurate approximation. Otherwise, one can j/i ¼ /ðxÞjxi j/ðxÞj2 : (67)
x x
choose one of the sparse Hamiltonian solvers described previously.
If the evolution is not unitary, then one can use the best available Then, to compute the amplitude, add an ancillary qubit to the system,
quantum linear differential equation solver algorithm (QLDSA).18,128 j/i j0i, and perform a specific choice of rotation, R/ , defined via
Simple statements of these algorithms might often an explicit forward
X
N pffiffiffiffi
Euler time step; however, this can clearly be improved to use higher R/ j/i j0i ¼ jxi /0 ðxÞj0i þ /ðxÞj1i = N (68)
order methods as well as implicit timestepping, using the QLSA. For x¼1
example, for each time step, one can consider performing a time split- ¼ cos ðhÞj/0 i j0i þ sin ðhÞj/i j1i: (69)
ting of the form
Finally, using the amplitude estimation subroutine, discussed in Sec.
½1 þ ihHDt=hwðt þ DtÞ ’ ½1 ið1 hÞHDt=
hwðtÞ: (63) III E 4, allows one to directly calculate the amplitude of the j1i state
Here, h represents a time splitting parameter that determines the to determine sin2 ðhÞ ¼ hOi=N with an amount of work that scales as
amount implicitness to be used; i.e., h ¼ 0 for fully explicit forward 1/accuracy rather than 1/accuracy2.
Euler, h ¼ 1=2 for implicit midpoint Crank–Nicolson, and h ¼ 1 for
fully implicit backward Euler. F. Quantum speedup
In order to handle a non-unitary evolution, the QLDSA must use For both the intrinsically quantum and the intrinsically classi-
post-selection. Care must be taken to be able to extract the information cal approaches, quantum Hamiltonian simulation of classical
of interest and this depends on the properties of the solution itself. If dynamics leads to large gains in memory and computational com-
the linear system is unstable and the magnitude of w is exponentially plexity. The Hamiltonian simulation algorithms are efficient as
unbounded, then the output will likely result in noise. If the linear sys- long as the Hamiltonian is sparse. It is a simple exercise to prove
tem is stable and the magnitude of w shrinks exponentially toward that the Hamiltonians of interest are sparse for N-body interac-
zero, then the output will also eventually result in noise. tions and for particles interacting with electromagnetic fields. This
For realistic applications of interest, the no-fast-forwarding theo- implies up to exponential savings relative to an Eulerian discretiza-
rem holds.117,147 The no-fast forwarding theorem implies that in the tion of the Liouville equation.
generic situation where the solution cannot be solved for semi- However, the best classical algorithms of this sort are time-
analytically, the solution must then be tracked on the shortest relevant dependent Monte Carlo methods.167 For example, particle-in-cell
time steps. This implies that the simulation will be restricted to short algorithms invented in the plasma physics community fall into this
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-19

C Author(s) 2023
V
class, as well as molecular dynamics algorithms invented in the con- A. Quantized dynamics
densed matter community. These probabilistic classical algorithms can If one begins with a Hamiltonian system, then a natural approach
also provide an exponential speedup over Eulerian methods when the is to use a quantized version of the classical physics model. There are
same restriction on the observables of interest are applied. In this case, many potential advantages to this approach. First of all, the quantized
one does not need to accurately sample the entire PDF, rather, one version may be the more accurate physical model, so important effects
only needs to maintain good accuracy for the most probable regions of such as interference, tunneling, and diffraction would be correctly
phase space.167 In this case, the equivalent quantum algorithm can described. Second of all, there are already many quantum algorithms
achieve up to a quadratic speedup over time-dependent Monte Carlo for the simulation of sparse Hamiltonians, and, this is especially true
(e.g., PIC) methods.168,169 for quantum Hamiltonians.
It is important to stress that it is the use of amplitude estimation Two examples of this approach, as well as their implementation
of physical observables and the use of quantum walks for state prepa- on present-day quantum hardware, will be discussed in detail in
ration and measurement transformations that enable the quadratic Sec. VI. In particular, Sec. VI B explores the physics of the sawtooth
speedup. Due to the central limit theorem, when classically sampling map,170 which can be viewed as a poor-man’s model for wave–particle
an observable with a finite mean and variance, the relative accuracy interactions that uses minimal numerical resources. Figure 9 compares
decreases as e 1=S1=2 where S is the number of samples. This the dynamics of the classical sawtooth map to that of the quantum
implies that the amount of computational work required to achieve a sawtooth map as the number of qubits is increased from 6 to 9 to 16.
certain accuracy grows as 1=e2 . The utility of the amplitude estimation To visualize the quantum dynamics in phase space, the Husimi-Q
algorithm is that it only requires an amount of computational work quasi-probability distribution is used to visualize the density matrix. In
that grows as 1=e. both cases, the initial state has uniform amplitude along the p=p
¼ 3=4 momentum surface, and, as time evolves, both the quantum
V. QUANTUM REPRESENTATION OF NONLINEAR wavefunction and the classical PDF explore similar regions of phase
DYNAMICS space. However, this comparison only works well once the number of
The majority of simulations used within FES today aim to solve quantum states becomes large enough. At n ¼ 9 and N ¼ 29 ¼ 512
initial-value and boundary-value problems: typically time-dependent states, for a total phase space resolution of 512 512 ’ 2:6 104 , the
simulations of nonlinear dynamics. While the system may often repre- quantum Husimi-Q distribution still appears quite “fuzzy” to the
sent partial differential equations (PDEs), such as the Navier–Stokes human eye. At n ¼ 16 and N ¼ 216 ’ 6:6 104 , for a total phase
or magnetohydrodynamics (MHD) equations of fluid mechanics, or space resolution of 232 ’ 4:3 109 , quantum artifacts are no longer
visible at this scale.
the kinetic equation for the probability distribution function, ulti-
This example makes it clear that, for the approach of quantiza-
mately, numerical discretization of the set of PDEs typically reduces
tion, it is not likely that one will have enough quantum computational
the system to time integration of a set of ODEs. The question then
resources to solve the exact quantum system underlying the classical
arises: How should quantum computers be used to efficiently simulate
problem. Instead, the basic quanta of action used to compute particle
nonlinear dynamics?
numbers, i.e., the effective Planck’s constant would be larger, effec-
The past few years have seen intensive research activity on
tively grouping large numbers of particles together in a manner that is
the ability to perform quantum simulations of nonlinear dynam-
similar to the way that particle-in-cell codes evolve macroparticles that
ics. As described above, these approaches generally rely on each represents large groups of physical particles. If the model is for-
embedding the nonlinear system into a large, in fact, infinite- mulated correctly, parameters of the problem could follow an appro-
dimensional linear system of equations.17,29,32,33 In this section, priate renormalization group flow, such that simulating a smaller
we simultaneously unify these approaches and systematically cate- number of groups still produces the correct results for a much larger
gorize their differences. system. Interestingly enough, this will generate numerical discretiza-
Classical dynamics is distinguished by the fact that the solution tion error of a type not commonly explored in the traditional applied
to a well-posed set of nonlinear differential equations typically has mathematics community.
both existence and uniqueness theorems that apply to almost all points Moreover, depending on the application, there can also be a
in phase space, aside from a set of measure zero. The uniqueness of number of disadvantages. First, taking the semiclassical limit requires
the individual trajectories then allows one to associate a unique evolu- very large quantum numbers for many parts of the system. As men-
tion of other objects, such as scalars, vectors, and tensors on phase tioned above, accuracy of the physical model must be balanced against
space. In particular, the Liouville equation expresses the fact that a vol- the numerical discretization errors, e.g., due to using a value for
ume form, which can be interpreted as a classical probability density, Planck’s constant that is likely to be too large. Second, the quantum
is advected by the flow in a specific manner that ensures that probabil- system is not equal to the classical system. Quantum dynamics shows
ity is conserved in time. The action of the evolution on phase space many well-known non-classical effects, such as interference, tunneling,
functions presents a linear but infinite-dimensional system of equa- and diffraction. If the application requires understanding the behavior
tions. Due to the Stone–von Neumann theorem, the vast majority of of a single classical trajectory, then these quantum effects may be
all possible representations are equivalent to this one, simply defined highly undesirable. If important decisions are to be made on the basis
by translations in phase space. of classical ODE calculations, then unwanted quantum effects could
From an algorithmic point of view, the next question that arises have catastrophic consequences.
is: How well does a finite truncation of the linear representation approx- The method of quantization is also not necessarily directly appli-
imate the underlying nonlinear dynamics? cable to non-Hamiltonian systems, e.g., those that include dissipation.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-20

C Author(s) 2023
V
FIG. 9. A sufficiently large number of

qubits captures the semiclassical limit of
the quantum sawtooth map corresponding
to the Hamiltonian in Eq. (93). Here, the
parameter K ¼ 0:1, so the system is in
the anomalous diffusion regime. (Left) The
wavefunction is initially in the state p=p
¼ 3=4 and evolved for 500 timesteps.
The colored plot shows the Husimi-Q qua-
siprobability function averaged over the
last 50 time steps for different qubit num-
bers n ¼ 6, 9, 16. (Right) Poincare plot of
the classical sawtooth map for 500 time
steps starting from the initial condition
p=p ¼ 3=4 and various angles.
Yet, reduced physical models that include dissipative processes, such precisely for the benefit of classical applications of the type that are of
as friction, diffusion, and collisions, are the de facto standard used interest to FES.
within high-performance FES simulations today.
In principle, dissipative dynamics can be simulated using an B. Semiclassical dynamics
open quantum system that interacts with its environment. Perhaps
1. Remember what to hold fixed
one route could be to use measurements to generate effective dissipa-
tive interactions with the physical environment. However, in order for Consider a set of ODEs for coordinates z ¼ fz i g on a manifold,
such interactions to be controlled, one must embed the dissipative m. In local coordinates, the ODEs are presented in the form
system into a larger ideal system that is simulated by the quantum
computer. In fact, many of the dissipative dynamical systems of dzi =dt ¼ V i ðt; zÞ: (70)
interest to FES are naturally realized as reduced order models of
In order to understand this equation, it is crucial to understand that
large ideal systems. Ideas such as simulating a many-body quan-
the initial conditions z0 :¼ zðt ¼ 0Þ are being held fixed during the
tum fluid to approximate a classical viscous fluid were explored
evolution of the trajectory. Thus, the solution is actually a multi-
early on in the quantum computing community.154,155 In princi-
variable function, zðt; z0 Þ, that depends on both time and on the set of
ple, simulating a modest-sized ideal system yet only measuring a
initial conditions, z0 . This simple fact implies that one is actually solv-
subset of this system—while controlling discretization errors—
ing a set of PDEs of the form
could be an efficient approach.
Due to the Schmidt decomposition, one can generate an arbitrary @t z i jz0 ¼ Vðt; zÞ rz i ; (71)
mixed state of dimension N N from the reduction of a pure state in
a bipartite quantum system of dimension N 2 . Thus, for a given system where here r :¼ f@zi g.
of interest, the size of the Hilbert space that is needed for the embed- It is also possible to derive the evolution of the initial conditions,
ding may need to be only be as small as the dimension of the system z0 ðt; zÞ, as a function of the final coordinates z. Using the chain rule,
squared, N 2 . The benefit is that would only require an embedding one finds that
with twice as many qubits, 2n, rather than a macroscopically large sys-
tem. Recently, there has been renewed interest in these methods153 @t jz0 ¼ @t jz þ @t z j jz0 @zj ¼ @t jz þ V r: (72)
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-21

C Author(s) 2023
V
Thus, the initial conditions satisfy the equation of motion to nonlinear PDEs, which are ubiquitous in FES. One approach is to
first discretize the PDEs to become nonlinear ODEs, and to then use
@t z0i jz ¼ Vðt; zðt; z0 ÞÞ rz0i : (73) the ODE embedding schemes discussed below. However, this
Notice that this equation has the opposite sign of the velocity so it approach yields a linear problem of high dimensionality and it remains
appears to go “backward in time.” However, this is absolutely consis- to quantify the quantum advantage expected for this case. An alterna-
tent with the definition of initial and final coordinates. tive approach has been developed for nonlinear hyperbolic PDEs. In
Consider a classical observable, Oðt; zÞ, that is a scalar function this case, they can be directly embedded in linear PDEs of higher
which is constant along the trajectory, so that Oðt; zðt; z0 ÞÞ dimension using the level-set embedding.36,175,176 It still remains to
¼ Oðt ¼ 0; z0 Þ. Applying the chain rule again implies that proven whether these embedding schemes will offer quantum advan-
tage for FES-relevant problems.
@t Ojz ¼ Vðt; zÞ rO: (74)
If, instead, one would like to express the observable as a function of 2. Prequantization
the initial coordinates via O0 ðt; z0 Þ :¼ Oðt; zðt; z0 ÞÞ, then again using The Liouville equation, which already linearizes the ODEs, can
the chain rule yields be further mapped to the Schr€odinger equation using the approach
developed by Koopman and von Neumann177–180 almost a century
@t O0 jz0 ¼ Vðt; zÞ rO0 : (75)
ago. By definition, the norm of the wavefunction squared represents
In the dynamical systems literature, the evolution equation of scalar the classical PDF, i.e., a volume form. Thus, the wavefunction itself
observables is the infinitesimal generator of the Koopman evolution must be treated as the square-root of a volume form, also known as a
operator. half-form in the language of geometric quantization. The law for the
The adjoint of the scalar equation (over the standard phase space evolution of a form of this type can be obtained by expressing the
measure dz) yields the evolution of a volume-form, pðt; zÞ, on phase wavefunction as
space. The Liouville equation for volume forms is
w ¼ p1=2 exp ðihÞ; (78)
@t pjz ¼ r ðVpÞ: (76)
where h is a scalar function and p is a volume form. If one assumes
In the dynamical systems literature, this is the infinitesimal generator that both h and p are simply advected along the flow, then this yields
of the Perron–Frobenius evolution operator. Again using the chain the standard Koopman–von Neumann (KvN) equation
rule, this equation can be expressed in the completely equivalent form
1
@t wjz ¼ ðV r þ r V Þw: (79)
@t p0 jz0 ¼ r ðVp0 Þ: (77) 2
The representations above converts ODEs that are finite- Crucially, this form is both linear and unitary. Once again, using the
dimensional but nonlinear into PDEs that are linear. The trade-off is chain rule, this can also be expressed as
that the PDE should be viewed as an infinite-dimensional set of ODEs.
1
The representation that focuses on observables is analogous to the @t w0 jz0 ¼ ðV r þ r V Þw0 : (80)
Heisenberg picture in quantum mechanics, and the representation 2
that focuses on the volume form is analogous to the Schr€odinger pic- More generally, let us assume that the complex phase of the
ture. These two pictures are complementary and equivalent. The wavefunction evolves locally along the trajectory, with a source func-
advection equation satisfied by observables in Eq. (74) is linear in O, tion, Lðt; zÞ that is a function of the coordinates of the trajectory. Since
which can be discretized, for example using upwinding or discontinu- the phase is a scalar, the evolution law must be of the form
ous Galerkin schemes, to a set of linear ODEs and solved using quan-
tum algorithms described in Sec. III E 7. Similarly, the continuity dh=dt ¼ @t hjz þ V rh ¼ L=h: (81)
equation satisfied by the volume form in Eq. (76) can be directly solved This assumption then yields a generalized form of the KvN equation
using quantum algorithms.
On classical computers, one may not want to solve nonlinear 1
@t wjz0 ¼ ðV r þ r V Þw iLw=h: (82)
ODEs using the detour of linear embedding, because discretizing the 2
PDEs could introduce numerical artifacts and the resultant linear In this case, time reversal requires using the chain rule and taking the
ODEs are of much higher dimensions than the original problem. complex conjugate which leads to
Nevertheless, for quantum computers, this detour offers a compromise
that enables nonlinear ODEs to be solved in the first place. The cost of 1
@t w0 jz0 ¼ ðV r þ r V Þw0 jt þ iLw0 =h: (83)
increasing the dimensionality and introducing discretization artifacts 2
may eventually be offset by quantum speedups that are not available If the equations of motion are Hamiltonian, with Hamiltonian
on classical computers. For example, Koopman methods are quite HðzÞ and canonical coordinates z ¼ fq; pg, then there is a special
popular method for using machine learning to develop reduced-order choice of the source function that matches Feynman’s prescription for
models of dynamical systems from both simulated and experimental the phase space path integral. This choice is simply the classical
data.171–174 Lagrangian
In this paper, we primarily focus on the embedding of nonlinear
ODEs, but it is worth pointing out that similar methods can be applied LðzÞ ¼ H p @p H: (84)
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-22

C Author(s) 2023
V
Differential operators of the form Once one understands that there is, in fact, a natural geometric
hlH w ¼ H p @p H þ ihfH; wgPB
i (85) setting in which to view the evolution of the wavefunction, it is natural
to try to compare this theory with quantum mechanics itself. Of
181
were first studied by van Hove in his Ph.D. thesis. Van Hove noted course, a typical quantum Hamiltonian, such as the Schr€odinger equa-
that these operators have the special Lie group property tion, is a PDE on configuration space, fqg, that is at least second order
in derivatives. As is clear from the path integral, the quantum evolu-
fihlA ; ihlB gPB ¼ ihlfA;BgPB (86)
tion operator is a weighted sum over all possible trajectories. Thus,
that was emphasized by Dirac in his study of canonical quantization. unlike the classical system, one cannot unambiguously invert the
In fact, this reproduces the classical Poisson algebra of observables on Eulerian coordinates to become a function of initial Lagrangian coor-
a Poisson manifold. Because these operators also represent the first dinates. This singles out the Eulerian picture as the approach that is
order asymptotic expansion of the Schr€odinger equation, Bertram most natural for comparing classical, semi-classical, and quantum
Kostant referred to the construction of these operators as the process dynamics.
of pre-quantization. Pre-quantization is the first step in the program of
geometric quantization which seeks to determine a mathematically 4. Integrable systems
consistent framework for defining quantization on arbitrary symplec-
The case of integrable classical dynamics is discussed in some
tic manifolds.166 Using the terminology introduced by a number of
detail in Ref. 29 and an explicit quantum algorithm for simulation of
authors,182,183 this evolution equation
Hamiltonian systems was derived in Giannakis et al.163 In this case, a
lH w ¼ 0: (87) set of action-angle variables, fJj ; hj g can be determined in which the
action variables are constant in time and the angle variables have a
is now called the Koopman–van Hove (KvH) equation. constant rate of change
Notice that although KvN and KvH equations may resemble the
j
Schr€odinger equations for quantum systems, they are completely J_ j ¼ 0 h_ ¼ xj ðJÞ: (88)
equivalent to the Liouville equation, which describes classical dynami-
cal systems. A crucial difference is that the KvN Hamiltonian only The total phase space is the tensor product of the domain of the action
involves first order derivatives while the non-relativistic Schr€odinger variables and the phase angle variables. An integrable system can
j 1=2 j
equation involves second-order derivatives. Although the Dirac equa- always be recast as a linear system for the variables z6 ¼ Jj e6ih ,
tion, which describes relativistic quantum systems, also only involves since
first order derivatives, each component of the Dirac spinor satisfies a
second-order hyperbolic PDE. The mathematical behavior of first- z_ j ¼ 6xj ðJÞz j : (89)
order and second-order hyperbolic PDEs is rather different and often
Hence, the quantum linear differential equation solver (QLDSA)
requires different types of numerical techniques for their solution.
would certainly work.
It is clear that the evolution of the angle variables can be encoded
3. Forward or backward?
in the evolution of the phase of a single qubit. Hence, it would be sim-
As we have shown in Secs. V B 1 and V B 2, classically, it does not ple to associate each trajectory in phase space with a single qubit.
matter whether the forward or backward form is used because both Then, phase space evolution could be accomplished by single-qubit
equations are completely equivalent. The real question is whether the RZ rotations. However, this would be an exponentially inefficient use
wavefunction is considered to be a function of the final coordinates, z, of resources.
i.e., the Eulerian picture, or of the initial coordinates, z0 , i.e., the In fact, an n-qubit register represents a linear system with N ¼ 2n
Lagrangian picture. From the point of view of the coordinates with states. Thus, it is actually possible to associate each trajectory in phase
fixed initial conditions, z0 , this evolution is forward, but from the point space with a single state in Hilbert space. If the trajectories are labeled
of view of scalar functions, with fixed coordinates, z, this evolution is by a finite set of action variables, Jj, that are held fixed during the evolu-
backward. This answers the questions posed by Ref. 33 on the convention, then each trajectory experiences its own evolution. The KvN
tions used for backward vs forward evolution. Hamiltonian is simply
The use of Koopman evolution operators in the dynamical sys- X
tems community171,172 follows from the desire to use data, collected H¼ xj ðJÞLj : (90)
j
either from experiment or simulation, to deduce the underlying equa-
tions of motion, and, then, given a data-driven model of the dynamics, where the operators Jj ¼ Jj are simply the action coordinates that
to control the dynamics of the system. Thus, this approach is usually index each trajectory and the Lj ¼ ih@=@hj are a set of angular
applied to the evolution of the coordinates themselves, using the momentum operators that are conjugate to the phase angles hj.
Lagrangian framework, rather than to the evolution of scalar functions. In order for an efficient quantum algorithm to exist for this sys-
Since the dynamical systems community often uses Carleman tem, the phase space must first be transformed to action-angle varia-
linearization, recent work on quantum algorithms for nonlinear bles. In general, this is an intractable problem. However, let us proceed
dynamics31,32,34 has used the Lagrangian framework. However, the by assuming that this generically nonlinear coordinate transformation
Carleman representation can be applied to the Eulerian formulation can be accomplished efficiently. Even in this case, it may not be easy to
equally well. In fact, the Eulerian framework might be more natural simulate the Hamiltonian in Eq. (90). The reason is that performing
from the point of view of quantum mechanics. an arbitrary phase shift of each state can be an exponentially difficult
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-23

C Author(s) 2023
V
problem. One of the most efficient algorithms for performing such a KvN Hamiltonian, there are issues associated with the fact that the
phase shift is given by Ref. 184, but has exponential complexity, 2N, in Hilbert space only exactly reproduces phase space in the limit that an
the worst case scenario. Alternatively, the method of qubitization and infinite number of basis functions is retained. For this reason, Ref. 29
QSP allows one to form an efficient polynomial approximation of discussed the potential need to smooth the phase space basis functions
smooth functions,22 so this method may be preferable for many and initial conditions, etc. In fact, a particular choice of smoothing
Hamiltonians of interest. Hence, it is clearly important to continue operation was advocated in Ref. 163 through their choice of RKHS.185
improving these algorithms.
Giannakis et al.163 improve this basic algorithm through the use
of an RKHS approach that acts to smooths the function space of inter- 2. Carleman, coherent states, and complex analytic
est.185 Their algorithm achieves an exponential advantage, in agree- RKHS
ment with Ref. 29, even when considering output. Note, however, that Recently, explicit quantum algorithms for quadratic differential
the cost of the full output of the dependence of an observable over the equations were developed using Carleman linearization,31,186 which
entire phase space would presumably reduce the speedup to quadratic was able to obtain an exponential speedup over the method of Ref. 28
at best. Moreover, they discarded the dependence on the fJj g, which when dissipation dominates nonlinearity. The algorithms are applied
occurs generically for Hamiltonian systems and, as discussed in the to simulate the logistic equation, the Burger’s equation, and reaction–
previous paragraph, can potentially make the problem much more
diffusion equations,187 none of which is Hamiltonian. Note that their
expensive. Quite impressively, they actually demonstrated this method
algorithm produced a quantum encoding; hence, they did not consider
for simulating integrable dynamics on the two-torus, i.e., a surface of
the complexity of measuring the output state nor did they compare
constant action with two phase angles, on the IBM quantum system
this complexity to that of Monte Carlo algorithms. Reference 34 has
one with four qubits per dimension, i.e., eight qubits for 256 states.
worked to improve this approach.
Carleman developed a linearization approach for embedding
C. Semiclassical representations classical dynamics188 just after Koopman’s seminal publication, which
1. Koopman–von Neumann, Koopman–van Hove, is today known as Carleman embedding or Carleman linearization.
and RKHS Carleman linearization is used today by many in the dynamical sys-
tems and control theory community for finding sparse approximations
The natural setting for linearization is the Liouville equation for
of nonlinear dynamical systems. The Carleman method has the advan-
the evolution of a classical PDF on the classical phase space.29 For
tage that, for linear systems of equations, the embedded linear system
Hamiltonian dynamics, the phase space flow is divergence-free and
exactly reproduces the eigenvalue spectrum of the Koopman evolution
the evolution operator, called the Perron–Frobenius operator, that is
operator. If the location of a fixed point of interest is known, it is useful
generated by the Liouville equation is unitary.177,178 For phase-space
to linearize a nonlinear dynamical system around a given fixed point.
flows that are not divergence-free, one must be careful to distinguish
between the evolution of a scalar function and the evolution of a Then, applying the Carleman approach ensures that the linear dynam-
volume-form or density. In order to determine a unitary evolution ics near the fixed point is faithfully reproduced.
operator, one must consider the evolution of the square-root of a vol- However, unless the original system is exactly linear, a truncation
ume form, called half-forms by the geometric quantization commu- of the linearized system will not generically converge to the evolution
nity. If there is no evolution of the phase of the half-form, then this operator in other regions of phase space nor will the linearized dynam-
evolution equation is called the Koopman–von Neumann (KvN) equa- ics capture all of the topological features of the original system. For
tion. Reference 29 referred to generalizations of the KvN equation that example, nonlinear systems often have multiple fixed points, but any
include the evolution of the complex phase as generalized KvN equa- finite-dimensional linearized system can only have a single fixed point.
tions. In fact, for Hamiltonian dynamics, there is one particular evolu- Thus, the linear system can never display truly chaotic dynamics, such
tion of the phase that corresponds to semiclassical dynamics called the as trajectories that are attracted to multiple fixed points. These issues
Koopman–van Hove (KvH) equation,182,183 as discussed in Sec. V B 2. are discussed further in Sec. V D.
Integrating the KvN or KvH equations in time leads to the gener- The Carleman approach was connected to Hilbert space methods
ation of a nonlinear map. Reference 30 explained that the technique and further generalized by Steeb189 and Kowalski190 and their
described above also applies to such maps in a similar sense that the approach is clearly described in Refs. 191 and 192. Carleman embed-
mapping between phase space functions is linear and unitary. ding is closely related to the use of number states and coherent states,
Concrete algorithms that fully exploit this method are still in i.e., the Segal–Bargmann space, which represent a type of complex-
their infancy. First, Engel et al.23 determined that there is a speedup analytic RKHS, defined over the entire complex plane. For this tech-
for solving the linearized Vlasov–Poisson system. They used the nique, assume there is a destruction operator a and a creation operator
method of qubitization22 to determine an efficient encoding of the a† that satisfy the CCR ½a; a† ¼ 1. In addition, assume there is a vac-
Hamiltonian. uum state j0i annihilated by the destruction operator a and that there
Next, Ref. 29 explained the general approach of simulating the is no subspace that is invariant under the action of a† . Then, there is a
KvN equation. However, Ref. 29 did not advocate a specific discretiza- countable infinity of orthogonal states jki ¼ ða† Þk j0i that are unnor-
tion for the PDF nor did it explicitly specify algorithms for the state malized versions of the usual number states, defined by the number
preparation required for the initialization and output subroutines. operator N ¼ a† a. The unnormalized coherent states are defined via
†
While it is clear that finite difference, finite volume, and finite element jzi ¼ eza j0i, and, can easily be shown to be eigenstates of the
approaches to discretization all lead to sparse representations of the destruction operator, ajzi ¼ zjzi. If the eigenvalues satisfy the
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-24

C Author(s) 2023
V
nonlinear equations of motion, dz=dt ¼ Vðz; tÞ, then the coherent wavefunction as a geometric object that represents the square root of a
states satisfy the general linear equations of motion190–192 real volume form, and, because this object corresponds to a quantum
pure state, the evolution is unitary. This implies that a quantum simu-
@t jzi ¼ Ajzi A :¼ a† Vða; tÞ: (91) lation can directly apply the techniques of Hamiltonian simulation.
In contrast, when working with a complex analytic RKHS, it is
In fact, due to the completeness of the coherent states, any complex desirable to ensure that the evolution remains complex analytic so that
analytic function, w, can be written as a linear combination of these it stays within the chosen RKHS. Thus, as shown in Appendix B, the
states. Due to the Stone–von Neumann theorem, one can make the Carleman embedding considers the evolution of a complex analytic
identifications a ¼ z and a† ¼ @z . Hence, complex analytic func- wavefunction that undergoes holomorphic dynamics. The direct pro-
tions of z satisfy the equations of motion jection of the complex-analytic KvN law to real space corresponds to
the evolution of a PDF rather than as a density. If the velocity, V, is
@t wjz0 ¼ r ðVðz; tÞwÞ: (92)
not divergence-free, then the evolution law for scalars and for volume
forms clearly differs from that of a half-form by terms proportional to
In other words, for this approach, the states appear to evolve as a PDF
r V. Hence, after projection, the resulting linear differential equa-
rather than as a wavefunction, as assumed for the KvN approach. As
tions may not be Hermitian and the resulting evolution operator may
explained in Appendix B, this is in fact the natural evolution law for a
complex analytic wavefunction on a holomorphic RKHS. not be unitary. In such cases, a quantum simulation must use the
Another important difference is that the evolution occurs on the quantum linear differential equation solver.
overcomplete space of coherent state eigenvalues. More generally, one
could use a potentially overcomplete multi-dimensional feature map D. Limitations of linear embedding
and assume that the dynamics takes place on the feature map rather 1. Finite-dimensional truncation
than evolving directly on the space of observables.
In order to generate a discrete set of equations, the next step is to The dynamical systems community has put much effort into
project onto a set of complex analytic functions. For example, one could understanding how to efficiently estimate the underlying dynamics
directly use the coherent states, which are overcomplete. This leads to experienced by a system from measurements along. For example, the
many different potential choices of bases. Another natural choice is to sparse identification of nonlinear dynamics (SINDy) algorithm173,174
† has become quite popular. The goal of these efforts is typically cast in
pffiffiffiffi states generated by N ¼ a a, which leads to the
use the number pffiffiffimono-
mials, z k = k!. In the position representation, q ¼ ða þ a† Þ= 2, these the following form: embed the observed data into a dynamical system
basis functions become the Hermite polynomials, Hk ðqÞ. that is almost linear. While the majority of the variables evolve linearly
Another interesting RKHS is the standard Bergman space of in response to the other variables, in order to exactly represent the
complex analytic functions inside the unit disk in the complex plane. nonlinear dynamics, a subset of variables must still experience nonlin-
This naturally generalizes to the unit disk within an n-dimensional ear equations of motion. When the SINDy algorithm is successful, the
pffiffiffiffiffiffiffiffiffiffiffi
complex space. The Bergman basis uses the monomials z k k þ 1. nonlinear terms are typically small for the majority of the evolution,
Using a standard conformal transformation called the Cayley trans- but they always eventually kick in within certain regions of phase space
form, this naturally maps to complex analytic functions on the upper to provide the nonlinearities that are necessary for truly chaotic
half plane (where the unit disk maps to the real line) and as well as to dynamics.173
the functions with positive real part. The generalized Bergman spaces The present version of quantum algorithms for nonlinear
are also defined by a weight function on the unit disk, and this general- dynamics is actually of a different character, where only exactly linear-
ization is described in Appendix C. ized systems are evolved. Unfortunately, a finite linear system can
The case where the weight function only has support on a curve never display truly chaotic dynamics, such as trajectories that are
in the complex plane, is called a Sz€ego kernel. In fact, the Carleman attracted to multiple fixed points. For the Carleman approach, there is
method uses the Hardy space of complex analytic functions inside the the benefit that the exact linearized dynamics near a fixed point can be
boundary of the unit disk, where the weight function lies entirely on followed exactly. However, because a linear system can only have a
the unit circle along the boundary. The Carleman basis uses the mono- single fixed point, trajectories cannot be attracted to other fixed points,
mials z k . Again using the Cayley transform, this can be conformally nor can there be changes in which basin of attraction a chaotic trajec-
mapped to a region of the upper half plane where the weight function tory is orbiting. Clearly, this is qualitatively quite different than truly
lies entirely along the real line. nonlinear dynamics, in which chaotic dynamics and bifurcations are
Engel et al.32 explored alternate approaches for generating new ubiquitous.
bases of functions for Hilbert space by exploring alternate choices of
operators that satisfy the CCR. As proven in Appendix C, this freedom 2. Domain of convergence
corresponds to exploring other choices of complex analytic RKHS.
Exactly which bases are more convenient likely depends on the appli- For all numerical methods, which are intrinsically finite-
cation of interest. dimensional, only a finite region of phase space can be described with
high accuracy. Clearly, both the initial conditions and the subsequent
evolution must fit within the domain. For example, for the discrete
3. Which form is best?
Fourier representation, the domain can be considered periodic. Note,
The approach advocated by Ref. 29 directly uses one of the gener- however, that if the actual region of interest is not periodic, then unde-
alized KvN representations, e.g., standard KvN or KvH. This treats the sirable interference effects will occur near the boundary.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-25

C Author(s) 2023
V
In contrast, for the case of the complex analytic Hilbert spaces, amount of information that needs to be passed between quantum-
e.g., the standard Bergman and Hardy spaces described above, the classical interface, namely, the classical inputs to the quantum algo-
domain is the inside of the unit disk in the complex plane. In fact, the rithm and the output measured from the quantum simulation, is rela-
solution will blow up if it attempts to cross outside the boundary of tively small.
the domain, i.e., the unit disk. Thus, in order for the computation to An important example of this is the ideal N-body problem in
remain valid, it is important for the solution to avoid traveling outside astrophysics, plasma physics, and molecular dynamics. In this case,
the boundary. Clearly, such a rescaling needs to be applied so that ini- the electromagnetic and/or gravitational fields can be considered the
tial conditions lie within the domain.32 Moreover, one must estimate reduced order model since it only has the three dimensions of configu-
the maximum excursion of the trajectory from the origin and rescale ration space. The linear evolution of the classical PDF/quantum wave-
the equations in a manner that ensures that the maximum excursion function can be considered to be the linear model for evolution on the
lies within the boundary. See the remarks at the end of Appendix C very high 6N-dimensional phase space/3N-dimensional configuration
for more discussion of these issues. space of the particles.
In fact, this is precisely the way that quantum co-processors are
3. Dissipation vs information scrambling expected to be used to accelerate calculations in chemistry, materials
science, and quantum hydrodynamics. The macroscopic degrees of
The KvN/KvH approaches present the dynamics as explicitly freedom could be simulated by a classical computer, while the micro-
occurring in phase space. Because the dynamics is unitary, there are scopic degrees of freedom could be simulated by a quantum computer.
no true Lyapunov exponents for the wavefunction itself. The basis If the microscopic degrees of freedom are assumed to respond linearly
functions of the KvN/KvH approaches are oscillatory functions are to the macroscopic degrees of freedom, as in an analog of the
analogous to the Fourier basis. Consequently, any truncation of the Born–Oppenheimer approximation, then the quantum computer
basis at finite order can potentially introduce Gibbs phenomena and could efficiently simulate their dynamics. Assuming the combined sys-
spurious oscillations. Yet, there is still a possibility of observing truly tem is sparse, the hybrid simulation will remain efficient even if there
chaotic behavior due to the fact that the dynamics does follow the clas- are exponentially more microscopic degrees of freedom in the high-
sical trajectories. order model.
However, there are limitations for modeling dissipative systems:
if information cannot be dissipated, then there are situations where VI. QUANTUM COMPUTING APPLICATIONS
the information must eventually become scrambled. In fact, Ref. 30 There are a number of quantum software simulators and actual
noted a certain scrambling of information that appears to generically hardware platforms available today for use through both commercial
occur for finite difference approximations of the KvN embedding for and government funding sources. This has led to a number of impres-
nonlinear maps. A similar phenomenon was seen in numerical studies sive demonstrations that now claim to have achieved quantum advan-
of Ref. 33 that targeted the dynamics near a stable sink. These authors tage or even quantum supremacy.10–12
introduced another approach that they called the chemical master In the literature, the phrase quantum simulator and quantum
equation approach, but this method had similar issues. emulator are often used to denote hardware platforms that can emu-
The scrambling effect appears to be related to the structure of the late the behavior of a given quantum mechanical system. Here, we
KvN/KvH eigenfunctions for standard finite difference approxima- shall use the phrase classical simulator to denote the use of classical
tions. Perhaps there are alternate numerical discretizations that would computer software to simulate a quantum mechanical system. Clearly,
solve this problem, but it seems likely that such discretizations would such classical simulators are less efficient than a future error-corrected
need to trade scrambling for dissipation, e.g., in the form of upwinding quantum computer should be. However, they are also much more effi-
the finite difference approximation of the advection operator. If dissi- cient than today’s quantum computing platforms. Today, using a clas-
pation is introduced at the numerical discretization stage, then one sical simulator it is relatively easy to simulate up to 15 qubits on a
would need to use the linear differential equation solver rather than laptop and up to 50 qubits on a supercomputer. Yet, in practice, it is
Hamiltonian simulation. difficult to perform high-precision calculations with more than a
handful of physical qubits on present-day quantum hardware.
E. New hybrid algorithms: Nonlinear classical For the classical dynamics application space, there have been a
dynamics coupled to linear quantum dynamics number of notable recent quantum software simulations including: the
linearized Vlasov equation,23 the diffusion equation,156 the Burgers
Perhaps another concept in hybrid classical-quantum dynamics
equation,31 and the linear mode conversion problem.130 In contrast,
would be interesting to study. Separate the nonlinear system into a
there are few demonstration of realistic algorithms on experimental
nonlinear classical reduced-order model with relatively few degrees of
quantum computing platforms. Notable experimental results include
freedom and a linearized quantum high-order model. The high-order
simulations of integrable classical dynamics,163 non-integrable quantum
model represents an effective bath that the reduced order model expe-
maps,193–195 which can be considered models of wave–particle interac-
riences. For example, the SINDy type algorithms discussed above can
tions, simulations of quantum wave–wave interactions,196 and simula-
sometimes represent the system as being completely linearized except
tions of the Dirac quantum walk for quantum hydrodynamics.197
for a single variable that depends nonlinearly on the overall system.
A classical computer can efficiently evolve the nonlinear reduced
dynamics, while a quantum computer can efficiently compute the lin- A. Quantum hardware platforms
ear bath dynamics. Assuming that the high-order model has exponen- Anyone can get started exploring laws of quantum mechanics
tially more degrees of freedom than the reduced-order model, then the by simulating their first quantum algorithms using the open-access
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-26

C Author(s) 2023
V
IBM-Quantum (IBM-Q) cloud-based superconducting quantum com- surfaces. Due to the effects of quantum interference, the insulating
puting platform. Today, this gives free but limited access to up to five regime persists even after the classical transition to chaos, namely, the
qubit devices. There are also many other companies developing quan- quantized map will still exhibit dynamical localization, even when the
tum hardware platforms including D-Wave, Google, Intel, IonQ, corresponding classical map is chaotic. Localization occurs when the
Microsoft, Quantinuum (Honeywell), Rigetti, Xanadu, etc. number of quantum states in the chaotic sea is too small. Since the
From a scientific perspective, having good experimental access to number of states roughly scales as the classical action in chaotic sea
the hardware platform and enough runtime to characterize the perfor- divided by h, the localization effect occurs when Planck’s constant is
mance of the system is crucial for improving performance. For this too large relative to the physical value.
reason, the U.S. Department of Energy’s (DOE) Advanced Scientific Dynamical localization was eventually realized to be closely
Computing Research (ASCR) program supports the activities of sev- related to Anderson localization in condensed matter theory.202 Thus,
eral research groups, called Quantum Testbeds for Science, for provid- in the generic case, the wavefunction is localized with exponential
ing open science, quantum hardware platforms for research on tails. Chirikov et al.201 realized that the effective localization length
experimental quantum computing architectures. can be estimated from a simple scaling argument: the number of
For superconducting qubits today,87 the typical gate times are momentum states included in the wavefunction is roughly the num-
10 100 ns while the typical relaxation and dephasing times are ber that can be reached by chaotic diffusion up to the Heisenberg
10 100 l s. This implies a minimum effective error rate on the time: ðh‘c Þ2 ¼ Dc sH . Estimating sH =s ’ Emax =DE ’ ‘c , then yields
order of 0.1% per gate application. Measured two-qubit entangling the estimate ‘c ’ Dc s=h2 .
gate error rates are typically ten times higher than this, on the order of The lack of a true Lyapunov exponent in quantum dynamics has
1%. Clearly, this implies the existence of additional sources of error, made the regime of quantum chaos more subtle to define than its classi-
such as imperfect gate design and crosstalk between neighboring cal counterpart. While the quantum-classical correspondence principle is
qubits. For ion-trap platforms,86 the gate times on the order of 1–100 known to hold for integrable systems, the meaning for chaotic systems
ls. Although relaxation and dephasing times are quite long when is more subtle and questions still remain. Clearly, for chaotic systems,
qubits are in isolation, 10–100 s, errors are introduced when gates are the correspondence principle only holds up until the Heisenberg time.
enacted and typical gate fidelities are better than but still comparable In fact, in the presence of noise, chaotic quantum dynamics is much eas-
to superconducting platforms. ier to reverse than the classical counterpart. After reversal, classical tra-
jectories will eventually begin to separate at the Lyapunov rate if there is
B. Wave–particle interactions any uncertainty present, but the quantized version will never do so
Some of the first quantitative studies of classical chaos within indefinitely. In other words, the quantum version often has a much
plasma physics were focused on wave–particle interactions. This led to stronger Loschmidt echo, i.e., overlap between the initial and final states,
highly simplified mathematical dynamical systems called nonlinear than the classical version. Although the relevant Heisenberg time might
maps that displayed all of the features of chaos in a form that was be longer than the lifetime of the universe for realistic values of Planck’s
numerically precise and independent of additional approximations, constant, this qualitative difference may be somewhat surprising.
e.g., introduced by numerical integration and/or nonlinear solvers. Peres realized that there is still a kind of quantum-classical corre-
Quantized versions of these classical nonlinear maps, called quantum spondence in the dependence of the Loschmidt echo on the dynam-
maps or Floquet systems, were developed for similar purposes: under- ics.203 For both the classical and quantum cases, the coarse-graining
standing quantum chaos in highly simplified yet precisely defined sys- introduced by lack of a precise knowledge of the exact system
tems that were relatively easy to simulate. Hamiltonian is responsible for increasing entropy with time.
In the 1970s and 80s Giulio Casati, Boris Chirikov, and Dima Eventually, it was discovered204,205 that the Loschmidt echo has three
Shepelyansky made some of the first fundamental discoveries about different regimes: a perturbative regime of slow quadratic decay,203 a
quantum chaos by studying quantum maps.199–201 They found that, if Fermi golden rule regime often observed in the microscopic quantum
the quantum map had chaotic dynamics, the wavefunction would rap- world,205 and a semiclassical regime where it is controlled by the classi-
idly spread out along the unstable manifolds of the semiclassical trajec- cal Lyapunov exponent.198,204,205
tories with a growth rate set by the fastest Lypanov exponent, k. This In fact, the decay of the fidelity of a noisy quantum simulation
leads to a rapid spreading of the wavefunction or density matrix on yields information about the nature of the dynamics being simu-
the Ehrenfest time, sE ¼ log ðNc Þ=k, where Nc is the number of eigen- lated.170 The protocol for measuring the so-called Loschmidt fidelity
states accessible within the chaotic sea. An effective chaotic diffusion echo is to simulate the dynamics for a certain number of time steps
in momentum or action space proceeds until the Heisenberg time, and then to simulate the dynamics backward in time for the same
sH ¼ h=DE, set by the average difference between energy levels; i.e., number of steps. In the presence of noise, the final result will not
2p=DE is the mean density of states at the energy scale of interest. exactly match the initial conditions. When the dynamics is chaotic
Since the quantum dynamics is linear, for time scales much longer and diffusive (conducting), as shown on the left of Fig. 10, then the
than the Heisenberg time, sH, one can observe discrete eigenfrequen- wavefunction spreads out over a large region of phase space. In this
cies, and the actual Lyapunov exponent of the quantum dynamics case, the fidelity decays exponentially, with a decay rate set by the
must vanish. Lyapunov exponents of the system.204–206 When the dynamics is
In the chaotic regime, where the dynamics is diffusive/conduc- approximately integrable and localized (insulating), as shown on the
tive, the typical wavefunction extends over the entire chaotic sea. In right of Fig. 10, then the wavefunction only spreads out over a surface
contrast, in the integrable regime, where the dynamics is localized/ of constant energy (action). In this case, the fidelity decays algebrai-
insulating, the wavefunction can only tunnel slowly through KAM cally as t d=2 where d is the dimension of phase space.198
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-27

C Author(s) 2023
V
FIG. 10. Here, the Loschmidt echo protocol measures the fidelity of the quantum sawtooth map [Eq. (93)] after simulating the dynamics for t steps forward in time and t steps
backward in time. (Left) If the dynamics is chaotic/diffusive/conducting, as in the sawtooth map at K ¼ 0.5, then the wavefunction spreads out over a large region of phase
space and the fidelity decays exponentially at the Lyapunov rate (center). (Right) If the dynamics is integrable/localized/insulating, as in the sawtooth map at K ¼ 0.1, then the
wavefunction only spreads out over the energy (action) surface and the fidelity decays algebraically (center).194,198
For Hamiltonians with separable kinetic and potential energies, a The classical sawtooth map (CSM) is defined by the Lagrangian
quantum computer can use QFT’s to go back and forth between posi- X
tion and momentum space exponentially faster than a classical com- LCSM ðpj ; qj Þ ¼ pjþ1 qj HSM ðpj ; qj Þ; (94)
j
puter. The only remaining question is whether the remaining phase
shifts associated with the kinetic and potential energies can be simu- where qj and pj are the coordinates at time t ¼ js. Hence, the CSM is
lated efficiently. In fact, for many energy functions, such as polyno- explicitly given by
mials or trigonometric functions, the phase shifts can be performed
efficiently.207 Hence, Benenti et al.144,170 and Georgeot and qj qj1 ¼ pj s; (95)
Shepelyansky207 proposed that this quantum factorization can be used pjþ1 pj ¼ Kqj s: (96)
to calculate many properties of the dynamical systems more quickly
on a quantum computer. For example, in the localized regime, this Despite the simplicity of the CSM, it still has very rich dynamics. It is
algorithm can speed up the calculation of the localization length, and, chaotic for K < 4 and K > 0, integrable for K ¼ 4; 3; 2; 1; 0,
in the chaotic regime, the chaotic diffusion coefficient (or, when trans- and generically displays anomalous diffusion without a positive
port has anomalous scaling, the anomalous transport coefficients) can Lyapunov exponent for 4 < K < 0. For any Hamiltonian that has a
also be measured with polynomial speedup.144,170 They also proposed quadratic kinetic energy and a periodic potential energy, the dynamics
that simply measuring the fidelity decay can be used to calculate the is periodic on p over the interval p ps < p.
Lyapunov exponent with polynomial or even superpolynomial The quantum version, the QSM, is defined by promoting the
speedup.170,206 coordinates fq; pg to operators fq; pg that satisfy the canonical com-
The idea that one can use a noisy quantum computer to measure mutation relation: ½q; p ¼ ih. Since there are no terms that mix the
the classical Lyapunov exponent is of profound importance to the clas- q’s and p’s, there is no operator ordering issue and there is a natural
sical computations of interest to FES. Following in their footsteps, we quantization of the sawtooth Hamiltonian given by HQSM
studied the sawtooth map, which is one of the easiest maps to simulate :¼ HSM ðq; pÞ. Because the quantum dynamics must be unitary, the
on quantum hardware and developed efficient implementations of the quantum sawtooth map evolution operator is then defined to be
map for the IBM-Q platform.194,195 The sawtooth map is relatively UQSM ¼ exp ðiHQSM s=hÞ.
easy to simulate because it uses a quadratic form for both the kinetic A finite discretization must have a finite number of degrees of
and potential energies freedom; e.g., N discrete points along q. For this choice, the momen-
tum operator has discrete eigenvalues, and is periodic with a range of
1 1 X twice the Nyquist frequency hN. Interestingly enough, for any map
HSM ðp; q; tÞ ¼ p2 Kq2 dðt jsÞ (93)
2 2 j
periodic in p, in order for the classical periodicity of p over the range
of 2p=s to consistently match the quantum periodicity of p on hN,
with time step s and kicking strength K. The potential energy is one must then choose the specific value for Planck’s constant
defined as above over the interval p q < p and is assumed to be h ¼ 2p=Ns. Because the Hamiltonian can be considered to consist of
periodic outside this interval, so that it is continuous but not smooth. separate applications of the kinetic and potential energy operators,
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-28

C Author(s) 2023
V
each time step naturally divides into four steps: QFT, phase shift, the highly entangled wavefunctions that are naturally generated by
inverse QFT, and another phase shift chaotic dynamics are quite sensitive to dephasing. The noise that gen-
erates dephasing acts to randomly alter the phases of each term in the
UQSM ¼ Ukinetic QFT† Upotential QFT: (97) wavefunction, and, hence, generates more fidelity loss when the wave-
Again, in this case, there are efficient algorithms for performing the function is extended over many states. Fitting the Lindblad model to
phase shift operations because each term is only quadratic in momen- three different types of noisy quantum circuit simulations consistently
tum/position operators. found that gate operation appears to induce a 3 increase in the
Figure 9 compares the result of simulating the quantum and clas- dephasing rate relative to the rates measured while the qubits are idle.
sical versions of the map in the regime of anomalous diffusion where If the gate-based Lindblad model is an accurate model of reality, then,
K ¼ 0:1. This choice is interesting because there is a lot of structure unfortunately, this pushes the operating window for measuring the
in phase space that is not quickly destroyed by chaotic diffusion. On Lyapunov regime to even higher qubit numbers.195
the right is a classical Poincare plot of the result, starting all trajectories
from p=p ¼ 3=4 and uniformly sampling the phase. On the left are C. Wave–wave interactions
the quantum results, shown by the Husimi-Q quasiprobability distri- Wave–wave interactions are important nonlinear processes in
bution function. In this case, the wavefunction is started in the pure plasma physics, fluid dynamics, nonlinear optics, and quantum field
momentum eigenstate with p=p ¼ 3=4. As the number of qubits is theory. For a general dynamical system in a weakly coupled regime,
increased from n ¼ 6 to 9 to 16, the Husimi-Q function begins to decomposing nonlinear dynamics into wave–wave interactions is a
resemble the Poincare plot more and more closely. ubiquitous method for perturbative analysis and can also yield highly
Our theoretical investigations showed that the QSM requires nontrivial non-perturbative predictions, as in the study of wave turbu-
more than five qubits to measure fidelity decay at the Lyapunov rate194 lence.208,209 When the system fluctuates near its fixed points with small
and our experimental investigations imply that much lower error rates amplitudes, the linearized system possesses a spectrum of eigenmodes,
are needed than are observed at present.195 The operating window for which are often called linear waves for hyperbolic systems. The linear
observing Lyapunov decay is limited for three reasons: (i) the need for waves form a convenient basis of the system, and the Hamiltonian is
the Lyapunov rate to be slower than the decay rate initially induced by determined by a rank-2 Hermitian tensor, which represents a second-
noise, so that an intermediate time asymptotic appears where semi- order coupling between waves. At larger wave amplitudes, nonlinear-
classical effects are observable, (ii) the need for the dynamics to be in ities manifest in perturbative analysis as higher order couplings
the diffusive/conducting chaotic phase, rather than the dynamically between waves. To lowest order, the nonlinear coupling typically
localized/insulating phase, and (iii) the need for the fidelity decay rate involves three waves, leading to what are known as three-wave interac-
to be slow enough to be able to measure multiple timesteps in the tions. Near stable fixed points, three-wave interactions are of decay
semiclassical intermediate time window. Based on our experimental type, where one pump wave decays to two daughter waves, or, con-
investigations on the IBM-Q superconducting platform as well as on versely, two daughter waves merge to excite the pump wave. Near
existing experimental data on the IonQ ion trap platform, we esti- unstable fixed points, three-wave interactions are of explosive type,
mated the reduction in error that would be required to measure the where the three waves are either spontaneous generated or destroyed.
Lyapunov exponent depending on potential hardware and software Couplings of higher order, such as four-wave coupling, emerge in
capabilities. In the best case scenario of all-to-all parallelization, the next-to-leading order perturbative analysis, giving rise to N-wave
gate depth is significantly reduced and then errors might only need to interactions.
be reduced by a factor of a few. However, in the worst case scenario of For our first attempt at quantum simulation, we investigated how
linear QPU topology, errors may need to be reduced by a factor of 100 to simulate three-wave interactions on quantum computers.196 The
or more.194 three-wave interaction is a nonlinear process, which in its classical
Hence, we were surprised to find that the experimental fidelity of form is difficult to map onto quantum hardware using embedding
quantum simulations on IBM-Q hardware using only three qubits methods discussed earlier. As an alternative approach discussed in Sec.
does in fact depend on the dynamics.195 In fact, relative to the IBM-Q V A, we chose to simulate a quantized version of the problem, which
reported error rate based on randomized benchmarking, the fidelity approaches the classical interaction in the limit of large photon num-
decays 3 faster for integrable dynamics and 5 faster for chaotic ber, where the statistics are typically Poissonian. This version is imple-
dynamics. While it is well known by experts in the field that such mentable on current quantum hardware platforms. For decay type
effects are possible, the most commonly used methods of characteriz- three-wave interactions, both classical and quantum processes are
ing quantum hardware platforms, randomized benchmarking, can described by what are known as the three-wave equations
only be used to estimate the effects of a simple depolarizing noise
model. Since the depolarizing noise model only generates a combina- da1 =dt ¼ þga2 a3 ; (98)
tion of the original density matrix, q, and an incoherent mixture, da2 =dt ¼ g a1 a†3 ; (99)
qincoherent ¼ 1=N, it can never display more than a single decay rate. da3 =dt ¼ g a1 a†2 ; (100)
Yet, for three qubits, this effect also does not appear using the unitary
noise models studied by Refs. 194 and 206. where d=dt ¼ @t þ V r is the convective derivative of waves at the
Instead, we discovered that a gate-based Lindblad model, in group velocity, V, and g is a complex-valued coupling coefficient. In
which single-qubit relaxation and dephasing rates are enhanced during the quantum problem, the aj are bosonic operators that satisfy canoni-
gate operation, can potentially explain this effect.195 The reason is that cal commutation relations ½aj ; a†k ¼ djk , whereas in the classical prob-
highly localized wavefunctions are insensitive to dephasing, whereas lem, these operators are reduced to commuting c-numbers,
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-29

C Author(s) 2023
V
X
½aj ; a†k ¼ 0. In fact, the coupling coefficient g for classical three-wave jwi ¼ wj js3 j; s2 s3 þ j; ji: (106)
interactions may be more easily derived from quantum theory.210 s2 s3 ;j
The action of the jth wave is proportional to the number operator
Therefore, it is sufficient to focus on each of the subspaces indepen-
nj ¼ a†j aj . For both the quantum and classical problems, the relative
dently. Moreover, since s2 ; s3 0, the values of j are constrained by
action invariants
0 j minðs2 ; s3 Þ. In other words, each subspace is finite-
S2 ¼ n1 þ n2 ; (101) dimensional. Notice that the total dimension of the problem is infinite,
S3 ¼ n1 þ n3 (102) but the system is a direct sum of disjoint irreducible subspaces each of
finite dimension. This special property makes artificial truncation
are conserved. For the quantum problem, these operators commute unnecessary, thereby avoiding the issues discussed in Sec. V D 1.
with the Hamiltonian and each other, while, for the classical problem, Notice that the special property only holds for decay type wave–wave
these quantities are conserved by the Hamiltonian and are in involu- interactions near an elliptic fixed point. For explosive type interactions
tion with each other. Because the operators Sj commute with the num- near a hyperbolic fixed point, the entire Hilbert space will become
ber operators nj , they can be simultaneously diagonalized, with connected.
eigenvalues sj and nj, respectively. In the Heisenberg picture, the sj are In each subspace, the nonlinear three-wave problem is mapped
constants of the motion, but the nj evolve in time. to a linear Hamiltonian simulation problem, where the Hamiltonian
The difference between the quantum and classical versions of the matrix is tridiagonal with zero diagonal elements. Explicitly, the
problem becomes manifest when computing observables. The quan- Schr€odinger equation ih@t jwi ¼ H3wave jwi becomes a system of lin-
tum problem satisfies ear equations for the expansion coefficients cj ðtÞ, which satisfies
@t2 hn1 i ¼ @t2 hn2 i ¼ @t2 hn3 i ih@t wj ¼ igHjþ12 cjþ1 ig Hj12 wj1 : (107)
2
¼ 2jgj s2 s3 ð2s2 þ 2s3 þ 1Þhn1 i þ 3hn21 i : This problem may be thought of as quantum walk on a one-
(103) dimensional lattice, where the quantum state hops from one lattice site
to its two nearest neighbors. What distinguishes this quantum walk
In contrast, the classical problem satisfies from others is the specific matrix element
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
@t2 hn1 i ¼ @t2 hn2 i ¼ @t2 hn3 i ¼ 2jgj2 s2 s3 2ðs2 þ s3 Þhn1 i þ 3hn1 i2 ; Hj12 ¼ jðs3 þ 1 jÞðs2 s3 þ jÞ; (108)
(104)
which encodes the nonlinear three-wave interaction. Notice that
where, now, the bracket denotes the expectation value averaged over a the transition matrix elements satisfy H12 ¼ HDþ12 ¼ 0, where
classical ensemble of initial conditions. Hence, the quantum equation D ¼ minðs2 ; s3 Þ þ 1 is the dimension of the problem. In other words,
of motion is almost the same as the classical equation except for the the quantum walk has a compact support, and the quantum state is
final two terms. First, notice the extra 1 in the parentheses of the mid- confined within a lattice of a finite number of sites. On ideal quantum
dle term Eq. (103). This term is due to spontaneous emission, which computers, the Hamiltonian simulation problem may be solved effi-
only occurs in the quantized version. However, it is subdominant in ciently using algorithms discussed in Sec. IV. However, on current
the large photon number limit s2 ; s3 1. Second, the variance of the noisy devices, we do not yet have the capability of implementing the
final term hn21 i often differs between the quantum and the classical necessary oracles.
versions, precisely because of the difference in statistics implied by the Instead of performing quantum Hamiltonian simulation, we
density matrix vs the classical PDF, as well as because of the different tested current quantum hardware platforms by performing a simple
evolution equations that the higher moments satisfy. For example, unitary evolution that can be represented using only two qubits. For a
for a single classical trajectory with a particular value of nj, then time step Dt, the unitary evolution under the time-independent
hn2j i ¼ hnj i2 ¼ n2j (and the same result holds true for any power of Hamiltonian H is given by the matrix exponential UðDtÞ
nj). This relation between the first two moments is similar to the statis- ¼ exp ðiHDt=hÞ. For small problem size, the exponential may be
computed analytically, and for larger problem sizes, the exponential
tics generated by the Poisson distribution.
can be calculated using classical computers. Of course, the point of
The three-wave Hamiltonian
quantum computing is that for very large problem sizes, the unitary
H3wave ¼ iga†1 a2 a3 ig a1 a†2 a†3 (105) can be realized to arbitrary precision without already knowing the uni-
tary matrix classically. However, for the few qubits that are available
governs the three-wave interaction, because its Heisenberg equations on present-day hardware, the problem size is not too large, so the
of motion are precisely the three-wave Eqs. (98)–(100). The quantum short-time propagator UðDtÞ can be computed classically. Then, the
problem can be readily solved in the Schr€odinger picture once a basis propagator is compiled into quantum gates, and the quantum com-
is chosen. A convenient basis is formed by eigenstates of the relative puter is responsible for repeating the gates to carry out the time
action invariants, S2 and S3 , in Eqs. (101) and (102). The basis is con- advance UðDtÞN for N time steps.
venient because these operators commute with both the Hamiltonian In this problem, compounding the unitary is somewhat
and the number operators. Since S2 and S3 are constants of motion, trivial because UðDtÞN ¼ UðNDtÞ can be fast forwarded. However,
the dynamics in subspaces with different eigenvalues s2 and s3 are this is a good benchmark for more complicated problems where
decoupled. Denoting the Fock states by jn1 ; n2 ; n3 i, then any quantum the total Hamiltonian is composed of non-commuting terms, such
state can be represented by as H ¼ HA þ HB . In this more typical non-commuting case,
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-30

C Author(s) 2023
V
each time step. Consequently, the unitary evolution can only reach a
depth of about 10, as we have observed in the test run.
In order to increase the simulation depth, we performed unitary
evolution of the three-wave problem on LLNL’s quantum design and
integration testbed (QuDIT) platform,65 which gives users white-box
access to control pulse design for a transmon qutrit. Instead of prepar-
ing a universal gate set and then compiling the desired unitary in terms
of these gates, we prepared a customized gate that directly realizes the
unitary as a single control pulse. In this way, each simulation time step
is realized by only one customized gate, instead of 20 universal gates.
At constant device performance, the compressed gate depth enables
roughly 10 longer effective simulation depth.
Optimal control uses standard optimization techniques, which
were first applied to manipulate quantum systems in the nuclear mag-
netic resonance community213 and were later applied to many other
quantum systems.214,215 For a given hardware Hamiltonian, H0 , and
FIG. 11. Experimental results of simulations of a three-level three-wave interaction control Hamiltonian, Hc , quantum optimal control searches for classi-
problem using different compilation techniques.146,196 On Rigetti Aspen-4-2Q-A (red cal wave form, c(t), such that by the end of the time evolution under
error bars), the unitary for a single time step is compiled into a circuit composed of H ¼ H0 þ cðtÞHc , the final state of the quantum device is mapped
CZ and Pauli gates. In contrast, on the LLNL QuDIT platform (blue dots) each time
step is compiled into a single optimized customized microwave control pulse. The
from its initial state by a unitary that is e-close to the desired unitary.
later approach better matches the exact results (black line), and the deviations are For the three-wave problem, control pulses for UðDtÞ can be readily
explained well by a phenomenological GKLS master equation (ME) model (purple generated using optimize_pulse_unitary in QuTIP,216,217 using H0
line) using experimentally measured decoherence rates. Reproduced with permis- and Hc that are calibrated immediately before the production runs on
sion from Phys. Plasmas 28, 042104 (2021). Copyright 2021 AIP Publishing. the QuDIT platform. Typical control pulses are shown in Fig. 12.
Repeating the single-time step pulses for N times, we achieve uni-
where HA and HB are easy to simulate separately but the total H is dif- tary evolution shown by blue points in Fig. 11. The experimental
ficult to simulate, product formula approximations such as results match the exact solutions much better, and the deviations can
ðUA ðDtÞUB ðDtÞÞN ¼ 6 UA ðNDtÞUB ðNDtÞ cannot be fast forwarded. be explained by a phenomenological decoherence model described
Quantum computers realize unitary evolutions by repeating pre- through simulations of the GKLS master equation (ME). The GKLS
compiled gate sequences, unlike classical computers which compute model is calibrated using experimentally measured decoherence rates
the total unitary by a series of matrix multiplications. We implement of the quantum hardware (and was not calibrated by fitting the three-
unitary evolution for a D ¼ 3 test problem, which can be embedded in wave experimental data). The optimal control technique allowed us to
two qubits, on the Rigetti superconducting hardware platform.211,212 extend the number of simulated time steps by approximately an order
The results are shown in Fig. 11 as red error bars. For comparison, the of magnitude, so that about five time periods could be simulated over-
exact answers are shown by the black line. As can be seen from the fig- all in this test problem. This result constituted the first interesting
ure, only about ten time steps can be performed before results become plasma simulation ever performed on quantum computing hardware.
indistinguishable from noise. Although the dynamics of this system is
integrable, results on current hardware are only able to capture about VII. CONCLUSIONS AND OUTLOOK
half of the nonlinear wave period. This unfortunate result is caused by Quantum information science holds great promise for accelerat-
the fact that hardware error rate is about 1% per gate, so after about ing scientific breakthroughs through the three pillars of quantum sens-
102 gate applications, the error becomes of order unity. For the two- ing, quantum communications, and quantum computing. Fusion
qubit problem, the unitary matrix has 24 ¼ 16 elements, so a brute- energy science (FES) stands to benefit, in particular, from quantum
force compilation to hardware native gates results in about 20 gates for computing, due to the importance of scientific computing in driving
FIG. 12. The optimized control pulse that realizes the single-step unitary UðDtÞ for the three-level three-wave problem on LLNL’s QuDIT platform. Probabilities of the ground,
first excited, and second excited states evolve during the control pulse application. The experimentally observed probabilities on the quantum hardware (dotted line) match the
prediction of master equation (ME) simulations on classical computers (lines). By the end of the control pulses, the desired unitary UðDtÞ is realized.
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-31

C Author(s) 2023
V
discoveries in the field. Many useful quantum algorithms have been ACKNOWLEDGMENTS
developed that could be put to good use including: efficient quantum
Fourier transforms, sparse linear solvers, sparse Hamiltonian simula- The authors would like to thank A. H. Boozer, as well as
positive interactions with the 2022 International Sherwood Fusion
tion, and variational eigensolvers. Intrinsically quantum simulations of
Theory Workshop community, for encouraging the writing of this
interest to FES that may stand to benefit include plasma chemistry,
review paper. This work, LLNL-JRNL-839545, was performed
fusion materials science, relativistic laser–plasma interactions, as well
under the auspices of the U.S. DOE by LLNL under Contract No.
as high energy quantum processes that take place in the most extreme
DE-AC52–07NA27344 and was supported by the DOE Office of
plasma environments in the universe, such as black holes and neutron
Fusion Energy Sciences “Quantum Leap for Fusion Energy
stars.
Sciences” Project No. FWP-SCW1680.
For a long time, it was not known whether intrinsically classical
simulations, of the type usually performed within FES computing
applications, would also stand to benefit in similar ways. Recent work AUTHOR DECLARATIONS
has proven that, in principle, Hamiltonian simulation of classical Conflict of Interest
dynamics can achieve an exponential speedup over Eulerian methods. The authors have no conflicts to disclose.
It is important to point out that classical probabilistic algorithms, such
as Monte Carlo methods, also have the ability to perform an effective Author Contributions
exponential reduction in complexity relative to direct solvers.167 For
generic problems, it is expected that quantum algorithms can achieve Ilon Joseph: Conceptualization (equal); Formal analysis (equal); Funding
up to a quadratic speedup over classical probabilistic algorithms. The acquisition (equal); Project administration (equal); Supervision (equal);
anticipated quantum speedup increases with the dimensionality of the Visualization (equal); Writing – original draft (equal); Writing – review &
problem and the lack of smoothness of the solution. If the overhead, editing (equal). Yuan Shi: Conceptualization (equal); Data curation
such as input and output between classical-quantum interface, and (equal); Formal analysis (equal); Investigation (equal); Software (equal);
physical costs can eventually be made sufficient small, then these Writing – original draft (supporting); Writing – review & editing (equal).
quantum algorithms may one day become the algorithms of choice. Maxwell Porter: Conceptualization (equal); Data curation (equal); Formal
This conclusion opens up the door to improving simulations of classi- analysis (equal); Investigation (equal); Methodology (equal); Software
cal N-body problems, molecular dynamics, and fluid and kinetic (equal); Validation (equal); Visualization (equal); Writing – review & edit-
plasma models that are central to FES applications. ing (supporting). Alessandro Roberto Castelli: Data curation (equal);
At present, multiple quantum computing hardware platforms Investigation (equal); Methodology (equal); Software (equal); Visualization
have developed high fidelity qubits that can enact multiple entan- (equal). Vasily I Geyko: Conceptualization (supporting); Visualization
gling gate operations. Unfortunately, the lack of error correction (equal); Writing – review & editing (supporting). Frank Reno Graziani:
severely limits the performance for practical applications—especially Conceptualization (supporting); Writing – review & editing (supporting).
for classical applications, such as classical logic and deterministic Stephen Libby: Conceptualization (equal); Writing – review & editing
simulations. If supremacy is defined as one computing platform (supporting). Jonathan L. DuBois: Conceptualization (equal); Data cura-
being able to perform calculations that another platform cannot, tion (equal); Formal analysis (supporting); Funding acquisition (equal);
then today is certainly the era of classical supremacy. Hence, much Investigation (equal); Methodology (equal); Resources (equal); Software
of today’s quantum computing applications research is focused on (equal); Supervision (equal); Validation (equal); Visualization (equal).
passive and active error reduction techniques, such as the quantum
optimal control methodology discussed in this work. Unfortunately, DATA AVAILABILITY
the estimated cost of fault-tolerant error correction is large enough that
the first useful quantum applications may need to offer speedups that The data that support the findings of this study are available
are beyond quadratic.91 within the article.
It is plausible that before perfect error correction becomes widely
available, quantum computing hardware could still provide a novel APPENDIX A: HILBERT SPACE PRIMER
and useful platform for the simulation of open many-body quantum
dynamics. It may even be the case that such platforms have significant 1. Hilbert space
potential for near-term quantum advantage. However, there is much
work to be done in order to understand the accuracy with which such Definition: A Hilbert space (HS) is a vector space with an inner
simulations can be actually performed and how this accuracy com- product that is complete under the metric induced by the inner
pares to classical algorithms that run efficiently on high performance product.
supercomputers. Simply verifying that the calculation achieves the • The inner product, denoted hji must be linear, skew-
desired fidelity becomes quite challenging once the complexity of the symmetric, and positive-definite
calculation pushes beyond the limits of classical computers. Finally,
since it is decoherence that controls the information confinement time, h/jaw þ bvi ¼ ah/jwi þ bh/jvi; (A1)
it can be hoped that the grand pursuit of a developing a functioning
quantum computer will lead physicists to gain a much deeper under- hwj/i ¼ h/jwi ; (A2)
standing of decoherence and, ultimately, of the fundamental workings
of the universe. hwjwi 0: (A3)
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-32

C Author(s) 2023
V
Theorem: All finite-dimensional Hilbert spaces are equivalent. Now, consider a complex analytic flow velocity Vðz; tÞ, which
satisfies @z V ¼ 0. The complex analytic scalar field, hðz; tÞ, satisfies
• A finite-dimensional Hilbert space is spanned by a finite set of
the constraint, @z h ¼ 0, so it evolves via
basis elements jji.
• Define the complex conjugate transpose operation hjj ¼ jji† . ð@t þ V @z Þhðz; tÞ ¼ 0: (B5)
• In this basis, the Hilbert Space metric has the representation
However, a volume form cannot be complex analytic because it
hjjki ¼ gjk ; (A4) must carry the weight of the full Jacobian, so it must still satisfy

where gjk ¼ gkj , i.e., as a matrix operator g ¼ g† . @t þ @z V þ @z V pðz; z ; tÞ ¼ 0: (B6)
• The Hilbert space metric, g, must be a positive definite matrix,
i.e., it must have positive definite real eigenvalues. Hence, g is A complex analytic flow keeps the solution complex analytic as
invertible and g1 exists. a function of time, so that zðz0 ; tÞ is independent of z0 , i.e.,
• Dual vectors can be defined that satisfy hjjki ¼ djk . The solution is @z=@z0 ¼ 0. Hence, the Jacobian, gðz; z Þ, on the full space,
† X fz; z g, is the norm squared of the complex analytic Jacobian
hjj ¼ g1 j ji ¼ hjjg1 ¼ gjk1 hkj: (A5) j ¼ @z0 =@z; i.e., the full Jacobian is
k
•
gðz; z Þ ¼ jðz ÞjðzÞ ¼ jjj2 : (B7)
An orthonormal basis exists in which the Hilbert space metric is
unity. The orthonormal basis vectors are defined via Because of this splitting of the Jacobian, one can still define a com-
1=2 plex analytic wavefunction via wðz; tÞ ¼ w0 ðz0 Þjðz; tÞ. Thus, the
j^ji ¼ g jji; (A6)
complex analytic KvN equation
^ 1=2
h jj ¼ h jjg : (A7)
ð@t þ @z V Þwðz; tÞ ¼ 0 (B8)
APPENDIX B: PREQUANTIZATION FOR HOLOMORPHIC corresponds to the evolution law for a complex analytic half-form—
HILBERT SPACES with an appropriate choice for the dynamics of the phase. Note how
it appears to be advected as a volume form on a space with half of
Due to the many important uses of complex analysis, it is nat-
the real dimension of the full system. Although this form of the evo-
ural to try to understand how to apply the Koopman–von
lution is not Hermitian, it is actually a special case of the Hermitian
Neumann (KvN) approach to nonlinear dynamics on complex
evolution law of a “wavefunction,” that is a z volume form and a z
manifolds. In order to pass from real to complex variables, one
scalar field
must combine two copies of phase space. Hence, given two copies

of coordinates for the original phase space, x; y 2 Rd , construct a @t þ @z V þ V @z wðz; tÞ ¼ 0 (B9)
combined list of coordinates, z ¼ fx; yg 2 Rd Rd ¼ R2d . By
defining the complex conjugation operator, : fx; yg ! fx; yg, subject to the initial condition @z w ¼ 0.
which is clearly an involution, the doubled phase space is isomor- Given the complex-analytic form of the KvN evolution law,
the norm of a complex analytic wavefunction is preserved in time
phic to a complexified version of phase space, z ¼ x þ iy 2 Cd .
Now, if the dynamics is to act naturally on the complex space, Ð Ð
hwðtÞjwðtÞi ¼ jw0 ðz0 ÞjðzÞj2 dzdz ¼ jw0 ðz0 Þj2 dz0 dz0 ¼ hw0 jw0 i:
the generator of the motion must commute with the complex con-
jugation operator, . Hence, the equations of motion for z and z (B10)
must be compatible in the following sense:
The complex analytic evolution law is still sensible when a
dz=dt ¼ Vðz; z ; tÞ dz =dt ¼ V ðz; z ; tÞ ¼ Vðz ; z; tÞ: (B1) nontrivial Hilbert space inner product is included. Assume the
inner product is defined by the function gðz; z Þ, which is fixed in
Thus, complex scalar fields evolve via time, via
ð
@t þ V @z þ V @z hðz; z ; tÞ ¼ 0 (B2)
h/jwi ¼ /ðz ; tÞwðz; tÞgðz; z Þdzdz : (B11)
and complex volume forms evolve via
For example, if gðz; z Þ represents a complex analytic coordinate
@t þ @z V þ @z V pðz; z ; tÞ ¼ 0: (B3)
transformation, then it has the form gðz; z Þ ¼ jjðzÞj2 . In this
A wavefunction that is defined as the square root of a volume form case, the inner product simplifies to
should evolve via the complex KvN equation ð
1 h/jwi ¼ /ðz ; tÞwðz; tÞjjðzÞj2 dzdz : (B12)
@t þ @z V þ V @z þ @z V þ V @z wðz; z ; tÞ
2
¼ iLðz; z ; tÞwðz; z ; tÞ=h: (B4) The relevant evolution law can be achieved by defining the diver-
gence operation as
Here, the function Lðz; z ; tÞ introduces the freedom of adding a
non-trivial evolution of the complex phase of the wavefunction. r V :¼ j1 ðzÞ@z ðjðzÞ Vðz; tÞÞ: (B13)
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-33

C Author(s) 2023
V
Then, the complex analytic KvN equation is generalized to X

1
z¼ j j 1ihjjj: (C7)
ð@t þ r V Þwðz; tÞ ¼ 0: (B14) j¼1
Again, this is equivalent to the Hermitian evolution law for a wave- An orthonormal set of number states is defined via
function that is a z volume form and a z scalar field
X 1=2
ð@t þ r V þ V r Þwðz; z ; tÞ ¼ 0 (B15) j^ji ¼ g1=2 j ji ¼ gkj jki: (C8)
k
subject to the initial condition @z w ¼ 0.
Given the assumptions, these states are equivalent to the standard
number states, with the standard destruction a, creation a† , and
APPENDIX C: HOLOMORPHIC REPRODUCING KERNEL number a† a operators. Using the definitions above
HILBERT SPACES AND CANONICAL COMMUTATION
RELATIONS (CCR) X
1
1=2 1=2
w¼ j^jih^kjgjm gm1k ; (C9)
192
Motivated by Kowalksi, Ref. 32 recommended exploring the j;k;m¼0
possibility that operators that satisfy the canonical commutation
relations (CCR) X
1
1=2 1=2
z¼ j^jih^kjgjm mgm1k : (C10)
½z; w ¼ 1 (C1) j;k;m¼0
but differ from the canonical creation and destruction operators The generalized coherent states are defined via
could lead to interesting representations that would be useful for X X 1=2
representing classical dynamics. The study of such generalized rep- jzi ¼ ezw j0i ¼ z j jji=j! ¼ z j gjk j^ki=j!: (C11)
resentations is the domain of the theory of reproducing kernel j jk
Hilbert spaces (RKHS).
These states are eigenstates of the destruction operator z because
The Stone–von Neumann theorem166 states that if these opera-
they satisfy zjzi ¼ zjzi. The inner product of these states is the
tors are sufficiently well-behaved, so that they satisfy an exponential
reproducing kernel for this space
form of the CCR, known as the Weyl commutation relations
X
(WCR) Ky ðzÞ ¼ hyjzi ¼ ðy Þk z j gjk1 =j!k!: (C12)
jk
ezs ewt ¼ est ewt ezs ; (C2)
then there is a unitary transformation that relates w and z to the The coherent space coordinates offer yet another equivalent param-
usual operators q and p ¼ i@=@q, respectively. In fact, for the eterization of the Hilbert space. Assume the Hilbert space inner
form above, they are analogous to the destruction and creation product in coherent state space is given by the definition
operators, a ¼ ðq þ ipÞ=21=2 and a† ¼ ðq ipÞ=21=2 , respectively. ð
In practice, the cases of interest often do satisfy the WCR, as was h/jwi ¼ / ðzÞwðzÞgðz; z Þdzdz =X; (C13)
the case for all examples explored in Kowalski192 and in Ref. 32.
Assume there is a vacuum state j0i annihilated by the destruc- where gðw; z Þ ¼ g ðz; w Þ is a positive Hermitian weight function
tion operator zj0i ¼ 0, and assume there is no finite-dimensional and X ¼ pd is a normalization constant. In order for the coherent
subspace that is left invariant by the action of w. Then, the infinite state inner product to be consistent, the Hilbert space metric must
tower of states be given by the moments of the weight function
ð
j ji :¼ wj j0i (C3) gjk ¼ ðz Þj z k gðz; z Þdzdz =j!k!X: (C14)
for j 2 N spans the Hilbert space. Let us define the notation
Thus, the choice of the coherent space weight function, gðz; z Þ, is
h jj :¼ jji† . In order to complete the definition of the Hilbert space,
tied to the choice of metric on Hilbert space, g.
one must specify the metric for these states
For example, in the simple yet important special case where
h jjki ¼ gjk : (C4) the metric is diagonal, i.e., the states jji ¼ w j j0i are orthogonal,
then the states can be normalized by the definition
Due to the fact that the Hilbert space inner product must be
1=2 1=2
Hermitian, g† ¼ g, the dual basis elements, hjj, which are defined to j^ji ¼ gjj jji ¼ gjj wj j0i (C15)
satisfy hjjki ¼ djk , are given by
X and one finds the representation
hjj :¼ h jjg1 ¼ hkjgjk1 : (C5)
k X
1
w¼ j^jihj ^ 1jðgjj =gj1j1 Þ1=2 ; (C16)
Then, using the CCR, one finds the representation j¼1
X
1 X
1
w¼ jjihj 1j; (C6) z¼ jj ^ 1ih^jjjðgjj =gj1j1 Þ1=2 : (C17)
j¼0 j¼1
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-34

C Author(s) 2023
V
For any quantity wj labeled by j, define series of orthogonal functions over a given weight function such as
j Fourier series, Hermite polynomials, Bessel functions, etc.
Y
wj ! :¼ wk : (C18) The interesting idea proposed by Kowalski192 is to encode the
k¼1 dynamics through translations on the space of coherent states
rather than through translations on the space of physical observ-
Then, the number states are related via
ables. While the action may be quite complicated in observable
jji :¼ ðgjj =gj1j1 Þ1=2 !j^ji: (C19) space, the Stone–von Neumann theorem guarantees that there is a
unitary transformation to coherent state space where the dynamics
For the standard Segal–Bargmann coherent states over the
is simply encoded in the usual manner: through trajectories in
complex plane, the weight function is gðjjwjjÞ ¼ exp ðjjwjj2 Þ coherent space.
which leads to gjk ¼ djk =j!. This implies that wj;j1 ¼ zj1;j ¼ j1=2 .
The number states are normalized jji ¼ j^ji, so that hzjji ¼ hzj^ji REFERENCES
¼ ðz Þj =ðj!Þ1=2 and the reproducing kernel is 1
M. Raymer and C. Monroe, “The U.S. National Quantum Initiative,”
Quantum Sci. Technol. 4, 020504 (2019).
Ky ðzÞ ¼ hyjzi ¼ ey z : (C20) 2
M. Tse et al., “Quantum-enhanced advanced LIGO detectors in the era of
For the standard Bergman space over the complex unit disk, gravitational-wave astronomy,” Phys. Rev. Lett. 123, 231107 (2019).
3
the weight function is 1 inside the unit disk and vanishes outside F. A. Narducci, A. T. Black, and J. H. Burke, “Advances toward fieldable atom
interferometers,” Adv. Phys.: X 7, 1946426 (2022).
the disk. This leads to qjk ¼ djk =ðj þ 1Þ!j! which yields wjþ1;j 4
A. Nordrum, “China demonstrates quantum encryption by hosting a video
¼ ðjðj þ 1ÞÞ1=2 and zj1;j ¼ ðj þ 1Þ1=2 . The number states are hzj^ji call,” in Proceedings of the IEEE Spectrum (2017).
5
A. Y. Kitaev, K. H. Shen, and M. N. Vyalyi, Classical and Quantum
¼ ðj þ 1Þ1=2 ðz Þj and the reproducing kernel is Computation (American Mathematical Society: Graduate Studies in
Mathematics, Providence, RI, 2002), Vol. 47.
Ky ðzÞ ¼ hyjzi ¼ 1=ð1 y zÞ2 : (C21) 6
M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum
Information (Cambridge University Press, New York, NY, 2010).
For Sz€ego kernels, the weight function has all of its support on 7
P. Kaye, R. LaFlamme, and M. Mosca, An Introduction to Quantum
the boundary of a domain in the complex plane, e.g., in the stan- Computing (Oxford University Press, New York, NY, 2010).
8
dard case, by the unit disk E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P.
ð Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, “Real-time dynamics
of lattice gauge theories with a few-qubit quantum computer,” Nature 534,
gjk ¼ ðz Þj z k gðz; z Þdðjjzjj 1Þdzdz =j!k!X: (C22)
516 (2016).
9
P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R.
For the standard Hardy space, the weight function is 1, which yields Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth et al., “Spectroscopic
qjk ¼ djk =j!2 . Hence, this leads to wj;j1 ¼ j and zj1;j ¼ 1. The un- signatures of localization with interacting photons in superconducting
qubits,” Science 358, 1175 (2017).
normalized number states are hzj^ji ¼ ðz Þj and the reproducing 10
F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas,
kernel is S. Boixo, F. G. S. L. Brandao, D. A. Buell et al., “Quantum supremacy using a
programmable superconducting processor,” Nature 574, 505 (2019).
Ky ðzÞ ¼ hyjzi ¼ 1=ð1 y zÞ: (C23) 11
Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H.
Deng, Y. Du, D. Fan et al., “Strong quantum computational advantage using a
This is the basis used by Carleman linearization. superconducting quantum processor,” Phys. Rev. Lett. 127, 180501 (2021).
The different choices of the Hilbert space metric are also tied 12
H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, L.-C. Peng, Y.-H.
to the domain over which the function spaces are defined. For Luo, D. Wu, S.-Q. Gong, H. Su et al., “Phase-programmable Gaussian boson
Segal–Bargmann space, the Hilbert space is composed of holomor- sampling using stimulated squeezed light,” Phys. Rev. Lett. 127, 180502
phic functions of the complex plane that vanish sufficiently rap- (2021).
13
J. Preskill, “Quantum computing in the NISQ era and beyond,” Quantum 2,
idly at complex infinity. In contrast, for the Bergman and Hardy
79 (2018).
spaces, the Hilbert space is composed of functions that are holo- 14
A. Abedi, N. Maitra, and E. K. Gross, “Exact factorization of the time-
morphic within the unit disk. For the standard coherent states, the dependent electron-nuclear wave function,” Phys. Rev. Lett. 105, 123002
eigenfunctions are localized with a Gaussian form factor near the (2010).
15
given point in phase space. In contrast, for the Bergman and A. Abedi, N. Maitra, and E. K. Gross, “Correlated electron-nuclear dynamics:
Hardy spaces, the eigenfunctions inside the unit disk are generated Exact factorization of the molecular wavefunction,” J. Chem. Phys. 137,
22A530 (2012).
by a singularity outside the unit disk. In these case, the kernel has 16
T. Schenkel, W. Dorland, A. Baczewski, M. Boshier, G. Collins, J. DuBois, A.
the form Ky ðzÞ ¼ ð1 y zÞa , and, so, the singularity actually Houck, T. Humble, N. Loureiro, C. Monroe, S. Parker et al., Fusion Energy
resides at the point z ¼ 1=y . Thus, the kernel becomes Sciences Roundtable on Quantum Information Science, 1–2 May 2018
unbounded if the coordinate y ever crosses over to the inside of (USDOE Office of Science, 2018).
17
the unit disk. I. Y. Dodin and E. A. Startsev, “On applications of quantum computing to
Another important point to keep in mind is that the observ- plasma simulations,” Phys. Plasmas 28, 092101 (2021).
18
D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, “Quantum algorithm
able state space, i.e., the corresponding number states, have not yet
for linear differential equations with exponentially improved dependence on
been specified. While the observable states could correspond to precision,” Commun. Math. Phys. 356, 1057 (2017).
complex analytic functions, they could also be defined in many 19
P. C. S. Costa, S. Jordan, and A. Ostrander, “Quantum algorithm for simulat-
other ways. For example, they could correspond to any infinite ing the wave equation,” Phys. Rev. A 99, 012323 (2019).
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-35

C Author(s) 2023
V
20 48
A. M. Childs, J.-P. Liu, and A. Ostrander, “High-precision quantum algo- P. Benioff, “The computer as a physical system: A microscopic quantum
rithms for partial differential equations,” Quantum 5, 574 (2021). mechanical Hamiltonian model of computers as represented by Turing
21
D. W. Berry and L. Novo, “Corrected quantum walk for optimal Hamiltonian machines,” J. Stat. Phys. 21, 563 (1980).
49
simulation,” Quantum Inf. Comput. 16, 1295 (2016). Y. I. Manin, Computable and Uncomputable (Sovetskoye Radio, Moscow,
22
G. H. Low and I. L. Chuang, “Optimal Hamiltonian simulation by quantum 1980) (in Russian).
50
signal processing,” Phys. Rev. Lett. 118, 010501 (2017). R. P. Feynman, “Simulating physics with computers,” Int. J. Theor. Phys. 21,
23
A. Engel, G. Smith, and S. E. Parker, “Quantum algorithm for the Vlasov 467 (1982).
51
equation,” Phys. Rev. A 100, 062315 (2019). R. P. Feynman, “Quantum mechanical computers,” Found. Phys. 16, 507
24 (1986).
I. Novikau, E. A. Startsev, and I. Y. Dodin, “Quantum signal processing for
52
simulating cold plasma waves,” Phys. Rev. A 105, 062444 (2022). J. Preskill, see http://theory.caltech.edu/preskill/ph229/ for “Lecture notes on
25 quantum computation, 2018.”
J. L. Park, “The concept of transition in quantum mechanics,” Found. Phys.
53
1(1), 23–33 (1970). H. Goldstein, C. P. Poole, and J. L. Safko, Classical Dynamics (Addison-
26 Wesley, New York, 2001).
D. Dieks, “Communication by EPR devices,” Phys. Lett. A 92, 271 (1982).
27 54
W. K. Wooters and W. H. Zurek, “A single quantum cannot be cloned,” E. C. G. Sudarshan and N. Mukunda, Classical Dynamics: A Modern
Nature 299, 802 (1982). Perspective (Addison-Wesley, New York, 1974).
28 55
S. K. Leyton and T. J. Osborne, “A quantum algorithm to solve nonlinear dif- H. D. Zeh, “On the interpretation of measurement in quantum theory,”
ferential equations,” arXiv:0812.4423 (2008). Found. Phys. 1, 69 (1970).
56
29
I. Joseph, “Koopman-von Neumann approach to quantum simulation of non- W. H. Zurek, “Decoherence, einselection, and the quantum origins of the
linear classical dynamics,” Phys. Rev. Res. 2, 043102 (2020). classical,” Rev. Mod. Phys. 75, 715 (2003).
57
30
I. Y. Dodin and E. A. Startsev, “Quantum computation of nonlinear maps,” W. H. Zurek, “Decoherence and the transition from quantum to classical—
arXiv:2105.07317 (2021). Revisited,” arXiv:quant-ph/0306072 (2003).
58
31
J.-P. Liu, H. O. Kolden, H. K. Krovi, N. F. Louriero, K. Trivisa, and A. M. W. F. Stinespring, “Positive functions on C -algebras,” Proc. Am. Math. Soc.
Childs, “Efficient quantum algorithm for dissipative nonlinear differential 6, 211 (1955).
59
equations,” Proc. Natl. Acad. Sci. 118, e2026805118 (2021). M.-D. Choi, “Completely positive linear maps on complex matrices,” Linear
32
A. Engel, G. Smith, and S. E. Parker, “Linear embedding of nonlinear dynami- Algebra Appl. 10, 285–290 (1975).
60
cal systems and prospects for efficient quantum algorithms,” Phys. Plasmas K. Krauss, States, Effects and Operations: Fundamental Notions of Quantum
28, 062305 (2021). Theory, Lecture Notes in Physics (Springer, New York, 1983).
61
33
Y. T. Lin, R. B. Lowrie, D. Aslangil, Y. Subaşi, and A. T. Sornborger, G. Lindblad, “On the generators of quantum dynamical semigroups,”
“Koopman-von Neumann mechanics and the Koopman representation: A Commun. Math. Phys. 48, 119 (1976).
62
V. Gorini, A. Kossakowski, and E. Sudarshan, “Completely positive dynamical
perspective on solving nonlinear dynamical systems with quantum com-
semigroups of n-level systems,” J. Math. Phys. 17, 821 (1976).
puters,” arXiv:2202.02188 [quant-ph] (2022). 63
34 D. Manzano, “A short introduction to the Lindblad master equation,” AIP
H. Krovi, “Improved quantum algorithms for linear and nonlinear differential
Adv. 10, 025106 (2020).
equations,” arXiv:2202.01054 (2022). 64
35 S. Heinrich and E. Novak, “Optimal summation and integration by determin-
F. Gaitan, “Finding solutions of the Navier–Stokes equations through quan-
istic, randomized, and quantum algorithms,” in Monte Carlo and Quasi-
tum computing: Recent progress, a generalization, and next steps forward,”
Monte Carlo Methods 2000 (Springer, 2002), pp. 50–62.
Adv. Quantum Technol. 4, 2100055 (2021). 65
36 X. Wu, S. L. Tomarken, N. A. Petersson, L. A. Martinez, Y. J. Rosen, and J. L.
S. Jin and N. Liu, “Quantum algorithms for computing observables of nonlin-
DuBois, “High-fidelity software-defined quantum logic on a superconducting
ear partial differential equations,” arXiv:2202.07834 (2022).
37 qudit,” Phys. Rev. Lett. 125, 170502 (2020).
A. Montanaro, “Quantum speedup of Monte Carlo methods,” Proc. R. Soc. A 66
A. M. Childs, “Universal computation by quantum walk,” Phys. Rev. Lett.
471, 20150301 (2015).
38 102, 180501 (2009).
B. Kacewicz, “Almost optimal solution of initial-value problems by random- 67
A. M. Childs, D. Gosset, and Z. Webb, “Universal computation by multiparti-
ized and quantum algorithms,” J. Complexity 22, 676 (2006). cle quantum walk,” Science 339, 791 (2013).
39
D. An, N. Linden, J.-P. Liu, A. Montanaro, C. Shao, and J. Wang, “Quantum- 68
E. Farhi, J. Goldstone, S. Gutman, and M. Sipser, “Quantum computation by
accelerated multilevel Monte Carlo methods for stochastic differential equa- adiabatic evolution,” arXiv:quant-ph/0001106 (2000).
tions in mathematical finance,” Quantum 5, 481 (2021). 69
E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, and A. Lundgre, “A quantum
40
G. Xu, A. J. Daley, P. Givi, and R. D. Somma, “Turbulent mixing simulation adiabatic evolution algorithm applied to random instances of an NP-
via a quantum algorithm,” AIAA J. 56, 687–699 (2018). complete problem,” Science 292, 472 (2001).
41
G. Xu, A. J. Daley, P. Givi, and R. D. Somma, “Quantum algorithm for the 70
S. Morita and H. Nishimori, “Mathematical foundation of quantum
computation of the reactant conversion rate in homogeneous turbulence,” annealing,” J. Math. Phys. 49, 125210 (2008).
Combust. Theory Modell. 23, 1090–1104 (2019). 71
R. Raussendorf, D. E. Browne, and H. J. Briegel, “Measurement-based quan-
42
E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimiza- tum computation on cluster states,” Phys. Rev. A 68, 022312 (2003).
tion algorithm,” arXiv:1411.4028 (2014). 72
D. M. Greenberger, M. A. Horne, and A. Zeilinger, “Going beyond Bell’s theo-
43
A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and rem,” arXiv:0712.0921 (2007).
J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for 73
J. S. Bell, “On the Einstein Podolsky Rosen paradox,” Phys. Phys. Fiz. 1, 195
small molecules and quantum magnets,” Nature 549, 242–246 (2017). (1964).
44 74
M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, “Variational quan- C. H. Bennett, “The thermodynamics of computation—A review,” Int. J.
tum algorithms for nonlinear problems,” Phys. Rev. A 101, 010301(R) (2015). Theor. Phys. 21, 905 (1982).
45 75
R. Steijl, “Quantum algorithms for fluid simulations,” in Advances in M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative
Quantum Communication and Information (IntechOpen Rijeka, Licko- dynamics in quantum optics,” Rev. Mod. Phys. 70, 101 (1998).
Senjska, Croatia, 2020). 76
J. Kempe, “Quantum random walks: An introductory overview,” Contemp.
46
R. Steijl, “Quantum algorithms for nonlinear equations in fluid mechanics,” Phys. 44, 307–327 (2003).
77
in Quantum Computing and Communications (IntechOpen, London, 2020). A. M. Childs, E. Farhi, and S. Gutmann, “An example of the difference
47
K. P. Griffin, S. S. Jain, T. J. Flint, and W. H. R. Chan, “Investigations of quan- between quantum and classical random walks,” Quantum Inf. Process. 1,
tum algorithms for direct numerical simulation of the Navier–Stokes equa- 35–43 (2002).
78
tions,” in Proceedings of the Center for Turbulence Research Annual Research F. Gaitan, Quantum Error Correction and Fault Tolerant Quantum
Briefs (Center for Turbulence Research, 2019), pp. 347–363. Computing (CRC Press, 2008).
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-36

C Author(s) 2023
V
79 106
I. Djordjevic, Quantum Information Processing and Quantum Error S. Heinrich, “Quantum summation with an application to integration,”
Correction: An Engineering Approach (Academic Press, 2012). J. Complexity 18, 1–50 (2001).
80 107
D. A. Lidar and T. A. Brun, Quantum Error Correction (Cambridge G. Brassard, F. Dupuis, S. Gambs, and A. Tapp, “An optimal quantum algo-
University Press, 2013). rithm to approximate the mean and its application for approximating the
81
A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, “Surface median of a set of points over an arbitrary distance,” arXiv:1106.4267 (2011).
108
codes: Towards practical large-scale quantum computation,” Phys. Rev. A 86, S. Lloyd, “Universal quantum simulators,” Science 273, 1073 (1996).
109
032324 (2012). D. S. Abrams and S. Lloyd, “Quantum algorithm providing exponential speed
82 increase for finding eigenvalues and eigenvectors,” Phys. Rev. Lett. 83, 5162 (1999).
S. J. Devitt, W. J. Munro, and K. Nemoto, “Quantum error correction for
110
beginners,” Rep. Prog. Phys. 76, 076001 (2013). G. Benenti and G. Strini, “Quantum simulation of the single-particle
83 Schr€ odinger equation,” Am. J. Phys. 76, 657 (2008).
B. M. Terhal, “Quantum error correction for quantum memories,” Rev. Mod.
111
Phys. 87, 307 (2015). M. Suzuki, “General theory of fractal path integrals with applications to many-
84 body theories and statistical physics,” J. Math. Phys. 32, 400 (1991).
E. T. Campbell, B. M. Terhal, and C. Vuillot, “Roads towards fault-tolerant
112
universal quantum computation,” Nature 549, 172 (2017). A. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, “Theory of Trotter error
85 with commutator scaling,” Phys. Rev. X 11, 011020 (2021).
J. Roffe, “Quantum error correction: An introductory guide,” Contemp. Phys.
113
60, 226 (2019). D. Aharonov, A. Ambainis, J. Kempe, and U. Vazirani, “Quantum walks on
86
C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, “Trapped-ion graphs,” in Proceedings of the Thirty-Third Annual ACM Symposium on
quantum computing: Progress and challenges,” Appl. Phys. Rev. 6, 021314 Theory of Computing (ACM, 2001), pp. 50–59.
114
(2019). R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals
87
P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. (McGraw-Hill, New York, 1965), p. 35.
115
Oliver, “A quantum engineer’s guide to superconducting qubits,” Appl. Phys. D. Aharonov and A. Ta-Shma, “Adiabatic quantum state generation and statis-
Rev. 6, 021318 (2019). tical zero knowledge,” arXiv:quant-ph/0301023 (2003).
116
88
J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. S. Allman, A. M. Childs, “On the relationship between continuous- and discrete-time
C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. quantum walk,” Commun. Math. Phys. 294, 581 (2010).
117
Neyenhuis, “Demonstration of the trapped-ion quantum CCD computer D. W. Berry, A. M. Childs, and R. Kothari, “Hamiltonian simulation with
architecture,” Nature 592, 209 (2021). nearly optimal dependence on all parameters,” in Proceedings of the 56th
89
E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, “Topological quantum mem- Annual IEEE Symposium on Foundations of Computer Science (FOCS) (IEEE,
ory,” J. Math. Phys. 43, 4452 (2002). 2015), p. 792.
118
90
S. Bravyi and A. Kitaev, “Universal quantum computation with ideal Clifford A. M. Childs and N. Wiebe, “Hamiltonian simulation using linear combina-
tions of unitary operations,” Quantum Inf. Comput. 12, 901 (2012).
gates and noisy ancillas,” Phys. Rev. A 71, 022316 (2005). 119
91 D. W. Berry, R. Childs, A. M. Cleve, R. Kothari, and R. D. Somma,
R. Babbush, J. R. McClean, M. Newman, C. Gidney, S. Boixo, and H. Neven,
“Exponential improvement in precision for simulating sparse Hamiltonians,”
“Focus beyond quadratic speedups for error-corrected quantum advantage,”
in Proceedings of the 46th ACM Symposium on Theory of Computing (STOC
PRX Quantum 2, 010103 (2021).
92 2014) (ACM, 2014), p. 792.
S. Krinner, N. Lacroix, A. Remm, A. Di Paolo, E. Genois, C. Leroux, C. 120
D. W. Berry, A. M. Childs, R. Cleve, and R. Kothari, “Simulating Hamiltonian
Hellings, S. Lazar, F. Swiadek, J. Herrmann et al., “Realizing repeated quan-
dynamics with a truncated Taylor series,” Phys. Rev. Lett. 114, 090502 (2015).
tum error correction in a distance-three surface code,” Nature 605, 669–674 121
L. Novo and D. W. Berry, “Improved Hamiltonian simulation via a truncated
(2022).
93 Taylor series and corrections,” Quantum Inf. Comput. 17, 0623 (2016).
C. Ryan-Anderson, J. G. Bohnet, K. Lee, D. Gresh, A. Hankin, J. P. Gaebler, 122
J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, “A grand unification
D. Francois, A. Chernoguzov, D. Lucchetti, N. C. Brown, T. M. Gatterman, S. of quantum algorithms,” Phys. Rev. X Quantum 2, 040203 (2021).
K. Halit, K. Gilmore, J. A. Gerber, B. Neyenhuis, D. Hayes, and R. P. Stutz, 123
A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear sys-
“Realization of real-time fault-tolerant quantum error correction,” Phys. Rev. tems of equations,” Phys. Rev. Lett. 103, 150502 (2009).
X 11, 041058 (2021). 124
A. Ambainis, “Variable time amplitude amplification and a faster quantum
94
D. Deutsch, “Quantum theory, the Church-Turing principle and the universal algorithm for solving systems of linear equations,” arXiv:1010.4458 (2010).
quantum computer,” Proc. R. Soc. London, Ser. A 400, 97 (1985). 125
A. M. Childs, R. Kothari, and R. Somma, “Quantum algorithm for systems of
95
D. Deutsch and R. Josza, “Quantum complexity theory,” Proc. R. Soc. linear equations with exponentially improved dependence on precision,”
London, Ser. A 439, 553 (1992). SIAM J. Comput. 46, 1920 (2017).
96
E. Bernstein and U. Vazirani, “Quantum complexity theory,” SIAM J. 126
A. Gilyen, Y. Su, G. H. Low, and N. Wiebe, “Quantum singular value transfor-
Comput. 26, 1411 (1997). mation and beyond: Exponential improvements for quantum matrix arithme-
97
D. R. Simon, “On the power of quantum computation,” SIAM J. Comput. 26, tics,” in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory
1474 (1997). of Computing (STOC) (ACM, 2019), p. 193.
98
P. W. Shor, “Polynomial-time algorithms for prime factorization and discrete 127
J. M. Martyn, Y. Liu, Z. E. Chin, and I. L. Chuang, “Efficient fully-coherent
logarithms on a quantum computer,” SIAM J. Comput. 26, 1484 (1997). Hamiltonian simulation,” arXiv:2110.11327 (2021).
99
P. W. Shor, “Scheme for reducing decoherence in quantum computer memo- 128
D. W. Berry, “High-order quantum algorithms for solving linear differential
ry,” Phys. Rev. A 52, R2493 (1995). equations,” J. Phys. A: Math. Theor. 47, 105301 (2014).
100 129
L. K. Grover, “Quantum mechanics helps in searching for a needle in a hay- G. H. Low and I. L. Chuang, “Hamiltonian simulation by qubitization,”
stack,” Phys. Rev. Lett. 79, 4709 (1997). Quantum 3, 163 (2019).
101 130
L. K. Grover, “Quantum computers can search rapidly by using almost any I. Novikau, E. A. Startsev, and I. Y. Dodin, “Quantum signal processing for
transformation,” Phys. Rev. Lett. 80, 4329 (1998). simulating cold plasma waves,” arXiv:2112.06086 (2022).
102 131
L. K. Grover, “A framework for fast quantum-mechanical algorithms,” A. M. Childs and W. van Dam, “Quantum algorithms for algebraic problems,”
arXiv:quant-ph/9711043 (1998). Rev. Mod. Phys. 82, 1 (2010).
103 132
G. Brassard, P. Høyer, M. Mosca, and A. Tapp, “Quantum amplitude amplifi- A. Montanaro, “Quantum algorithms: An overview,” npj Quantum Inf. 2,
cation and estimation,” arXiv:quant-ph/0005055 (2000). 15023 (2010).
104 133
G. Brassard, M. Mosca, and A. Tapp, “Quantum amplitude amplification and M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R.
estimation,” Quantum Computation and Information (AMS Contemporary McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, “Variational quantum
Mathematics, 2002), Vol. 305, p. 53. algorithms,” Nat. Rev. Phys. 3, 625 (2021).
105 134
D. S. Abrams and C. P. Williams, “Fast quantum algorithms for numerical J. B. Parker and I. Joseph, “Quantum phase estimation for a class of general-
integrals and stochastic processes,” arXiv:quant-ph/9908083 (1999). ized eigenvalue problems,” Phys. Rev. A 102, 022422 (2020).
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-37

C Author(s) 2023
V
135 163
N. Shenvi, J. Kempe, and K. B. Whaley, “Quantum random-walk search algo- D. Giannakis, A. Ourmazd, P. Pfeffer, J. Schumacher, and J. Slawinska,
rithm,” Phys. Rev. A 67, 052307 (2003). “Embedding classical dynamics in a quantum computer,” Phys. Rev. A 105,
136
A. M. Childs and J. Goldstone, “Spatial search by quantum walk,” Phys. Rev. 052404 (2022).
164
A 70, 022314 (2004). N. Aronszajn, “Theory of reproducing kernels,” Trans. Am. Math. Soc. 68,
137
M. Szegedy, “Quantum speedup of Markov chain based algorithms,” in 337 (1950).
165
Proceedings of the 45th Annual IEEE Symposium on Foundations of Computer V. I. Paulsen and M. Raghupathi, An Introduction to the Theory of
Science (FOCS) (IEEE, 2004), p. 32. Reproducing Kernel Hilbert Spaces (Cambridge University Press, New York,
138
R. D. Somma, S. Boixo, H. Barnum, and E. Knill, “Quantum simulations of 2016).
166
classical annealing processes,” Phys. Rev. Lett. 101, 130504 (2008). B. Hall, Quantum Theory for Mathematicians, Graduate Texts in Mathematics
139
P. Wocjan and A. Abeyesinghe, “Speedup via quantum sampling,” Phys. Rev. (Springer, New York, 2013), Vol. 267.
167
A 78, 042336 (2008). M. H. Kalos and P. A. Whitlock, Monte Carlo Methods (Wiley-VCH Verlag
140
P. Wocjan, C.-F. Chiang, D. Nagaj, and A. Abeyesinghe, “Quantum algorithm GmbH and Co. KGaA, Weinheim, 2008).
168
for approximating partition functions,” Phys. Rev. A 80, 022340 (2009). R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles
141
A. Ambainis, “Quantum walk algorithm for element distinctness,” SIAM J. (Taylor and Francis, New York, 1998).
169
Comput. 37, 1700169 (2007). C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation,
142
D. W. Berry and A. M. Childs, “Black-box Hamiltonian simulation and uni- Institute of Physics: Series in Plasma Physics (Taylor and Francis, New York,
tary implementation,” Quantum Inf. Comput. 12, 29 (2012). 2004).
143 170
E. Hairer, C. Lubisch, and G. Wanner, Geometric Numerical Integration: G. Benenti, G. Casati, S. Montangero, and D. L. Shepelyansky, Quantum Inf.
Structure-Preserving Algorithms for Ordinary Differential Equations, Springer Process. 3, 273 (2004).
171
Series in Computational Mathematics (Springer, New York, 2008), Vol. 31. I. Mezić and A. Banaszuk, “Comparison of systems with complex behavior,”
144
G. Benenti, G. Casati, S. Montangero, and D. L. Shepelyansky, Phys. Rev. Lett. Physica D 192, 101 (2004).
172
87, 227901 (2001). I. Mezić, “Spectral properties of dynamical systems, model reduction and
145
Y. Shi, J. Xiao, H. Qin, and N. J. Fisch, “Simulations of relativistic quantum decompositions,” Nonlinear Dyn. 41, 309 (2005).
173
plasmas using real-time lattice scalar QED,” Phys. Rev. E 97, 053206 (2018). S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations
146
Y. Shi, H. Qin, and N. J. Fisch, “Plasma physics in strong-field regimes: from data by sparse identification of nonlinear dynamical systems,” Proc.
Theories and simulations,” Phys. Plasmas 28, 042104 (2021). Natl. Acad. Sci. 113, 3932 (2016).
147 174
D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, “Efficient quantum algo- E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear
rithms for simulating sparse Hamiltonians,” Commun. Math. Phys. 270, 359 dynamics for model predictive control in the low-data limit,” Proc. R. Soc. A
(2007). 474, 20180335 (2018).
148 175
A. M. Childs and J.-P. Liu, “Quantum spectral methods for differential equa- S. Jin and S. Osher, “A level set method for the computation of multi-valued
tions,” Commun. Math. Phys. 375, 1427 (2020). solutions to quasi-linear hyperbolic PDE’s and Hamilton-Jacobi equations,”
149
A. Montanaro and S. Pallister, “Quantum algorithms and the finite element Commun. Math. Sci. 1, 575–591 (2003).
176
method,” Phys. Rev. A 93, 032324 (2016). H. Liu, S. Osher, and R. Tsai, “Multi-valued solution and level set methods in
150
N. Guo, K. Mitarai, and K. Fujii, “Nonlinear transformation of complex ampli- computational high frequency wave propagation,” Commun. Comput. Phys.
tudes via quantum singular value transformation,” arXiv:2107.10764 (2022). 1, 765–804 (2006).
151 177
Z. Holmes, N. Coble, A. T. Sornborger, and Y. Subasi, “On nonlinear transfor- B. O. Koopman, “Hamiltonian systems and transformations in Hilbert space,”
mations in quantum computation,” arXiv:0812.4423 (2022). Proc. Natl. Acad. Sci. 17, 315 (1931).
152 178
I. Y. Dodin and E. A. Startsev, “On applications of quantum computing to B. O. Koopman and J. von Neumann, “Dynamical systems of continuous
plasma simulations,” arXiv:2105.07317 (2020). spectra,” Ann. Math. 18, 255 (1932).
153 179
S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, J. von Neumann, “Zur operatorenmethode in der klassischen mechanik,” Ann.
and T. Palmer, “Quantum algorithm for nonlinear differential equations,” Math. 33, 587 (1932).
180
arXiv:2011.06571 (2020). J. von Neumann, “Zusatze zur arbeit, zur operatorenmethode…,” Ann. Math.
154
B. M. Boghosian and W. E. Taylor, “Quantum lattice gas models for the 33, 789 (1932).
181
many-body Schr€odinger equation,” Int. J. Mod. Phys. C 8(4), 705 (1997). L. van Hove, On Certain Unitary Representations of an Infinite Group of
155
J. Yepez, B. M. Boghosian, and W. E. Taylor, “An efficient and accurate quan- Transformations (World Scientific, NJ, 1951).
182
tum lattice-gas model for the many-body Schr€ odinger wave equation,” D. I. Bondar, F. Gay-Balmaz, and C. Tronci, “Koopman wavefunctions and
Comput. Phys. Commun. 146, 280 (2002). classical-quantum correlation dynamics,” Proc. R. Soc. A 475, 20180879
156
F. Gaitan, “Finding flows of a Navier–Stokes fluid through quantum (2019).
183
computing,” npj Quantum Inf. 6, 61 (2020). C. Tronci and I. Joseph, “Koopman wavefunctions and Clebsch variables in
157
B. Z. Kacewicz, “Optimal solution of ordinary differential equations,” Vlasov-Maxwell kinetic theory,” J. Plasmas Phys. 87, 835870402 (2021).
184
J. Complexity 3, 451 (1987). J. Welch, D. Greenbaum, S. Mostame, and A. Aspuru-Guzik, “Efficient quan-
158
N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. tum circuits for diagonal unitaries without ancillas,” New J. Phys. 16, 033040
Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn et al., “Quantum optimization (2014).
185
using variational algorithms on near-term quantum devices,” Quantum Sci. S. Das and D. Giannakis, “On harmonic Hilbert spaces on compact Abelian
Technol. 3, 030503 (2018). groups,” arXiv:1912.11664 (2022).
159 186
J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. C. Xue, Y.-C. Wu, and G.-P. Guo, “Quantum homotopy perturbation method
Rungger, G. H. Booth, and J. Tennyson, “The variational quantum eigensolver: for nonlinear dissipative ordinary differential equations,” New J. Phys. 23,
A review of methods and best practices,” Phys. Rep. 986, 1–128 (2021). 123035 (2021).
160 187
K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, D. An, D. Fang, S. Jordan, J.-P. Liu, G. H. Low, and J. Wang, “Efficient quan-
M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.- tum algorithm for nonlinear reaction-diffusion equations and energy
C. Kwek, and A. Aspuru-Guzik, “Noisy intermediate-scale quantum algo- estimation,” arXiv:2205.01141 (2022).
188
rithms,” Rev. Mod. Phys. 94, 015004 (2022). T. Carleman, “Application de la theorie des equationa integrales lineaires aux
161
K. Kubo, Y. O. Nakagawa, S. Endo, and S. Nagayama, “Variational quantum systems d’equations differentielles non lineaires,” Acta Math. 59, 63 (1932).
189
simulations of stochastic differential equations,” Phys. Rev. A 103, 052425 W.-H. Steeb, “Embedding of nonlinear finite dimensional systems in linear
(2021). infinite dimensional systems and Bose operators,” Hadronic J. 6, 68 (1987).
162 190
L. Bittel and M. Kliesch, “Training variational quantum algorithms is NP- K. Kowalski, “Hilbert space description of classical dynamical systems I,”
hard,” Phys. Rev. Lett. 127, 120502 (2021). Physica A 145, 408 (1987).
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-38

C Author(s) 2023
V
191 205
K. Kowalski and W.-H. Steeb, Nonlinear Dynamical Systems and Carleman P. Jaquod, P. G. Silvestrov, and C. W. J. Beenakker, “Golden rule decay versus
Linearization (World Scientific, Singapore, 1991). Lyapunov decay of the quantum Loschmidt echo,” Phys. Rev. E 64,
192
K. Kowalski, Methods of Hilbert Spaces in the Theory of Nonlinear Dynamical 055203(R) (2001).
206
Systems (World Scientific, Singapore, 1991). G. Benenti and G. Casati, “Quantum-classical correspondence in perturbed
193
A. Pizzamiglio, S. Y. Chang, M. Bondani, S. Montangero, D. Gerace, and G. chaotic systems,” Phys. Rev. E 65, 066205 (2004).
207
Benenti, “Dynamical localization simulated on actual quantum hardware,” B. Georgeot and D. L. Shepelyansky, “Exponential gain in quantum comput-
Entropy 23, 654 (2021). ing of quantum chaos and localization,” Phys. Rev. Lett. 86, 2890 (2001).
194 208
M. D. Porter and I. Joseph, “Observability of fidelity decay at the Lyapunov V. E. Zakharov, V. S. L’vov, and G. Falkovich, Kolmogorov Spectra of
rate in few-qubit quantum simulations,” arXiv:2110.07767 (2022). Turbulence. I. Wave turbulence, Springer Series in Nonlinear Dynamics
195
M. D. Porter and I. Joseph, “Impact of dynamics, entanglement, and (Springer-Verlag, New York, 1992).
209
Markovian noise on the fidelity of few-qubit digital quantum simulation,” S. Nazarenko, Wave Turbulence, Lecture Notes in Physics (Springer, New
arXiv:2206.04829 (2022). York, 2011).
196 210
Y. Shi, A. R. Castelli, X. Wu, I. Joseph, V. Geyko, F. R. Graziani, S. B. Libby, J. B. Y. Shi, H. Qin, and N. J. Fisch, “Three-wave scattering in magnetized plas-
Parker, Y. J. Rosen, L. A. Martinez, and J. L. DuBois, “Simulating nonnative cubic mas: From cold fluid to quantized Lagrangian,” Phys. Rev. E 96, 023204
interactions on noisy quantum machines,” Phys. Rev. A 103, 062608 (2021). (2017).
197 211
J. Zylberman, G. Di Molfetta, M. Brachet, N. F. Loureiro, and F. Debbasch, P. J. Karalekas, N. A. Tezak, E. C. Peterson, C. A. Ryan, M. P. da Silva, and R.
“Hybrid quantum-classical algorithm for hydrodynamics,” arXiv:2202.00918 S. Smith, “A quantum-classical cloud platform optimized for variational hybrid
(2022). algorithms,” Quantum Sci. Technol. 5, 024003 (2020).
198 212
P. Jaquod and C. Petitjean, “Decoherence, entanglement and irreversibility in R. S. Smith, M. J. Curtis, and W. J. Zeng, “A practical quantum instruction set
quantum dynamical systems with few degrees of freedom,” Adv. Phys. 58, 67 architecture,” arXiv:1608.03355 [quant-ph] (2016).
213
(2009). N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbr€ uggen, and S. J. Glaser,
199
G. Casati, B. V. Chirikov, F. M. Izrailev, and J. Ford, Stochastic Behavior in “Optimal control of coupled spin dynamics: Design of NMR pulse
Classical and Quantum Hamiltonian Systems, Lecture Notes in Physics, edited sequences by gradient ascent algorithms,” J. Magn. Reson. 172, 296–305
by G. Casati and J. Ford (Springer, New York, 1979), Vol. 93. (2005).
200 214
B. V. Chirikov, F. M. Izrailev, and D. L. Shepelyansky, “Dynamical sto- S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. K€ ockenberger, R.
chasticity in classical and quantum mechanics,” Sov. Sci. Rev. 2C, 209 Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbr€ uggen et al.,
(1981). “Training Schr€ odinger’s cat: Quantum optimal control,” Eur. Phys. J. D 69,
201
B. Chirikov, F. Izrailev, and D. Shepelyansky, “Quantum chaos: Localization 279 (2015).
215
vs. ergodicity,” Physica D 33, 77–88 (1988). N. A. Petersson, F. M. Garcia, A. E. Copeland, Y. L. Rydin, and J. L. DuBois,
202
S. Fishman, D. R. Grempel, and R. E. Prange, “Chaos, quantum recurrences, “Discrete adjoints for accurate numerical optimization with application to
and Anderson localization,” Phys. Rev. Lett. 49, 509 (1982). quantum control,” arXiv:2001.01013 (2020).
203 216
A. Peres, “Stability of quantum motion in chaotic and regular systems,” Phys. J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP: An open-source Python
Rev. A 30, 1610 (1984). framework for the dynamics of open quantum systems,” Comput. Phys.
204
R. A. Jalabert and H. M. Pastawski, “Environment-independent decoherence Commun. 183, 1760–1772 (2012).
217
rate in classically chaotic systems,” Phys. Rev. Lett. 86, 2490 (2001). J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP 2: A Python framework for
the dynamics of open quantum systems,” Comput. Phys. Commun. 184,
1234–1240 (2013).
Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-39

C Author(s) 2023
V

Quantum Computing For Fusion Energy Science Applications

Uploaded by

Copyright:

Available Formats

Quantum Computing For Fusion Energy Science Applications

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Quantum Computing For Fusion Energy Science Applications

Uploaded by

Copyright:

Available Formats

Quantum computing for fusion energy

I. Joseph, Y. Shi, M. D. Porter, et al.

This paper was selected as Featured

ARTICLES YOU MAY BE INTERESTED IN

Editorial: Preface to the 30th volume of Physics of Plasmas

Phys. Plasmas 30, 010501 (2023); https://doi.org/10.1063/5.0123765 30, 010501

Quantum computing for fusion energy science

I. Joseph,a) Y. Shi, M. D. Porter, A. R. Castelli, V. I. Geyko, F. R. Graziani, S. B. Libby,

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-1

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-2

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-3

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-4

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-5

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-6

Given the deﬁnition above, when measurements are made along

Rn^ ðhÞ ¼ eih^n ~r =2 ¼ cos ðh=2Þ1 i sin ðh=2Þ^

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-7

FIG. 2. The Hilbert space of multiple

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-8

required to generate the whole group. Another important example is

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-9

• Directly sampling a quantum state to estimate the probability

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-10

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-11

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-12

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-13

FIG. 5. Circuit diagram for the quantum

FIG. 6. Circuit diagram for quantum phase

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-14

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-15

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-16

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-17

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-18

D. Simulation E. Estimation of physical observables

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-19

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-20

FIG. 9. A sufﬁciently large number of

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-21

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-22

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-23

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-24

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-25

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-26

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-27

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-28

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-29

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-30

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-31

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-32

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-33

Then, the complex analytic KvN equation is generalized to X

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-34

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-35

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-36

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-37

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-38

Phys. Plasmas 30, 010501 (2023); doi: 10.1063/5.0123765 30, 010501-39

You might also like

Rn^ ðhÞ ¼ eih^n ~r =2 ¼ cos ðh=2Þ1 i sin ðh=2Þ^