Quantum computing for fusion energy science applications

The three pillars of quantum information science (QIS)—quantum sensing, quantum communications, and quantum computing (QC)—promise to have transformative impact on science, engineering, and technology as we know it.¹ This article presents a pedagogical introduction to quantum computing and reviews recent research to develop and apply quantum algorithms and utilize quantum computing hardware platforms for fusion energy science (FES) applications. The interesting results obtained so far make it hopeful that QIS may one day lead to game-changing capabilities for FES.

Before diving into quantum computing (QC), let us briefly mention the other two pillars of QIS. First, quantum sensing techniques are already being used today to improve measurement sensitivity by orders of magnitude. Instead of being limited by the central limit theorem to yield a noise-to-signal ratio of $1/S$⁠, where S is the number of samples, using intrinsically quantum entangled states, such as squeezed states, as sensitive probes can reduce the detection threshold to the Heisenberg limit $1/S$⁠. Such techniques have been used to improve the sensitivity of the Laser Interferometer Gravitational-wave Observatory (LIGO) gravitational wave detector by a factor of 2 for nearly a decade.² New detectors based on nitrogen vacancy centers (NV centers) in diamond have provided unprecedented sensitivity for measurements of magnetic and electric fields as well as temperature and pressure. Advances in atom interferometry³ have led to a revolution in gravitational and inertial (gyroscopic) sensing, now transitioning to real world applications, such as navigation without the Global Positioning System (GPS), rapid passive sensing of mass distributions, and underground structure discovery, as well as basic science applications, such as gravitational wave detection in the frequency regime between LIGO and the Laser Interferometer Space Antenna (LISA).

Second, quantum communications offer the possibility of secure and intrinsically parallel information transfer.⁴ The goal is to transform and transport quantum information over long distances and to transduce information between different quantum hardware platforms. This field generated some of the earliest rigorous proofs that combining entangled states with quantum communication protocols could surpass their classical counterparts in their ability to carry and manipulate information. The application to network security is considered to be so important that many researchers are already working on the quantum internet.

Finally, quantum computing, which is the focus of this review article, promises to efficiently perform quantum algorithms that have the power to achieve polynomial to exponential improvements in complexity while manipulating exponentially large quantum memory resources.^5–7 The recent growth in the capabilities of today's quantum computing devices has spurred great interest in quantum simulation, and they have been used to perform key demonstrations of quantum calculations.^8–10 The potential impact of this so-called quantum advantage is so exciting that many government and private industry laboratories around the world are working to perfect the technology. In fact, in recent years, there have been claims that some demonstrations have already passed the point of quantum supremacy,^10–12 the point at which quantum computers can surpass even the world's best supercomputers at performing certain computational tasks.

While there has certainly been enormous progress over the last 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 correction. Because present-day and near-term quantum devices offer large quantum memory registers but lack error correction, the present era has been dubbed the noisy intermediate-scale quantum (NISQ) era.¹³ The search for the ultimate physical platform for realizing a quantum computer spurs research onward to invent new quantum materials, new hardware platforms, and new algorithms that can surpass today's limitations.

Fusion energy science (FES) is defined by the mission of Achieving a safe sustainable fusion energy source for the foreseeable future. Hence, it encompasses all fields of science required to achieve this goal, notably including plasma physics, nuclear physics, materials science, and chemistry. Thus, the FES mission requires developing accurate predictive models within all of these disparate fields, to form combined whole-device models that ultimately span a wide range of physical regimes from classical to quantum.

To stress the difficulty of this endeavor, FES requires calculations of strongly correlated quantum materials such as superconductors at temperatures below 1 meV to the collective motions of plasma at 1–10 keV as well as nuclear reactions at 10s of MeV. Whereas magnetic confinement fusion uses low-density plasmas on the order of 10²⁰ m⁻³, inertial confinement fusion seeks to obtain a compression of particle density on the order of 1000 times that of solid matter, for particle densities beyond $1030 m−3$⁠. Thus, this represents a range of 9 orders of magnitude in energy and 10 orders of magnitude in density (without considering the density of nuclear matter).

Note, however, that, today, the vast majority of scientific computations in the fusion energy sciences are classical in nature. The meaning of this statement is that such calculations seek to predict the time-dependent evolution of nonlinear sets of both ordinary differential equations (ODEs) and partial differential equations (PDEs). Ensembles of systems are typically evolved using classical probability distribution functions (PDFs), which are simply advected along the trajectories. While this is certainly true in the traditionally classical realm of plasma physics, it is also quite true in the traditional quantum realm of chemistry and materials science. The computational workhorses of density functional theory (DFT) and molecular dynamics crucially rely on the ability to efficiently evaluate nonlinear functions and to evolve nonlinear differential equations. The same is true for evolving classical kinetic equations that describe the evolution of nuclear and chemical species with reaction rate theory.

An essential limitation of the standard quantum computing paradigm is that, aside from measurements, only linear unitary operations can be performed efficiently in the quantum Hilbert space. If one believes that quantum mechanics (QM) is the correct physical theory of nature, then linear dynamics is all that is needed to describe the evolution of the universe. (A notable exception to this rule is that some physicists believe that the measurement process might require nonlinear dynamics to be described completely.) However, according to the exact factorization theorem^14,15 the process of generating reduced order models typically generates nonlinear interactions between subsystems. This is the reason that essentially all useful physical models that achieve significant complexity reduction, e.g., fluid dynamics, density function theory (DFT), and kinetic theory, are nonlinear.

Many scientific computations of interest to FES can be solved more efficiently using quantum algorithms (see the excellent overviews given in Refs. 16 and 17). Quantum algorithms promise to accelerate many of the algorithms that are central to scientific computing, including the Fourier transform, sparse linear solvers, and sparse Hamiltonian simulation. Quantum algorithms also exist for linear ODEs¹⁸ and linear PDEs^19,20 that can be cast in the form of initial value problems (IVPs), boundary value problems (BVPs), and eigenvalue problems (EVPs).

For example, for plasma and materials science applications, some IVPs are naturally posed as being Hamiltonian and, hence, can naturally be presented as a unitary evolution resulting from a Hermitian Hamiltonian. If the Hamiltonians of interest have a sparse structure, then they can be solved using efficient quantum algorithms for Hamiltonian simulation,^21,22 which we will discuss in detail in Sec. IV. Notably, Ref. 23 developed a quantum algorithm to solve the linearized Vlasov–Poisson system and Ref. 24 developed a quantum algorithm to solve the problem of cold plasma wave propagation and linear mode conversion. Yet, the majority of IVPs for reduced models include dissipation. Since the dynamics is not Hamiltonian, they need to be treated with specialized methods, for example, the general purpose quantum linear differential equation solver,¹⁸ discussed in Sec. III E 7.

For many fields of science and engineering, including biology, chemistry, and physics, a large share of today's computational resources are used for the simulation of classical nonlinear dynamics. Hence, it is important to understand whether quantum computers can provide similar gains in efficiency for the simulation of nonlinear problems. Nonlinear operations can be performed by forming the tensor product of identical states. However, due to the no-cloning theorem,^25–27 it is not possible to make copies of a quantum state; rather, one must prepare multiple identical states from scratch. Thus, if a nonlinear operation is to be iterated, e.g., in time, then this will require preparing a number of identical states that is exponentially large in the number of iteration steps.²⁸ 

Efficient quantum algorithms for simulating nonlinear physical processes have only recently begun to be invented and improved^17,29–34 and some of these efforts were initiated for FES applications.¹⁶ The new algorithms all use the strategy of embedding the nonlinear system within an infinite-dimensional linear system that is then truncated to finite dimension. For systems of nonlinear ODEs, the general method of using a quantum computer to solve for the evolution of the wavefunction corresponding to a classical PDF was first clearly described in Refs. 17 and 29, using what is known as the Koopman–von Neumann embedding. Using an alternative approach, Ref. 31 developed concrete quantum algorithms for dissipative differential equations with quadratic nonlinearity based on an approach known as Carleman embedding, which is applied to simulate the Burgers equation. Connections between these embedding techniques are discussed in Ref. 32, and comparisons of these approaches are elucidated using numerical examples in Ref. 33. Building upon ODE solvers, quantum algorithms for solving PDEs, including the Navier–Stokes equations, have also been explored.^35,36

Apart from deterministic methods, algorithms based on random processes have also been developed. For example, quantum versions of Monte Carlo algorithms can generally be performed with a quadratic speedup over their classical counterparts.³⁷ Similar concepts can then be used to speed up calculations of stochastic differential equations (SDEs).³⁸ General multi-level Monte Carlo methods for stochastic differential equations can also been accelerated.³⁹ In fact, quantum algorithms based on these concepts have been proposed for computing turbulent mixing and the reaction rate within turbulent flows.^40,41

The algorithms mentioned above requires fault-tolerant quantum computers to implement various oracles that are required by quantum subroutines. Alternatively, on NISQ devices, which have been shown to be capable of performing classical-quantum hybrid optimization problems,^42,43 advancing nonlinear equations may be treated as an optimization problem using variational algorithms.⁴⁴ In this case, evaluating the nonlinear functional only requires a fixed number of copies per time step. Hence, there have been recent proposals for simulating the Navier–Stokes equations using the variational technique.^45–47 

Even though simulating nonlinear problems is perhaps one of the most useful applications of scientific computation, relevant quantum algorithms are still in their infancy and there is much work to be done in terms of improving their speed and accuracy and verifying their performance on quantum computing hardware platforms. Once fully vetted and understood, they can be used as subroutines to provide a foundation for quantum programs that perform more useful and more complex tasks.

Section II will provide a brief introduction to the principles of quantum information and explain the reason that quantum memory registers can be considered to be much larger than their classical counterparts. Then, Sec. III will discuss the digital circuit model of universal quantum computation and review the key quantum subroutines that are useful for accelerating classical computations and for building quantum programs. Next, Sec. IV will discuss quantum approaches to simulating differential equations, the type of application that is used most by FES researchers today. Then, Sec. V will unify several recent proposals for using quantum algorithms to simulate nonlinear dynamics, including generalized approaches to Koopman–von Neumann (KvN) and Carleman embedding, as well as explain their differences. The penultimate section, Sec. VI will discuss recent findings in applying several leading quantum computing hardware platforms to quantum simulation problems of interest to FES.  Appendix B extends the KvN framework to the complex-analytic setting relevant to Carleman embedding.  Appendix C proves that different choices of the relevant Hilbert space bases, i.e., reproducing kernel Hilbert spaces (RKHS), depend on the choice of Hilbert space metric. The concluding section summarizes our findings and offers perspectives on the future outlook for quantum computing for FES applications.

Benioff,⁴⁸ Manin,⁴⁹ and Feynman,^50,51 were some of the first scientists to think seriously about the potential of using machines that satisfy the laws of quantum mechanics for performing scientific calculations. Feynman also left us with one of the great quotes in the field:⁵⁰ Nature is not classical, dammit, and if you want to make a simulation of nature, you'd better make it quantum mechanical, and by golly it's a wonderful problem, because it does not look so easy. Understanding this statement is the key to understanding why scientists believe that quantum computers can, in principle, be so powerful.

To see why quantum computers are fundamentally different from classical computers, let us begin by recalling the postulates of quantum mechanics^6,7,52 and compare them with postulates of classical mechanics (CM). A primer on Hilbert spaces is given in  Appendix A. In the following, we assume that the reader is familiar with Dirac's bra-ket notation.

In the following formulas, the Hilbert space will be treated as if it were discrete (which is typically true for digital models of quantum computation), but if it is continuous, then the sums should be replaced by integrals.

1. Physical states.

   • Pure physical states, denoted by the wavefunction ψ, are rays in Hilbert space.

     • Pure states are vectors in Hilbert space, generically, a linear superposition of basis elements that span the Hilbert space.

     • Pure states are elements of projective Hilbert space; i.e., the overall complex constant is unimportant, and the normalization can be set to unity, $||ψ||=1$⁠.

     • The Hilbert space is typically a space of functions that is associated with a configuration manifold, $m$⁠, but it can also be a discrete set.

   • Mixed physical states are probability distributions over pure states and, hence, are represented by a positive self-adjoint operator on Hilbert space, $ρ†=ρ$⁠, called the probability density matrix.

     • Due to Hermiticity, mixed states can be diagonalized

       $ρ=∑jpj|ψj⟩ ⟨ψj|.$

       (1)

       Because the eigenvalues, $pj≥0$⁠, must be non-negative this represents a probability distribution over a set of pure states, $|ψj⟩$⁠, with probability $pj$⁠.

     • Probability is normalized to unity,

       $Tr (ρ)=∑jpj=1.$

     • Proper mixed states must satisfy $Tr (ρa)<1$ for a > 1 and $Tr (ρa)>1$ for $0<a<1$⁠.

     • Mixed stated are required to represent the outcome of measurements.

   • The Hilbert space of a composite system is the tensor product of the Hilbert spaces of the subsystems.

     • The Hilbert space dimension of the composite system is the product of the dimensions of the subsystems.

     • Entangled states are states that cannot be represented as tensor products of the states of the individual subsystems; i.e., states that are a nontrivial superposition of the states of the subsystems.

2. Measurement of physical observables.

   • Physical observables are self-adjoint operators $O†=O$ on Hilbert space.

   • Measurements yield eigenvalues of the operator corresponding to the observable

     $O|ψα⟩=α|ψα⟩.$

     (2)
   • The probability of observing a particular eigenvalue, α, is given by the overlap between the eigenstate and the physical state

     $Prob(α)=⟨ψα|ρ|ψα⟩=∑jpj||⟨ψα|ψj⟩||2.$

     (3)
     Therefore, the expectation value of measuring an observable is given by

     $⟨O⟩=Tr (ρO)=∑jpj⟨ψj|O|ψj⟩.$

     (4)

   • Immediately after a measurement of eigenvalue α, the wavefunction collapses to the state $|ψα⟩$⁠.

3. Temporal evolution.

   • General open systems

     • Physical states evolve via linear superoperations that act on the density matrix

       $ρ(t)=e(t,t0)ρ(t0).$

       (5)

       The evolution must preserve the Hermiticity, positivity, and trace of any density matrix.

     • Continuous time evolution is generated by the generator of a linear superoperator

       $dρ/dt=l(t)ρ(t).$

       (6)

       The evolution must preserve the Hermiticity, positivity, and trace of any density matrix.

   • Ideal closed systems

     • Ideal closed systems evolve via a linear unitary, $U−1=U†$⁠, transformation of Hilbert space,

       $ψ(t)=Uψ(t0) ρ(t)=Uρ(t0)U†.$

       (7)
     • Continuous time evolution is generated by a Hermitian Hamiltonian operator, H

       $iℏ∂tψ=Hψ iℏ∂tρ=[H,ρ],$

       (8)
       where the commutator is defined via

       $[A,B]:=AB−BA.$

       (9)

The usual statement of Newton's laws^53,54 discusses the kind of equations of motion that are allowed rather than the underlying mathematical assumptions. This does not make it easy to compare the underlying mathematical assumptions of quantum and classical mechanics. Perhaps, this is one of the reasons that quantum mechanics took some time to be developed and understood. If one were to attempt to write the postulates of classical mechanics in a form that parallels that of quantum mechanics, then one might arrive at the following.

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 replaced by sums.

1. Physical states.

   • Pure physical states are points in a set, $z∈m$⁠, called the phase space.

     • Phase space is typically a continuous and smooth space such as a differentiable manifold, but it can also be a discrete set.

   • Mixed physical states are represented by a probability distribution function, $p(z)$⁠, over phase space.

     • Probability is normalized to unity $∫pdz=1$⁠.

   • The phase space of a composite system is the tensor product of the phase spaces of the subsystems.

     • The dimension of the composite phase space is the product of the dimensions of the subsystems.

2. Measurement of physical observables.

   • Physical observables are real functions of coordinates in phase space, $O:m→ℝ$⁠.

   • In principle, measurements of the coordinates, z, can be made with arbitrary precision, and, so, a measurement of an observable at the point z yields its value, $O(z)$⁠.

   • The probability of measuring a particular value of the observable, $O(z)$⁠, is given by the probability, $p(z)$⁠, of observing the pure state at the point z. Therefore, the expectation value of an observable is given by

     $⟨O⟩=∫Opdz.$

     (10)

3. Temporal evolution.

   • General open systems

     •
     • Physical states evolve via (generically nonlinear) coordinate transformations, $ζ:m×i→m$⁠,

       $z(t)=ζ(z(t0),t).$

       (11)
     • Continuous time evolution is generated by (generically nonlinear) differential equations

       $dz/dt=V(z,t).$

       (12)

   • Ideal closed systems

     • Ideal closed systems are Hamiltonian systems that evolve via Hamilton's equations of motion, which yield infinitesimal symplectic transformations of phase space.

     • In canonical coordinates, $z={q,p}$⁠, composed of equal numbers of configuration space coordinates, q, and momentum space coordinates, p, Hamilton's equations of motion for Hamiltonian $H(z)$ are

       $dq/dt=∂pH dp/dt=−∂qH.$

       (13)

     • Evolution over a finite time interval yields a (generically nonlinear) symplectic transformation of phase space, $ζ:m×i→m$⁠.

Many physicists and mathematicians are familiar with the temporal evolution of ideal closed systems, e.g., in the classical CM Postulate 3b and quantum QM Postulate 3b contexts. The more general formulation for open classical systems in CM Postulate 3a is also quite natural. This form applies to non-ideal open classical systems, which typically experience various forms of dissipation, such as friction, diffusion, and collisions.

In contrast, the discussion of quantum dynamics in QM Postulate 3a was far too brief to do justice to this important topical area. While the mathematical justification is on solid footing today, in terms of applications to real-world quantum hardware platforms, this field is still under active development.

For closed ideal quantum systems, QM Postulate 3b requires time evolution to be unitary and, hence, reversible. However, similar to the way that classical chaos can scramble classical information, quantum dynamical processes can scramble quantum information. Classically, the combination of the butterfly effect (positive Lyapunov exponents) and coarse-graining due to the finite precision of measurements leads to irreversibility. However, quantumly, the unitary dynamics only has vanishing Lyapunov exponents, and, so, while the evolution may track the semiclassical dynamics for long periods of time, any apparent exponential instability associated with the Lyapunov exponent must always eventually cease.

The paradox is resolved by the fact that most quantum systems are actually open; i.e., the observer does not have access to information to all parts of the system. The parts of the system that either are not or cannot be measured are usually referred to as the environment. The existence of this unmeasured or hidden information has effects that are similar to coarse-graining. Thus, open quantum systems can display many different types of non-ideal and non-unitary behavior.

In fact, the act of measurement is the primary example of a non-ideal but linear evolution process. Because the measurement apparatus in QM Postulate 2 is assumed to be classical, the probability of each eigenvalue is described by classical probability theory. Thus, even if one begins with a pure state, after the measurement is performed, the resulting knowledge about the system must become a mixed state: a probability distribution over possible eigenfunctions of the observable operator that was measured.

Decoherence refers to the way in which small but persistent interactions with the environment cause an almost irretrievable loss of information for open quantum systems. The phrase decoherence was introduced by Zeh⁵⁵ and has been explored in great detail by Zurek^56,57 and collaborators. Physicists and mathematicians, such as Choi, Gorini, Krauss, Kossakowski, Lindblad, Redfield, Stinespring, and Sudarshan, explored possible types of evolution laws, called master equations, that an open quantum system can undergo.⁵² The goal is to develop a master equation that describes the loss of information in a manner similar to the way the Fokker–Planck equation describes the loss of information due to diffusive processes for classical probability theory.

If one were to use an accurate quantum model of the environment, then the state of the system can interact in many different ways with the environmental degrees of freedom. For simplicity, it is customary to attempt to determine a reduced order model that has the expected qualitative behavior. For example, it is customary to model the interactions with the environment as a quantum process: a linear transformation of the density matrix that has no memory, i.e., a Markovian process. Note, however, that using a more realistic reduced order model of the environment that has nontrivial dynamics will typically generate a nonlinear model that is not Markovian.

Linear operators that act on the density matrix, rather than on pure states, are referred to as quantum processes, quantum channels, or superoperators. If the dimension of the Hilbert space is N, then there are $N2$ unitary operations that describe ideal evolution. However, there are many more, $N4−N2$⁠, quantum processes that describe decoherence.

If the density matrix is to remain a positive Hermitian operator, then the superoperator must itself be Hermitian and positive. In fact, it is common to demand a stronger criterion called completely positive (CP)^58,59 that results from considering the evolution of the system when coupled to a reservoir that represents the environment with a Hilbert space of arbitrary dimension, d. In this case, the state is $ρ⊗1d$⁠, and evolution of the subsystem with superoperator, $e$⁠, while the reservoir is unchanged results in the state $e(ρ)⊗1d$⁠. The map, $e(ρ)$⁠, must be completely positive in order for the coupled map, $e(ρ)⊗1d$⁠, to be positive for the total system ⊗ reservoir Hilbert space. It turns out that complete positivity is a stronger condition than positivity alone as there are examples of maps, e.g., the transpose operation, that are positive but not completely positive. In order to conserve probability, the superoperator must also be trace-preserving (TP). Such superoperators are called completely positive trace-preserving (CPTP) operators.

A Hermiticity preserving superoperator must have the form^52,59

where ${Lj}$ defines a complete basis of operators of size $N×N$ and the matrix χ_ij is Hermitian. In order for the map to be CP, the matrix χ_ij must also be positive semi-definite. (Note that completely positive refers to the action on the total Hilbert space rather than the definiteness of the eigenvalues.) In order to be trace preserving (TP), the condition

must hold. The trace preservation criterion represents $N2$ constraints. Since there are $N2$ unitary operations, this implies that there are $N4−2N2$ non-unitary decoherence operations.

Because the $χjk$ matrix is Hermitian, this form can be simplified by diagonalizing this matrix with a unitary transformation. Applying this transformation to the operator basis allows one to find a new basis, ${L̂α}$⁠, that simplifies the equation to a sum over the eigenvalues of the $χjk$ matrix, which we denote as $χα$

Because the eigenvalues $χα$ are positive, the so-called Krauss form⁶⁰ can be achieved by absorbing the square root of the eigenvalue into the definition of the operators $L̂α→χα1/2L̂α$⁠.

Note that completely positive only requires the eigenvalues to be positive semi-definite $χα≥0$⁠. If the map is positive but not completely positive, then it must have at least one negative eigenvalue.

The Gorini–Kossakowsi–Lindblad–Sudarshan (GKLS) master equation^61,62 represents the most general form of a continuous CPTP evolution of the density matrix when the system is separable from the environment and when their interactions are Markovian. Just as the Fokker–Planck equation describes processes such as diffusion due to collisions for the Liouville evolution of the classical PDF, the Lindblad equation adds nonunitary CPTP processes to the von Neumann equation for the density matrix. The GKLS equation is given by

where the anticommutator is defined via ${A,B}:=AB+BA$⁠. Thus, in addition to unitary evolution determined by the Hamiltonian H, one adds a general CP evolution as well as the TP damping term

required to preserve the trace. In this case, the collapse operators, ${Lj}$⁠, define a basis over trace-free matrices of size $N×N$⁠. (The freedom of including the unit matrix, which possesses a trace, was already used in the construction of the damping operator A.) This evolution is CPTP as long as the rate matrix, ν_ij, is Hermitian and positive semidefinite. Information loss occurs whenever the rate matrix has positive eigenvalues.

Again, because the rate matrix is Hermitian, this form can be simplified by diagonalizing $νjk$ with a unitary transformation. Applying this transformation to the collapse operators allows one to find a new basis, ${L̂α}$⁠, that simplifies the equation to a sum over the eigenvalues of the χ_ij matrix, which we denote as $να$

Since the eigenrates $να$ are positive, the square root of the rates can be absorbed into the definition of the collapse operators, as in the Krauss form above. While the diagonal form offers simplicity, it hides the complexity associated with the myriad ways in which the environment can potentially interact with a quantum system to form an arbitrary rate matrix. A pedagogical introduction to the GKLS master equation is given in Ref. 63.

There are a few important environmentally induced decoherence processes that affect almost all quantum systems: (1) relaxation: of higher energy levels to lower energy levels, (2) excitation: of lower energy levels to higher energy levels, and (3) dephasing: the diffusion of the complex phases within a linear superposition. For a harmonic oscillator, the relaxation collapse operator is the destruction operator, a, the excitation collapse operator is the creation operator, $a†$⁠, and the dephasing collapse operator is the number operator, $N=a†a$⁠. While excitation is important in general contexts it is often much slower than the other two processes for typical quantum hardware platforms.

Just as the bit is the smallest unit of classical information, the 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 states: 0 or 1. In the language of quantum information, we might call these classical states $|0⟩$ and $|1⟩$⁠. However, it would be clearer to label these states as density matrices $|0⟩⟨0|$ and $|1⟩⟨1|$ to signify the fact that, in the classical world, superpositions are not allowed and, hence, these two states are mutually exclusive. In other words, in the classical world, only one of the two possible pure states $|0⟩⟨0|$ and $|1⟩⟨1|$ is allowed to exist at any given moment in time.

To make a physical bit from a physical device, one must devise a physical system that is subject to a potential that has two deep wells, as shown in Fig. 1. Then, one must ensure that the temperature of the system and, more generally, the possible range of energy exchanges with the environment remains far below the energy barrier separating the two wells.

If one constructs a similar potential with a quantum mechanical device, then the wavefunction can potentially be in a superposition of all accessible states, including the lowest energy level, called the ground state or vacuum state, as well as the higher energy levels, called the excited states, also shown in Fig. 1. To make a physical qubit from a physical device, one must confine the wave function to two states. For this two-state system, $|0⟩$ is the ground state and $|1⟩$ is the first excited state. Again, one must control the temperature and the interactions with the environment so that the probable range of energy exchanges with the environment is much less than the energy differences between these two states as well as the energy differences between these two states and all other states. Moreover, one must also be careful to control environmental interactions that cause decoherence and destroy superposition and entanglement.

The wavefunction of a qubit lives in a two-dimensional complex Hilbert space

Due to the fact that the wavefunction is a ray in Hilbert space, i.e., it actually lives in projective Hilbert space, the overall complex constant is unimportant and is normalized to unity. Note, however, that for describing a normalized vector in Hilbert space, and for describing unitary transformations of physical states, one must include the overall complex phase factor, $eiχ|ψ⟩$⁠.

Thus, the qubit is described by two real parameters that form the surface of a sphere, called the Bloch sphere. By convention, the $+ẑ$ direction that points toward the top of the sphere corresponds to the $|0⟩$ state while the $−ẑ$ direction that points toward the bottom of the sphere corresponds to the $|1⟩$ state. Thus, an arbitrary qubit state can also be defined by the direction that it points in; if the unit direction vector is $n̂$⁠, then the state is denoted $|n̂⟩$⁠. For example, $|ẑ⟩:=|0⟩$ and $|−ẑ⟩:=|1⟩$⁠, but remember that these states are orthogonal. Similarly, one can define the states

Given the definition above, when measurements are made along the $ẑ$ axis, the probability of $|0⟩⟨0|$ is $p0=cos2(θ/2)=(1+cos θ)/2$ and the probability of $|1⟩⟨1|$ is $p1=sin2(θ/2)=(1−cos θ)/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 $p1S$ and variance $p0p1S$⁠, so the variance of the mean decreases as the central limit theorem would predict, $p0p1/S$⁠. In agreement with QM Postulate 2: immediately after a measurement, the state is known with certainty, so the variance must vanish. As $S→∞$ for fixed probability, 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→∞$⁠, but the product $λ=p1S$ is held fixed, then the binomial distribution tends toward the Poisson distribution with mean and variance λ.

There is also the possibility of rotating the measurement basis to point along the $x̂$ or $ŷ$ directions. In this case, the probability is controlled by the angle $ϕ$ rather than θ. 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 $±x̂$ axis, simply rotate the state counterclockwise around the $ŷ$ axis by $∓90°$⁠.

It is useful to explicitly define the effect of a rotation as a unitary operation on the qubit. Define the Pauli vector as the vector of Pauli matrices $σ→=(σx,σy,σz)T$⁠. A rotation around the unit vector $n̂$ by the angle θ is defined by

Note that a rotation by $θ=2π$ multiplies the qubit by the overall phase $e−iπ=−1$⁠. This explicitly demonstrates the fact that the special unitary group, SU(2), is actually a double cover of the special orthogonal group of rotations SO(3).

A classical probability distribution over qubit wavefunctions is equivalent to a positive Hermitian matrix called the density matrix

Due to Hermiticity, the off diagonal elements must be complex conjugates, $ρ10*=ρ01$ and the diagonal elements must be real. These diagonal elements represent the probability of measuring the two states, $|0⟩⟨0|$ and $|1⟩⟨1|$⁠, and are guaranteed to be positive due to the positivity of the density matrix. In order for the probabilities to sum to one, there is the constraint $Tr ρ=ρ00+ρ11=1$⁠. The off diagonal elements represent non-classical superposition states that have no analog in the classical world.

Again, due to Hermiticity, the density matrix can be diagonalized, so that it has the form

where $|±n̂⟩$ are orthogonal states and $p±$ are the classical probabilities of each of these states occurring. Any 2 × 2 density matrix can also be written in the form

Thus, the qubit density matrix can be visualized as an arrow of length $|p+−p−|$ pointing in the $n̂$ direction, called the Bloch vector. A proper mixed state must have $|p+−p−|<1$⁠. A pure state can only have eigenvalue 1 for one entry and eigenvalue 0 for the other, so that $Tr ρa=1$ for any a > 0. In contrast, for a proper mixed state, both eigenvalues must be less than 1, so that $Tr (ρa)<1$ for any power a > 1 [and $Tr (ρa)>1$ for any power $0<a<1$].

Given the discussion above, a classical probability distribution function (PDF) over classical bits is simply a diagonal density matrix

Because it is diagonal, it can also simply be considered to be a function of the classical states alone, $p:{0,1}→ℝ$⁠, with the two values $p(0)=p0=ρ00$ and $p(1)=p1=ρ11$⁠. Again the normalization condition, $Tr (p)=1$⁠, ensures that the probabilities sum to unity.

In the following, we assume that individual qubit states can be physically prepared and measured in the same basis each time, which is referred to as the computational basis. Two qubits live within a four-dimensional Hilbert space, which is the tensor product of the two individual qubit Hilbert spaces. Hence, basis states that span the two qubit Hilbert space can be denoted $|q1q2⟩=|q1⟩⊗|q2⟩$⁠, where $q1,2$ refer to the state of qubit 1 and 2, respectively. There are four possible measurement outcomes, which correspond to four possible classical bit strings. However, there are many more quantum states that are not classical, corresponding to three complex degrees of freedom.

In addition to single-qubit superpositions, the key resource to be exploited is the ability to create entangled states. Entangled states are defined as states that cannot be represented as tensor products of the states of the individual subsystems. Thus, a state such as $|x̂x̂⟩=(|00⟩+|01⟩+|10⟩+|11⟩)/2$ is not entangled, even though each qubit is in a superposition state itself. In contrast, a two-qubit state such as the Bell state $(|00⟩+|11⟩)/21/2$ is truly entangled. When quantum information is entangled, but only part of the system can be measured, this leads to kind of information hiding that causes decoherence. Note that the concept of entanglement explicitly refers to consideration of the tensor product structure of the Hilbert space and, hence, the basis that is being used.

In order to find the size of a multi-qubit wavefunction, simply apply QM Postulate 1c: multiply the dimensions of the Hilbert spaces of the subsystems and apply the ray condition. As illustrated in Fig. 2, for n qubits, there are $N=2n$ computational basis states, which is equivalent to the number of possible measurement outcomes, and, hence, to the number of possible classical bit strings. For n qubits, the wavefunction is specified by $N=2n−1$ complex numbers. Although classically, the state is only ever in one of the possibilities, a classical probability distribution (PDF) of n bits is specified by $N=2n−1$ real numbers. Hence, the quantum wavefunction holds exactly twice the information of the classical PDF due to the fact it has a complex phase associated with each probability amplitude. Yet, the classical PDF should really be compared to the quantum density matrix, which is specified by $N2−1=22n−1$ real numbers and, thus, has quadratically more information.

In many ways, these simple facts are the reason for all of the excitement about quantum computing. This capacity of a quantum memory register for holding an exponentially large amount of information and the potential ability of a quantum computer to efficiently manipulate this exponentially large amount of information is amazing from the standpoint of classical computation.

This also implies that performing a direct classical simulation of a multi-qubit quantum system is exponentially hard! Although many ingenious algorithms have been invented for approximately solving these kinds of problems, progress with direct solvers is limited due to the exponentially large size of the required Hilbert space. Classical computations much use drastically reduced models to even attempt the simulation of large systems. Now we can understand the first part of Feynman's statement: If you want make a (quantum-mechanical) simulation, you'd better make (the computational engine) quantum mechanical… so that it has a sufficiently large Hilbert space and so that this Hilbert space can be manipulated efficiently.

Note that classical algorithms that can efficiently approximate the classical PDF are also quite powerful. From the simple scaling of the amount of information in the classical PDF vs the quantum density matrix, one might expect that quantum algorithms can generically outperform classical probabilistic algorithms by this same quadratic factor. As will be discussed in detail later, this is indeed the case.^37,64

Just as classical information with d possible values is referred to as a dit, the quantum analog of a wavefunction with d possible states is called a qudit. For both the classical and quantum cases, systems with multiple bits and dits can be used to simulate each other with polynomial efficiency, so there is no change in computational complexity class. In the classical world, transistors are so cheap that essentially all modern computing technology is based on bits. In the quantum world, experimentalists are still striving to reduce the error budget to the point at which error correction can be applied. For this reason, a number of researchers have proposed using the higher energy levels of a quantum system to encode information, as illustrated in Fig. 3. For certain quantum device platforms, these higher energy states are naturally present, so it could be a good idea to effectively pack more qubits into a single device, potentially allowing less error-prone communication between these states.⁶⁵ However, since decoherence rates are often faster for higher energy levels, there is clearly a trade-off in performance that must be investigated for each experimental hardware platform.

Universal quantum computation^5–7,52 refers to models of quantum computing that can efficiently simulate a universal classical computer, i.e., a Turing machine. There are a number of different models of universal quantum computation that have been proven to be equivalent to one another in computational capabilities; i.e., they can simulate each other with polynomial complexity. Popular models include the digital circuit model,^5–7 the quantum walk model,^66,67 the adiabatic model,⁶⁸ the quantum annealing model,^69,70 and the measurement-based model.⁷¹ The digital circuit model, also known as the digital quantum computational model, has power and simplicity and is perhaps closest to digital classical computation based on classical circuits, which is familiar to most scientists and engineers today.

The basic idea of the digital circuit model is that quantum states can be transformed efficiently via linear unitary operations: $ψ=Uψ0$⁠. In fact, according to QM Postulate 3b, this is the only operation that can be performed without making a measurement. After making a measurement, the state generically becomes a density matrix $ρ$⁠, but in between measurements, linear unitary operations can be performed: $ρ=Uρ0U†$⁠. The fact that a quantum computer can do this efficiently is unbelievable from the classical standpoint, because it allows one to efficiently perform matrix–vector operations with exponentially large matrices and exponentially large vectors.

While there are a large number, $N2$⁠, of independent unitary operations, they are generated by a small, $∼O(n)−O(n2),$ number of basic operations called a gate set. Single qubit operations can be performed efficiently with just two standard gates that perform rotations of the Bloch sphere. Rotations by angle θ around the $x̂, ŷ$⁠, and $ẑ$ axes of the Bloch sphere, are denoted $RX(θ), RY(θ)$⁠, and $RZ(θ)$⁠, respectively. A rotation by angle $θ=π$⁠, simply generates the Pauli operations, denoted $X=σx=NOT, Y=σy$⁠, and $Z=σz$⁠. A standard result in Lie group theory says that only two such operators are required to generate the whole group. Another important example is the Hadamard gate

Thus, starting from the $|ẑ⟩:=|0⟩$ state pointing along the $ẑ$ axis, this generates the uniform superposition state $|x̂⟩$⁠. Similarly, starting from the state $|−ẑ⟩=|1⟩$⁠, this generates the superposition state $|−x̂⟩$⁠.

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

Another useful operation is the SWAP operation that simply exchanges information between qubits

If $U2×2$ is a single-qubit unitary and $12×2$ is the single-qubit unit matrix, then a controlled U operation has the form

This shows that the operation U is performed if the first qubit is in the $|1⟩$ state. To make this clearer, use the notation $CUjk$ to indicate the fact that the state of the kth qubit is controlled on the state of the jth qubit.

A sequence of unitary operations can be visualized through a circuit diagram. Operations are applied from left to right in a circuit diagram—which is opposite to the usual right to left operator ordering convention in mathematical expressions. A circuit diagram has a number of lines going from left to right that each indicate individual qubits. Blocks in the middle of the diagram represent single-qubit operations. Blocks that are controlled on another qubit have a line from the block to the control qubit. For example, consider the creation of the entangled Greenberger–Horne–Zeilinger (GHZ) state⁷² on three qubits defined by $(|000⟩+|111⟩)/21/2$ (an analog of the Bell state for two qubits). This state can be created with the circuit diagram shown in Fig. 4. The instructions are to (i) apply H to the first qubit; (ii) perform a $CNOT1,2$ operation, where the state of the second qubit is controlled by the state of the first qubit; (iii) perform a $CNOT1,3$ operation, where the state of the third qubit is controlled by the state of the first qubit; and (iv) finally, the quantum state can be verified by measuring all qubits. If one rotates the states before measurement, then one can test the consequences of Bell's theorem.⁷³ 

The fundamental discovery of early research on quantum algorithms is that many useful computations can be performed using relatively few $O(poly (n))$ of these basic gate operations.^5–7 Key resources are the ability to create superpositions and entanglement and to use interference to one's advantage. This is often called quantum parallelism. The second fundamental discovery is that any reversible classical computation can also be performed.^5–7 This requires first computing and storing the result and then uncomputing or reversing the order of operations on the original memory register.⁷⁴ 

However, emulating a classical computer usually comes with non-negligible overhead, and so, one might prefer to use a classical computer to perform certain operations such as classical logic. For this reason, the designers of quantum computer platforms usually consider a hybrid architecture where classical processing units (CPUs) sit beside as well as drive quantum processing units (QPUs) much like today's heterogenous CPU and graphical processing unit (GPU) architectures. Like CPU/GPU platforms, there will be a large overhead in input/output (I/O) operations that exchange data between CPUs and QPUs. In fact, the I/O will be exponentially large if one attempts to exchange the entire quantum random access memory (QRAM). In today's world, QPUs are also far more expensive than CPU's so hybrid architectures are an obvious choice. While it may be difficult to envision today, perhaps at some point in the future, quantum computers will become useful enough and cheap enough that they will become ubiquitous.

Now it is very important to discuss some key limitations of the universal quantum computing paradigm.^5–7 

• Approximating an arbitrary unitary is exponentially hard.

  This implies that only certain classes of unitaries, precisely those that are composed of $O(poly )$ basic gate operations, can be constructed “easily.”

• Initializing all quantum information is exponentially hard.

  Initializing quantum information requires applying an arbitrary unitary operation (usually starting from the ground state). So, this also implies that initializing all information in quantum memory is exponentially hard.

• Measuring all quantum information is exponentially hard.⁶ 

  Measuring all quantum information requires making an exponential number of measurements. Moreover, each measurement generically requires a different unitary transformation to be applied.

• 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 measurements of the state of a qubit is given by the classical binomial distribution. 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, the relative 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 exponentially hard.⁶ 

  Since these values would only be observed with exponentially small probability, one must make an exponentially large number of measurements to observe the values with any accuracy.

• Nonzero failure probability.⁷ 

  Quantum algorithms are probabilistic and, as such, they have a certain probability of failure. This failure probability can be made quite small by batching many runs together (see the clear explanation of the Powering Lemma in Ref. 37). However, generically, the failure probability does not vanish and this feature is quite different from classical deterministic algorithms.

Thus, it is misleading to say that the universal quantum computation model offers a method of efficiently simulating any quantum system. In fact, in order to remain efficient, this model can only approximately simulate certain quantum systems and it can never actually exploit the entire Hilbert space. Although we began with the notion that there is an exponential amount of information that can be stored and manipulated, we now see that the actual amount of such information that can be manipulated and measured is much smaller than one might have hoped.

In fact, it is useful to compare the power of the quantum computing model to classical probabilistic models of computation, such as Monte Carlo methods. For probabilistic models, one never manipulates all of the information required by a direct calculation; rather, one relies on the central limit theorem to provide a slow but steady increase in relative accuracy as the number of individual samples is increased.

The quantum computational model is certainly a probabilistic model and many runs of the same algorithm are needed to yield a final answer. After the required measurements are made for each run, one can associate the outcome with a quantum trajectory that is consistent with the outcome.⁷⁵ It is only by assembling a large number of quantum trajectories that an accurate answer is achieved. In a certain sense, the key difference is the fact that quantum trajectories can exploit superposition states and interference effects rather than classical trajectories which are limited to computational basis states alone.

Nonetheless, there are many elegant mathematical theorems that prove that the quantum computational model can exceed even the classical probabilistic counterpart.^5–7 Generically speaking, the ability of the quantum density matrix to store quadratically more information than the classical PDF enables up to a quadratic speedup over classical probabilistic methods.⁷⁶ However, there are certain classes of problems for which superpolynomial and even exponential speedup has been proven to hold.⁷⁷ 

From the inception of quantum computing, it was recognized that the ability to correct quantum errors is indispensable for building large-scale quantum computers. Quantum error correction (QEC) protocols^5–7 seek to undo errors that are induced by the inability to enact the desired unitary operations with perfect precision as well as the uncontrollable interactions with the environment that cause decoherence. The general strategy is to use an encoding that repeats information so that it can be checked against various types of error syndromes. This implies that there is a distinction between the number of physical qubits that comprise the quantum hardware platform and the number of logical qubits that can encode quantum information for a given time period with negligible error rates. The subject of QEC is vast and we cannot do justice to the wide range of ideas necessary to understand the contemporary approaches to the field. The interested reader can begin by consulting standard textbooks,^5–7 books on quantum error correction,^78–80 and relatively recent review articles.^81–85 

In general, quantum error correction protocols seek to undo errors that are induced by the inability to enact the desired unitary operations with perfect precision as well as the uncontrollable interactions with the environment that cause decoherence. As discussed in Sec. II B, if the system of interest has N states, then, each time one attempts to enact a quantum gate, there are $N4−N2$ ways, including both decoherence operations and accidental unitary operations, that the intended gate can go awry. If one is attempting to obtain an exponential speedup by utilizing a large number of qubits, so that $N=2n$⁠, then there are an exponential number, $∼N4=24n$ of gate features that must be controlled. Just obtaining enough measurements to fully characterize the operation of the gate then becomes an intractable problem. Now we can understand the second part of Feynman's quote: by golly it's a wonderful problem because it doesn't look so easy.

For example, from the point of view of gate design, it is very efficient to have all-to-all connectivity between qubits, as is the case for ion trap platforms.⁸⁶ This eliminates the need for many extra SWAP gates that would otherwise be needed to shuttle quantum information between different parts of the memory register. Because these SWAP gates must often be performed in serial, eliminating them also increases the potential for parallelization of multiple gate operations.

However, from the point of view of error correction, it is far superior for qubits to be well-separated both in terms of their resonant frequencies and physical locations. For example, this can be accomplished by carefully choosing and tuning the resonant frequencies and by physically arranging qubits into a lattice, as is the case for superconducting platforms.⁸⁷ Another example is the quantum charge coupled device (QCCD) ion-trap architecture,⁸⁸ which manages its error budget by moving inactive ions to distant locations from active gate regions. If one can limit the amount of significant crosstalk between qubits to m nearest neighbors alone, then the number of gate features that need to be controlled is reduced to $24m$ per qubit. The key point is that now m is independent of the size of the memory register.

Known quantum error correction algorithms can be proven to be successful, allowing one to store and manipulate quantum information indefinitely, as long as the error rates are below a specific limit depending on the protocol, typically at $10−3−10−5$ per operation. For the best performing approaches known today, e.g., the topological surface code⁸⁹ with magic state distillation,⁹⁰ the number of physical qubits necessary to represent a logical qubit is estimated⁸¹ to be on the order of $103−104$⁠. Hence, the number of qubits required to perform a useful calculation is estimated to be on the order of $105−106$ or more.⁸⁴ Unfortunately, this large overhead cost has been investigated and interpreted⁹¹ as implying that the first successful error-corrected quantum applications should offer speedups that are beyond quadratic, e.g., at least quartic.

After decades of theoretical development of error correction codes and experimental demonstrations of primitive steps, a surface code was recently demonstrated in a superconducting circuit,⁹² and a color code was recently demonstrated on a QCCD ion-trap quantum computer.⁹³ These recent demonstrations are starting to approach the threshold where logical qubits outperform physical qubits.

Since each operation takes a certain physical time period, the limit on the error rate can also be translated to a limit on the norm of the effective decoherence rate matrix, $||ν||$⁠. Thus, the basic idea of quantum error correction is to increase the information confinement time until this limit can be achieved. Perhaps surprisingly, this goal is very similar to the Lawson criterion, which requires a net energy-producing fusion reactor to achieve a specific energy confinement time. Within the field of fusion energy, in addition to inherent physical limitations of the reactor design, it is anomalous and often turbulent transport processes that typically control the energy confinement time. Similarly, within the field of quantum computing, in addition to inherent physical limitations of the hardware platform, it is decoherence due to anomalous interactions with the environment that controls the information confinement time. Thus, it will be interesting to gauge the relative progress in the two fields by measuring the rate of improvement of the respective confinement times.

The field of quantum algorithms has seen a rapid development in interest over the last two decades and it is not possible for a short introduction such as this to do justice to the wide range of brilliant ideas that have begun to bear fruit. Nonetheless, when compared to classical algorithms, the number of essentially quantum subroutines are actually relatively few. In part, it is the vast number of ways that these subroutines can be configured as well as the number of applications that these quantum subroutines can speed up that leads to the considerable complexity in the field. Moreover, since the field is still undergoing rapid development, there may be many new concepts waiting to be discovered.

The challenge of performing computations using an exponentially large quantum Hilbert space was reconceived as a potential opportunity by Feynman^50,51 and Manin⁴⁹ in the early 1980s. Feynman reported that his interest in the subject was piqued by Charles Bennett, who had worked to understand reversible classical computation⁷⁴ and who wanted to understand how the physics would work when using the laws of quantum mechanics.⁵¹ In fact, it was Charles Bennett and Giles Brassard that made the first inroads into proving rigorous theorems that demonstrated that quantum communication channels had more information processing capacity than classical communication channels.^6,7

In the 1980s, Benioff⁴⁸ and Deutsch⁹⁴ defined quantum Turing machines and showed that that they were as powerful as classical Turing machines. In 1992, the Deutsch–Josza algorithm⁹⁵ was discovered, a concrete example of a toy problem that could be solved faster with a quantum computer than a classical one. Soon after, the Bernstein–Vazirani algorithm⁹⁶ and Simon's algorithm,⁹⁷ were discovered as well. It turns out that these simple algorithms rely on the ability to perform the Fourier transform more efficiently on a quantum computer.

In 1994, Peter Shor's algorithm⁹⁸ for factoring integers and taking discrete logarithms exponentially faster than a classical computer surprised the world with the realization that quantum computers can break classical security protocols relatively easily. The key step is estimating the periodicity of a sequence with unknown period and the key discovery is that there is an exponentially efficient factorization of the discrete Fourier transform on quantum memory registers—today called the quantum Fourier transform (QFT). After using the QFT, the period can be estimated using phase estimation protocol to estimate the eigenvalue of a Hamiltonian, given an eigenvector. Shor then went on to proving that there were quantum error correction protocols that could be used to protect quantum information,⁹⁹ which gave the abstract theory a sense of physical realizability, perhaps for the first time.

The next breakthrough was Lov Grover's search algorithm^100,101 which speeds up unstructured search by the square root of the number of items. Key elements of the search algorithm were then separated into algorithms for amplitude amplification, quantum counting, and amplitude estimation.^102–104 These subroutines form the basis for quantum algorithms for computing sums and integrals^64,105–107 and forms the basis for quadratically speeding up classical Monte Carlo algorithms.³⁷ 

In the late 1990s, Seth Lloyd and Daniel Abrams discovered efficient algorithms for simulating $k−$local Hamiltonian interactions.^108,109 The basic idea is to split the Hamiltonian into a relatively low number of terms that do not commute with each other and to then use the Lie–Trotter–Suzuki decomposition, aka Trotterization, to approximate the exact unitary as a product of terms. For example, for Hamiltonians defined on a 1D lattice with nearest neighbor interactions, one can split the Hamiltonian into a sum of three pieces, wherein all terms within each piece commute: a diagonal piece, an even index off diagonal piece, and an odd index off diagonal piece. As another example, the single particle Schrödinger equation can be simulated efficiently¹¹⁰ in a manner that mimics the phase space path integral by splitting the kinetic and potential energy terms into two pieces. Then one can repeatedly use the QFT and inverse QFT to switch back and forth between the momentum and position basis, so that each term is diagonal in the respective basis. While operator splitting methods, such as Strang splitting and higher order Suzuki methods,¹¹¹ are familiar to applied mathematicians, recent work¹¹² has led to new understanding of error convergence for these methods.

Quantum walks¹¹³ are a vast generalization of the amplitude amplification paradigm that allow discrete state changes based upon the outcome of flipping a quantum coin or rolling quantum dice; i.e., on the essential use of ancillary qubit registers. Quantum walks were originally constructed for the purpose of understanding simple models of quantum physics. Notably, Richard Feynman developed a quantum walk model for how Dirac fermions traverse a discrete lattice.¹¹⁴ By the early 2000s, the essential ideas behind quantum walks were codified into a powerful set of algorithms¹¹³ and were then proven to offer a universal computing model of their own.^66,67

Eventually, general-purpose Hamiltonian simulation algorithms (HSAs) that exploit sparsity in the Hamiltonian itself were developed.¹¹⁵ These algorithms were then improved by using quantum walks,^116,117 new representations of Hamiltonians as a linear combination of unitaries^118,119 (LCoU), and as unitary approximations of the truncated Taylor series of the propagator.¹²⁰ The latest development of these algorithms eventually achieved nearly optimal dependence on all parameters.^{21,22,121,122}

Hamiltonian simulation powers the Harrow, Hassidim, and Lloyd (HHL) sparse quantum linear solver algorithm¹²³ (QLSA) as well as improved versions with superior performance.^124–127 Next, sparse quantum linear differential equation solver algorithms (QLDSA) were developed that use linear solvers in combination with sparse Hamiltonian simulation algorithms.^18,128

Today, there are new classes of optimal Hamiltonian simulation methods powered by quantum signal processing (QSP)²² and qubitization.¹²⁹ The term qubitization refers to the use of an ancillary qubit, as in quantum walks, while the term quantum signal processing refers to the way in which this ancillary qubit is manipulated to perform useful calculations. The combination of these methods can generate nonlinear quantum eigenvalue transformations (QEVT) of a block-encoded operator that is a normal matrix, $[H,H†]=0$⁠. This approach was then generalized to generate nonlinear quantum singular value transformations (QSVT) for an arbitrary matrix.¹²⁶ The use of QEVT and QSVT was shown to provide a grand unification of many quantum algorithms,^122,126 including search, amplitude amplification. matrix multiplication, matrix inversion, phase estimation, and QFT. QSP has already been put to good use within quantum algorithms designed to solve plasma physics problems.^23,130

There are a number of useful textbooks^5–7 and review articles^131–133 on quantum algorithms. Although the quantum computational model might take some time to become familiar with and new quantum algorithms with improved complexity are under rapid development, it is helpful to realize that there are actually only a few essentially quantum subroutines that power the vast majority of quantum algorithms.

In the 1990s, it was discovered that the Fourier transform can be factorized in an extremely efficient manner on a quantum memory register. Many Hamiltonians of interest can be simplified and, hence simulated efficiently, if the Fourier transform is exponentially cheap and if one does not need to extract intermediate information. In fact, the key subroutine that powers Shor's algorithm for factoring integers and taking discrete logarithms is the Fourier transform.

There are many variations of the Fourier transform that arise from the character theory of commutative groups. The most familiar of these is the Walsh–Hadamard transform, which arises from the group $(ℤ/2ℤ)n$⁠. The Walsh–Hadamard transform on multiple qubits, $Hn$⁠, is simply the tensor product of Hadamard gates

The subscript n is not necessary, because the number of qubits to be acted upon should be clear from the context. Evidently, this transform factorizes so well that it can be performed with single-qubit rotations alone. It is simple to show by induction that this generates a uniform superposition of states from vacuum

Then, a quantum subroutine can work on the entire superposition at once: $ψ=U|s⟩$⁠. This simple initialization step is used by many quantum algorithms.

The discrete Fourier transform arises from the modular group $ℤ/Nℤ$ and is defined via

For the binary version of the algorithm, $N=2n$⁠, and, one can express the index in binary notation via

The quantum Fourier transform (QFT), which is nothing but the usual Fourier transform, factorizes into the following form:

Now the state repeated within the tensor product is simply $R(2πj2−ℓ)H|0⟩=R(2n−ℓxθ)H|0⟩$⁠, where $θ=2π/N=2π/2n$ and the single qubit phase rotation operator is defined as

Thus, one can also write

The effect of the factor $2n−k$ is to truncate the phase to the last n–k binary digits of x, i.e.,

so that

Now, it is easy to see that each rotation should itself factorize as the product of $n−k−1$ controlled rotations, $CRℓ,k$⁠, between the kth qubit and all higher $ℓ$ th qubits. The operator

is almost correct, but it applies the phase to the kth qubit rather than the n–kth qubit and the desired rotations were supposed to be implemented in the opposite order within the tensor product. This is easy to correct by first computing the result for both k and n–k and then swapping the results. In other words, simply add a final stage of SWAP gates at the end defined by

The final factorized expression is

The main QFT stage has n single-qubit Hadamard gates and $n(n−1)/2$ controlled two-qubit rotations, for a total of $n(n+1)/2$ operations. The final SWAP stage consists of $⌊n/2⌋$SWAP operations. The swap stage can often be eliminated with a logical reordering of qubits in the next subroutine or by reordering the measurements that read out the ancillary register. As shown in Fig. 5, this circuit has a simple graphical structure.

Phase estimation allows one to determine the eigenvalue of a unitary operator, given the ability to prepare an eigenvector. If the application of the unitary, e.g., through Hamiltonian simulation, is inexpensive, then phase estimation is also inexpensive. Assume that $|ψα⟩$ is an eigenvector of the unitary operator U, so that $U|ψα⟩=e−iα|ψα⟩$⁠. Then repeated applications of the unitary yields information about the less significant digits of α, i.e., $U2k=e−2kiα|ψα⟩$⁠. This information can be stored in an ancillary quantum memory register by performing the operations

Now performing the QFT on the ancillary register yields the binary digits of the eigenphase

The circuit diagram for quantum phase estimation (QPE) is shown in Fig. 6.

Thus, if the eigenvector can be prepared efficiently and the unitary can be simulated efficiently, then the QFT can be used to extract the eigenphase efficiently. Note that there is nothing intrinsically quantum about this algorithm other than these assumptions. In fact, Kitaev et al.⁵ presented a version of this algorithm that simply uses direct measurements of the phase $cos (2kα)$ and $sin (2kα)$ dispenses with the ancillary qubit register altogether.

Phase estimation for certain classes of generalized eigenvalue problems, of the form commonly found in Stürm–Liouville problems, magnetohydrodynamics (MHD), and kinetic theory was developed in Ref. 134.

Understanding the key components of Grover's search algorithm led to the discovery of amplitude amplification and amplitude estimation. While the utility of Grover's search is often debated because the cost of data entry (for an exponentially large phone book) must be amortized over an exponentially large number of searches, the utility of amplitude amplification and estimation is not in question. Amplitude amplification can be used to amplify the overlap between two wavefunctions quadratically faster than this can be performed classically. This is the key part of Grover's search that allows it to obtain the correct answer with high probability with quadratically fewer queries that the classical counterpart.

Assume that one has a function oracle operator that marks chosen states by flipping their sign, i.e.,

where $f(x)$ is a binary function that returns the values 0 or 1. In particular, the oracle that marks a single particular state $|y⟩$ will be denoted $O|y⟩$ if $f(x)=−δx,y$⁠. (Note the minus sign.)

The idea behind amplitude amplification is to generate a rotation that amplifies the amplitude of the subspace of interest, called the good subspace. Now, the product of two reflections is a rotation around the point of intersection of the two lines of reflection. The rotation angle is twice the size of the angle between the two lines of reflection. Clearly, the oracle is a reflection operator, so, if one can find another suitable reflection, then it might be possible to generate a useful rotation. In fact, there are an infinite number of suitable choices.

Let S denote a unitary transformation that creates the state $|s⟩=S|0⟩$ from the vacuum state $|0⟩$⁠. This state should be a superposition of all states to be searched over; for example, the choice $S=H$⁠, which generates the uniform superposition state $|s⟩=H|0⟩$⁠, is often a useful choice. The amplitudes and phases of the superposition may be chosen almost arbitrarily, but only amplitudes that are not exponentially small will have a reasonable probability of being searched. The search oracle operator based on this state is

where $O|0⟩=2|0⟩⟨0|−1$ marks the vacuum state $|0⟩$⁠. The generalized Grover walk operator, $G$⁠, alternates between the search oracle and the function oracle

Now assume one is given an initial state $|ψ0⟩$⁠. Clearly, this can be decomposed into the sum of two parts: a part in the marked subspace, denoted $|good⟩$⁠, and a part in the unmarked subspace, denoted $|bad⟩$⁠, i.e.,

The action of the generalized Grover walk operator is to rotate the amplitude of the good and bad states via

The amplitude amplification process is illustrated in Fig. 7.

For example, for Grover's search, choose the walk operator to be defined by the uniform superposition state generated by $S=H$ acting on the vacuum. Assume that one is guaranteed that the function oracle $Of$ marks M out of N states. If one begins with an initial uniform superposition, then $sin (θ/2)=M/N$⁠. By choosing the appropriate value of k, i.e.,

one will generate a state that has the probability of measuring good states amplified to a high level. Then, simply performing measurements will yield these marked states with high probability. The final limiting form occurs in the limit that the dimension of the good subspace is much smaller than the overall dimension $M/N→0$⁠.

Once amplified, the less significant digits of the amplitude can be measured efficiently. Hence, the amplitude can be measured with quadratically fewer steps than it can by directly sampling the probability.¹⁰⁴ In fact, whenever the answer that one seeks in encoded in a probability amplitude, then the only way to extract this answer more efficiently than by direct sampling is by using amplitude estimation. Algorithms for quantum counting,¹⁰⁴ which can be used to compute sums and integrals up to quadratically faster than classical probabilistic algorithms^37,64 fall into this class. Clearly, these algorithms are essential for many scientific computing applications.

Let the unitary S generate the state of interest $|s⟩=S|0⟩$ whose overlap with the good subspace, marked by $Of$⁠, is to be measured. Assuming that both S and $S†$ can be performed efficiently, the search oracle $O|s⟩$ and the Grover walk operator, G, can be constructed efficiently.

The essential idea of amplitude estimation is to determine the least significant bits of θ through the action of powers of the Grover walk operator $Gk$⁠. For example, using phase estimation on an eigenvector of G would allow one to directly readout the value of the amplitude, $cos (θ)$⁠. In this case, we cannot construct an exact eigenvector, but we can construct the initial state $|ψ0⟩=|s⟩$⁠, which is a sum of two eigenvectors, as seen from Eqs. (54) and (55). The result of applying phase estimation to $|s⟩$⁠, as shown in Fig. 8, yields an equal superposition of the states $|θ⟩$ and $|π−θ⟩$ in the ancillary register.¹⁰⁴ Clearly, both states yield the same value for the amplitude.

Quantum walks on graphs are a vast generalization of the amplitude amplification paradigm. In its discrete-time version, quantum walks allow state changes based upon the outcome of flipping a quantum coin¹¹³ (or by rolling quantum dice). Quantum walks accelerate the solution of various problems, e.g., the search hitting time, typically quadratically. A heuristic explanation is that the ability to explore superposition states increases the effective speed of the walker.

If one includes Grover's search, amplitude amplification, and estimation within this group, then the majority of quantum algorithms that achieve a quadratic speedup today use quantum walks in an essential manner.¹³¹ Quantum walks have been used to speed up search on graphs,^135,136 accelerate simulation of Markov chains,¹³⁷ and Chebeyshev cooling schedules for quantum annealing,^138–140 as well as to solve the element distinctness problem.¹⁴¹ Notably, they are useful for Hamiltonian simulation^116,142 as well as for quantum signal processing.²² 

Although quantum walks are deterministic, they share similarities with classical random walks which are used as the basis for classical probabilistic Monte Carlo algorithms that simulate the evolution of a classical PDF in time. Due to the Feynman–Kac formula, this evolution can often be ascribed to Langevin dynamics defined by a stochastic forcing of a deterministic system, i.e., a stochastic differential equation (SDE). This temporal evolution of the PDF can be used as the basis of search algorithms that seek to find marked states within a certain hitting time. Another important goal is to evolve the PDF toward a desired statistical distribution which may physically represent a canonical ensemble, such as the Boltzmann distribution, the ground state of a wavefunction, or the solution to an optimization problem. After a few mixing times, the desired PDF is achieved with high accuracy and one can generate samples from the PDF and/or calculate the partition function. As the temperature is reduced, the PDF corresponds more and more closely to the ground state, which leads to the simulated annealing approach to optimization and ground state preparation.

Szegedy introduced¹³⁷ a quantum walk analog for any Markov process defined by a stochastic transition matrix. Like Grover's walk, it is composed of two reflections, but now on the tensor product of two Hilbert spaces. He proved that the hitting and mixing times of the quantum walk are generically improved quadratically because the spectral gap for the quantum walk is the square root of that for the random walk. This quadratically improves the complexity of using simulated annealing to find the solution to combinatorial optimization problems in addition to the quadratic speedup in accuracy obtained through amplitude estimation.^138–140 The improvement in convergence also improves the Chebeyshev cooling schedule.³⁷ The complexity of general multi-level Monte Carlo methods for SDEs can also generically be improved quadratically.³⁹ Quantum algorithms based on these concepts have recently been explored for simulating turbulent mixing and for computing the reaction rate within turbulent flows.^40,41

Simulation of many-body quantum systems was the inspiration for the development of quantum algorithms for quantum computers. The first efficient algorithms proposed were based on the realization that many Hamiltonians of interest have a sparse decomposition into a sum of non-commuting parts, where each part is composed of commuting terms, and, hence, each part can be computed individually as a product of terms.¹⁰⁸ Then, one can use operator splitting techniques, such as the Trotter–Suzuki decomposition and Lie group decomposition formulas,¹⁴³ to approximate the unitary corresponding to the full Hamiltonian as a product of terms with sufficiently small time steps. In fact, while these techniques are well known to the operator splitting community, a rigorous formulation of the approximation theory was only recently developed.¹¹² 

For example, for the Schrödinger equation for particles interacting with few-body potentials, the Hamiltonian breaks into a sum of the kinetic energy of all particles and the potential energy of all particles. One can then efficiently switch back and forth between the momentum and position basis in order to evaluate each as a simple phase shift.^110,144 There is also a similar decomposition for lattice Hamiltonians with nearest neighbor coupling. These techniques can also be applied to both bosonic and fermionic systems as well as to systems of distinguishable particles.¹⁰⁹ Lattice simulations may be directly applied to model relativistic quantum plasmas.^145,146

Quantum algorithms for general sparse Hamiltonians were then developed by a number of authors.^115,147 Today, optimal strategies are based on Childs' generalization of Szegedy's approach to define a quantum walk for an arbitrary Hamiltonian.^116,142 Using strategies such as the linear combination of unitaries¹¹⁸ and Taylor series approximations,¹²⁰ the algorithms were improved to have nearly optimal dependence on all key parameters, including condition number, precision, and failure probability, in Refs. 117, 119, and 120 and, eventually, to optimal dependence in Refs. 21 and 121. Optimal dependence on all parameters was actually first achieved for a large class of Hamiltonians using the methods of qubitization and quantum signal processing (QSP).^22,129 This method works by using QSP to correct the eigenvalue distribution of the Childs walk so that it generates a more accurate approximation of the underlying Hamiltonian.

Harrow, Hassidim, and Lloyd (HHL)¹²³ invented the first general quantum algorithm for solving sparse linear equations, $Ax=b$⁠. The HHL algorithm estimates the eigenvalues of A and replaces them by their inverse using an ancillary qubit as well as a technique called post-selection, i.e., measurement of the state of the ancillary qubit. The method was subsequently improved using the methods of variable time amplitude amplification¹²⁴ and linear combinations of unitaries.¹²⁵ General matrix operations can be performed using the QSP and qubitization techniques.¹²² 

After the introduction of linear solvers and advanced Hamiltonian simulation techniques, methods for solving linear differential equations were developed.^18,128 For hyperbolic partial differential equations, such as general wave equations, Ref. 19 found similar results. More recent linear differential equation solvers have much improved dependence on precision by using a linear combination of operations¹⁸ and by using spectral methods.^20,148

The quantum linear solver algorithm (QLSA) and quantum linear differential equation algorithm (QLDSA) have an exponential improvement in the complexity of solving for a quantum encoded version of the solution. Note, however, if the amplitudes of the entire solution are required, then the speedup is reduced to quadratic. For example, Ref. 149 found that, after properly accounting for precision, the improvement for solving elliptic PDEs using finite elements is only quadratic.

The QM Postulates require linear evolution of the state, and, if the evolution is ideal, the linear operation must be unitary. Hence, nonlinear operations of the entire state cannot be performed with a speedup over classical algorithms. As discussed previously, nonlinear quantum eigenvalue transformations (QEVT) of a normal matrix¹²⁹ and nonlinear quantum singular value transformations (QSVT)¹²⁶ can be performed efficiently for a block-encoded matrix. This approach has begun to yield efficient methods for performing nonlinear operations on blocks of the wavefunction.^150,151 An alternate approach, based on considering a nonlinear operation on a given set to be a linear operation on the Hilbert space of functions on the set was explored by Ref. 152. This method is the analog of the Koopman–von Neumann approach, discussed below, for finite time intervals.

There is good potential for these methods to develop into the ability to efficiently perform arbitrary nonlinear operations. If this is the case, then the ability to perform a nonlinear operation, combined with the ability solve linear equations, allows one to develop methods for solving nonlinear equations using either fixed-point Picard iteration or Newton's method.

The first concrete quantum algorithm to simulate nonlinear differential equations was proposed in Ref. 28, but had an exponential complexity in time. This is because the entire wavefunction needed to be stored for each time step in order to compute the required nonlinear functions. After this result, the quantum algorithms community began to focus on solvers for linear differential equations, as described in Sec. III E 7.

The next conceptual breakthrough was provided by Ref. 29, which proposed an algorithmic approach to solve nonlinear differential equations by mapping the nonlinear problem to the linear Koopman–von Neumann representation of classical dynamics. After this step, Hamiltonian simulation could be directly applied to the system. Then, Liu et al.³¹ derived an explicit algorithm for differential equations with quadratic nonlinearities based on the closely related technique of Carleman linearization. The connections between these techniques will be elucidated in detail in Sec. V.

Another algorithmic approach for solving the dissipative nonlinear differential equations that arise from many-body physics was also recently advocated by Ref. 153. The basic idea, which goes back to initial considerations of quantum computing,^154,155 is that such a description is an approximation of many-body quantum physics. Hence, a simulation of the many-body quantum physics Hamiltonian should lead to an approximate solution of the classical system in the appropriate limit.

Yet another approach for solving the Navier–Stokes equations of fluid dynamics was explored by Gaitan¹⁵⁶ based on a method for integrating differential equations developed by Kacewicz.³⁸ In turn, the time-integration method developed by Kacewicz relies on quantum algorithms developed by Heinrich and Novak for computing sums and integrals.^37,64 As explained previously, these methods can obtain up to a quadratic quantum speedup over Monte Carlo methods for high-dimensional systems with non-smooth right hand sides. However, a key point is that this speedup occurs when only partial information is given to specify the right hand side,¹⁵⁷ i.e., only a certain number of derivatives exist and the position of possible discontinuities in the last derivative are not known in advance. Thus, such gains can only achieved if there are non-smooth initial conditions whose evolution is too difficult to track in time or if there is a stochastic forcing of the differential equations. Stochastic forcing is in fact useful for studies of fluid and plasma turbulence as well as for stochastic differential 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 with analytic coefficients and initial conditions.

Variational quantum algorithms^{133,158–160} are typically formulated as hybrid quantum-classical algorithms for solving optimization problems. All of the quantum subroutines discussed above require large-scale fault-tolerant quantum computers which will not be available in the near term. As a middle ground for the NISQ-era, hybrid quantum-classical algorithms have been proposed where quantum devices are envisioned as co-processors that accelerate the otherwise classical calculations. It is hoped that these will serve as NISQ algorithms for difficult many-body problems in quantum chemistry and quantum materials science, typically involving fermions. They may even be useful as methods for solving nonlinear problems.⁴⁴ For example, variational simulation of the Navier–Stokes equations has recently been nvestigated.^45–47 Variational simulation of nonlinear stochastic differential equations has also been explored.¹⁶¹ 

The idea is that there is a cost function of high complexity that is evaluated by a quantum computer which depends on a number of parameters that is minimized by using a classical optimization algorithm. For example, the parameters may be in the form of phase angles associated with the circuits that are to be optimized to minimize the cost of the objective function. Quantum resources are useful because they allow the use of Hamiltonians with many degrees of freedom and of ansatzes that are classically expensive. Thus, they have a few key independent parts: (i) constructing the cost function via Hamiltonian simulation, (ii) constructing a suitable ansatz for the initial wavefunction, (iii) developing efficient measurement protocols to determine the cost function and, potentially, its derivatives with respect to optimization parameters, and (iv) the classical optimization algorithm.

While it is often hoped that these algorithms will be useful for achieving at least quadratic quantum advantage on near-term NISQ hardware platforms, there are a number of outstanding challenges that must be considered and these challenges can make it difficult to obtain rigorous proofs of their complexity. Notable issues are associated with developing efficient ansatzes for the initial wavefunction, developing efficient measurement techniques, avoiding barren plateaus with little information on the optimization problem, and the presence of many local minima, which make the problem intrinsically difficult. In fact, the classical-quantum hybrid scheme may not have any advantage when the classical optimization problem is NP-hard.¹⁶² 

The quantum simulation application space is closest to the typical applications used within FES. Several authors^{17,29,31,32,163} have outlined a framework for the quantum simulation of nonlinear dynamics. As discussed in Sec. III, there are a number of subroutines to be employed that have steadily improved in computational complexity and accuracy. A summary of the key steps is discussed in this section.

The key conceptual steps are as follows:

1. Determine an embedding of the dynamics into a linear differential system of equations.

2. Numerically discretize the linear system.

If the system is intrinsically linear, like a quantum Hamiltonian, then the first step is relatively straightforward, but the choice of representation can still have an impact on complexity, which is especially important when resources are limited. If the system is nonlinear, then a linearization must be determined. In either case, it is important to choose the basis so that the resulting equations have a sparse structure or an efficient Trotterization that can be utilized by the quantum simulation algorithms.

The key computational steps are as follows:

1. Prepare the initial state.

2. Evolve the system in time using an efficient linear differential equation solver, or if the system is Hamiltonian, an efficient Hamiltonian simulation algorithm.

3. Estimate the observable to be measured efficiently, typically by using an ancillary quantum memory register to encode the observable and performing amplitude estimation.

If the linearized system is sparse, then using the sparse linear differential equation solver ensures that the key steps are efficient.³¹ If the linearized system is both sparse and Hamiltonian, then using the sparse Hamiltonian simulation algorithm ensures that the simulation step is efficient.²⁹ State preparation is required for the initialization step and is often required for the final measurement stage. For an efficient algorithm, one must choose initial conditions that can be prepared using an efficient method such as Hamiltonian simulation or a quantum walk.²⁹ 

Many works on quantum computing algorithms consider the output of a subroutine to be a quantum encoding of the solution, i.e., wavefunction or density matrix that can then be used as input for another subroutine. While this is quite reasonable for intermediate stages of a quantum program, for almost all scientific calculations of interest, the output usually consists a set of real floating point numbers that quantify the system of interest. Whenever these data are encoded in the amplitude of the wavefunction, then one requires a final stage where amplitude amplification and estimation is used to extract the quantities of interest.

Choosing a finite representation of the wavefunction is clearly required for numerical simulation. Choosing a basis that keeps the operations sparse is clearly an important consideration. In general, a complete basis allows one to represent the wavefunction, $ψ(x)$⁠, as a linear combination of basic elements, ${ϕk(x)}$ so that

The variables x are indices representing the computational basis. If the computational basis represents position, then it might be convenient to use the position basis $ϕk(x)=δkx$ or the momentum basis, $ϕk(x)=exp (2πikx)/N1/2$⁠. More generally, the discretized Hilbert space could span any finite collection of linearly independent functions, such as orthogonal polynomials, Bessel functions, etc. One could also use finite difference or finite element basis functions, i.e., orthogonal polynomials (or, more generally, orthogonal functions) over subdomains with finite support. If the differential equation to be solved has an efficiently computable sparse representation in terms of these basis functions, i.e., s terms, then the resulting system of equations will be s-sparse.

One can also use an overcomplete basis, such as the coherent states. Coherent states are an important example of a reproducing kernel Hilbert space (RKHS), which is defined by the fact that the pointwise evaluation of functions is continuous.^164,165 In this case, the Riesz representation theorem implies that there is a continuous linear kernel function, $Ky(z)$⁠, that allows one to evaluate a function at the point y via

This is called the reproducing property of the kernel. The definition of the bilinear kernel function in terms of the Hilbert space inner product, $K(y,z)=⟨Ky|Kz⟩$⁠, demonstrates that the kernel function must be symmetric and positive definite. The Moore–Aronszajn theorem states that every symmetric positive definite kernel defines an RKHS.

There are actually many types of RKHS' that are used within the machine learning community for efficient representations of data in higher-dimensional spaces. One may expect that RKHS' should play a similar role in developing efficient representations for quantum computing and should certainly be of interest for quantum machine learning (QML). Examples of RKHS' are bandwidth-limited functions, finite collections of orthogonal polynomials, and complex analytic function spaces, such as Segal–Bargmann coherent states, Bergman spaces, and Hardy spaces. For example, as proven in  Appendix C, the Carleman representation is based on the standard Hardy space, a complex analytic RKHS that uses a Szëgo kernel. This RKHS is based on the collection of monomials, $zk$⁠, where evaluation is given by the Szëgo kernel for the boundary of the unit disk in the complex plane. Determination of which representations are more or less efficient for representing different applications is an active research area.

The Stone–von Neumann theorem implies that all representations that satisfy the (Weyl exponential form of the) canonical commutation relations (CCR) are unitarily equivalent.¹⁶⁶ While there are important exceptions, the vast majority of the representations used in practice are unitarily equivalent to the position and momentum operators; i.e., there is a unitary transformation that relates the two different representations. For example, the number basis defined by normalized Hermite polynomials, $Hk(x)$ are related to the position representation by a unitary transformation. Similarly, the coherent states are related to the number states via a unitary transformation. The combination of these two transformations generates another unitary transformation from position space to coherent state space known as the Segal–Bargmann transform.

There have already been some interesting uses of RKHS' within quantum computing for nonlinear dynamical systems. Reference 163 used a real RKHS that used smoothed basis functions to develop a well-defined Hamiltonian simulation algorithm for integrable systems that displays smooth pointwise convergence in the limit $N→∞$⁠. Reference 31 used the Carleman embedding, which is based on the standard Hardy space, in developing a quantum algorithm for dissipative systems, such as the logistic equation and the Burgers equation. Reference 32 considered the possibility of exploring other representations of the CCR. As proven in  Appendix C, all of the concrete examples discussed in Ref. 32 are well-known examples of complex analytic RKHS'.

In the typical models of quantum computing, the initial wavefunction must start in a standard state, such as the ground state $|0⟩$⁠. This is useful from the point of view of a physical device, because, if the temperature is much lower than the energy difference between the first excited state and the ground state, then there is an exponentially small likelihood of finding the initial wavefunction in an excited state. Due to dissipative decoherence processes such as relaxation, simply waiting long enough allows one to re-initialize the quantum program from $|0⟩$⁠. Active restart techniques intentionally use decoherence to speed up the process of resetting to vacuum.

Then, to actively begin the program, one must typically initialize the wavefunction in a specific state. For simple programs, such as search, one might use the Walsh–Hadamard transform to prepare a uniform superposition of a quantum memory register. For a simulation, one would typically initialize the memory register with a state corresponding to the initial condition of interest. As discussed previously, this initial condition should be efficiently computable, perhaps a smooth function with some random noise which might be useful for triggering the dynamics of interest, e.g., as in simulations of turbulence. As discussed previously, many useful functions can be computed efficiently by using quantum walks or Hamiltonian simulation.

In order to obtain output from a given memory register, a projection step is usually required to obtain the data of interest. In practical models of quantum computing, there is typically only a possibility of performing measurements in the computational basis. Hence, in order to extract other properties of the wavefunction, e.g., measurements along another axis of the Bloch sphere, one must perform judiciously chosen unitary transformations before the measurement is taken. This final projection step can either be thought of as a transformation of the data by some unitary operation U or as a transformation of the state to be measured against with the inverse operation $U†$⁠.

After defining an efficient representation, the simulation step has now reduced the problem to solving for the evolution of a linear set of ODEs. In components, these ODEs can be written as

If the evolution is unitary, then one can directly use the best available sparse Hamiltonian simulation algorithm (HSA).^108,109,117 This constructs an approximation of the form

where $t$ is the time ordering operator. The result will typically be approximate due to the need to simulate the Hamiltonian efficiently. If the Hamiltonian can be broken into a small number of pieces where 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 decomposition to generate an accurate approximation. Otherwise, one can choose one of the sparse Hamiltonian solvers described previously.

If the evolution is not unitary, then one can use the best available quantum linear differential equation solver algorithm (QLDSA).^18,128 Simple statements of these algorithms might often an explicit forward Euler time step; however, this can clearly be improved to use higher order methods as well as implicit timestepping, using the QLSA. For example, for each time step, one can consider performing a time splitting of the form

Here, θ represents a time splitting parameter that determines the amount implicitness to be used; i.e., θ = 0 for fully explicit forward Euler, $θ=1/2$ for implicit midpoint Crank–Nicolson, and θ = 1 for fully implicit backward Euler.

In order to handle a non-unitary evolution, the QLDSA must use post-selection. Care must be taken to be able to extract the information of interest and this depends on the properties of the solution itself. If the linear system is unstable and the magnitude of ψ is exponentially unbounded, then the output will likely result in noise. If the linear system is stable and the magnitude of ψ shrinks exponentially toward zero, then the output will also eventually result in noise.

For realistic applications of interest, the no-fast-forwarding theorem holds.^117,147 The no-fast forwarding theorem implies that in the generic situation where the solution cannot be solved for semi-analytically, the solution must then be tracked on the shortest relevant time steps. This implies that the simulation will be restricted to short time steps, of the order of $δt∼ℏ/||H||$⁠. Hence, the total number of time steps, $Nt=Δt/δt∼||H||Δt/ℏ$⁠, will typically be large. For example, to observe temporal variations on the longest relevant timescale would require $Δt∼ℏ/min(eigenvalues(H))$⁠, so the number of steps will typically scale as the condition number $Nt∼κ(H)$⁠. Hence, while the QLDSA is more costly than Hamiltonian simulation, particularly with regard to condition number, the QLDSA is similar to that of Hamiltonian simulation in terms of complexity in time steps. Thus, which choice is optimal may depend on other desired features of the solution.

In fact, there are some relatively simple methods for extracting observables from a set of measurements.²⁹ To simplify the discussion, here we consider the case where all observables commute with each other (are in involution). The expectation value of a local observable

can be found by performing a reversible classical computation of the functions

First, form the states

Then, to compute the amplitude, add an ancillary qubit to the system, $|ϕ⟩⊗|0⟩$⁠, and perform a specific choice of rotation, $Rϕ$⁠, defined via

Finally, using the amplitude estimation subroutine, discussed in Sec. III E 4, allows one to directly calculate the amplitude of the $|1⟩$ state to determine $sin2(θ)=⟨O⟩/N$ with an amount of work that scales as 1/accuracy rather than 1/accuracy².

For both the intrinsically quantum and the intrinsically classical approaches, quantum Hamiltonian simulation of classical dynamics leads to large gains in memory and computational complexity. The Hamiltonian simulation algorithms are efficient as long as the Hamiltonian is sparse. It is a simple exercise to prove that the Hamiltonians of interest are sparse for N-body interactions and for particles interacting with electromagnetic fields. This implies up to exponential savings relative to an Eulerian discretization of the Liouville equation.

However, the best classical algorithms of this sort are time-dependent Monte Carlo methods.¹⁶⁷ For example, particle-in-cell algorithms invented in the plasma physics community fall into this class, as well as molecular dynamics algorithms invented in the condensed matter community. These probabilistic classical algorithms can also provide an exponential speedup over Eulerian methods when the same restriction on the observables of interest are applied. In this case, one does not need to accurately sample the entire PDF, rather, one only needs to maintain good accuracy for the most probable regions of phase space.¹⁶⁷ In this case, the equivalent quantum algorithm can achieve up to a quadratic speedup over time-dependent Monte Carlo (e.g., PIC) methods.^168,169

It is important to stress that it is the use of amplitude estimation of physical observables and the use of quantum walks for state preparation and measurement transformations that enable the quadratic speedup. Due to the central limit theorem, when classically sampling an observable with a finite mean and variance, the relative accuracy decreases as $ε∼1/S1/2$ where S is the number of samples. This implies that the amount of computational work required to achieve a certain accuracy grows as $1/ε2$⁠. The utility of the amplitude estimation algorithm is that it only requires an amount of computational work that grows as $1/ε$⁠.

The majority of simulations used within FES today aim to solve initial-value and boundary-value problems: typically time-dependent simulations of nonlinear dynamics. While the system may often represent partial differential equations (PDEs), such as the Navier–Stokes or magnetohydrodynamics (MHD) equations of fluid mechanics, or the kinetic equation for the probability distribution function, ultimately, numerical discretization of the set of PDEs typically reduces the system to time integration of a set of ODEs. The question then arises: How should quantum computers be used to efficiently simulate nonlinear dynamics?

The past few years have seen intensive research activity on the ability to perform quantum simulations of nonlinear dynamics. As described above, these approaches generally rely on embedding the nonlinear system into a large, in fact, infinite-dimensional linear system of equations.^17,29,32,33 In this section, we simultaneously unify these approaches and systematically categorize their differences.

Classical dynamics is distinguished by the fact that the solution to a well-posed set of nonlinear differential equations typically has both existence and uniqueness theorems that apply to almost all points in phase space, aside from a set of measure zero. The uniqueness of the individual trajectories then allows one to associate a unique evolution of other objects, such as scalars, vectors, and tensors on phase space. In particular, the Liouville equation expresses the fact that a volume form, which can be interpreted as a classical probability density, is advected by the flow in a specific manner that ensures that probability is conserved in time. The action of the evolution on phase space functions presents a linear but infinite-dimensional system of equations. Due to the Stone–von Neumann theorem, the vast majority of all possible representations are equivalent to this one, simply defined by translations in phase space.

From an algorithmic point of view, the next question that arises is: How well does a finite truncation of the linear representation approximate the underlying nonlinear dynamics?

If one begins with a Hamiltonian system, then a natural approach is to use a quantized version of the classical physics model. There are many potential advantages to this approach. First of all, the quantized version may be the more accurate physical model, so important effects such as interference, tunneling, and diffraction would be correctly described. Second of all, there are already many quantum algorithms for the simulation of sparse Hamiltonians, and, this is especially true for quantum Hamiltonians.

Two examples of this approach, as well as their implementation on present-day quantum hardware, will be discussed in detail in Sec. VI. In particular, Sec. VI B explores the physics of the sawtooth map,¹⁷⁰ which can be viewed as a poor-man's model for wave–particle interactions that uses minimal numerical resources. Figure 9 compares the dynamics of the classical sawtooth map to that of the quantum sawtooth map as the number of qubits is increased from 6 to 9 to 16. To visualize the quantum dynamics in phase space, the Husimi-Q quasi-probability distribution is used to visualize the density matrix. In both cases, the initial state has uniform amplitude along the $p/π=3/4$ momentum surface, and, as time evolves, both the quantum wavefunction and the classical PDF explore similar regions of phase space. However, this comparison only works well once the number of quantum states becomes large enough. At n = 9 and $N=29=512$ states, for a total phase space resolution of $512×512≃2.6×104$⁠, the quantum Husimi-Q distribution still appears quite “fuzzy” to the human eye. At n = 16 and $N=216≃6.6×104$⁠, for a total phase space resolution of $232≃4.3×109$⁠, quantum artifacts are no longer visible at this scale.