{'fraction_non_alphanumeric': 0.0305092976986825, 'fraction_numerical': 0.009947832625757836, 'mean_word_length': 5.004750482103382, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum algorithms for quantum chemistry and quantum materials science

As we begin to reach the limits of classical computing, quantum computing has emerged as a technology that has captured the imagination of the scientific world. While for many years, the ability to execute quantum algorithms was only a theoretical possibility, recent advances in hardware mean that quantum computing devices now exist that can carry out quantum computation on a limited scale. Thus it is now a real possibility, and of central importance at this time, to assess the potential impact of quantum computers on real problems of interest. One of the earliest and most compelling applications for quantum computers is Feynman's idea of simulating quantum systems with many degrees of freedom. Such systems are found across chemistry, physics, and materials science. The particular way in which quantum computing extends classical computing means that one cannot expect arbitrary simulations to be sped up by a quantum computer, thus one must carefully identify areas where quantum advantage may be achieved. In this review, we briefly describe central problems in chemistry and materials science, in areas of electronic structure, quantum statistical mechanics, and quantum dynamics, that are of potential interest for solution on a quantum computer. We then take a detailed snapshot of current progress in quantum algorithms for ground-state, dynamics, and thermal state simulation, and analyze their strengths and weaknesses for future developments.


I. INTRODUCTION
In this section, we briefly introduce the concept of quantum computer and current quantum computing architectures.

A. What is a quantum computer? How is it relevant to quantum simulation?
A quantum computer is a device which computes by manipulating quantum information [1][2][3][4][5][6]. The basic unit of quantum information, called a qubit, is synonymous with a two-level quantum system. Denoting the two basis states of a qubit as |0 , |1 , the general single-qubit state may be a superposition c 0 |0 + c 1 |1 where c 0 , c 1 are complex numbers. For n qubits, there are 2 n basis states, which can be enumerated as the bitstrings |00 . . . 0 , |00 . . . 1 , etc. To specify a general quantum state, one must specify a complex coefficient (amplitude) for each basis state, e.g. c 00...0 , c 00...1 etc. The exponential number of amplitudes needed to specify the state of n qubits stands in contrast to the linear amount of information needed to encode a single bitstring, the state of n classical bits.
Measuring the qubits (in this basis, often referred to as the computational basis), collapses the state onto the measurement outcome, which is a single (random) bitstring x that appears with the probability |c x | 2 . Despite this probabilistic interpretation, the fact that quantum amplitudes can be complex prevents one from treating quantum states as classical probability distributions, which is at the heart of the quantum phenomenon of entanglement, the correlations in a system that cannot be generated by a classical distribution of states. The possibility of creating entangled states in a space of exponentially large dimension and manipulating these states by exploiting their constructive or destructive interference is the source of the computational power of a quantum computer.
A convenient way to describe the manipulation of quantum states in the context of quantum computation is via few-qubit operations called gates. The action of each gate can be viewed as a unitary time evolution of the n-qubit system under a suitable Hamiltonian that acts non-trivially only on a few qubits (usually one or two). Any quantum computation can be expressed by a sequence of elementary gates, called a quantum circuit, applied to an initial basis state, e.g. |00 . . . 0 , and followed by the measurement of some set of the qubits.
The exponential separation between quantum and classical information does not simply mean that quantum computers can compute answers to problems "in parallel" with an exponential speedup, for example, by storing multiple different solutions in the many amplitudes. This is because reading out from a quantum computer destroys the state, and thus to harness the power of quantum information, a quantum algorithm must coordinate interference between the amplitudes such that useful information can be read out with high confidence without many measurements.
The interest in quantum computing for quantum simulations of molecules and materials stems from the fact that in many cases, the chemistry and physics of molecules and materials is best described using quantum mechanics. Thus, the state of a many-particle molecule also encodes quantum information, and as the number of atoms increases, similarly can require an exponentially large number of classical bits to describe. This means that in the worst case, quantum simulation is exponentially hard on classical computers. This is the motivation for Feynman's famous observation that "Nature isn't classical, dammit, and if you want to make a simulation of nature, you'd better make it quantum mechanical" [7].
A moment's reflection, however, suggests that the potential quantum advantage for a quantum computer in quantum simulation is nonetheless subtle. For example, if it were indeed impossible to say anything about how atoms, molecules, or materials behave, without using a quantum computer, there would be no disciplines of chemistry, condensed matter physics, or materials science! Decades of electronic structure and quantum chemistry simulations suggest that reasonable, and in some cases very high accuracy, solutions of quantum mechanics can be obtained by classical algorithms in practice. Quantum advantage in quantum simulation is thus problemspecific, and must be tied both to the types of questions that are studied, as well as the accuracy required.
FIG. 1. Schematic of a quantum simulation of quantum dynamics. A quantum simulator (bottom) is prepared in an initial state |ψ(0) , that is a representation of the initial state |φ(0) of the actual system of interest (top). The simulator is manipulated by a unitary transforma-tionÛ that is an encoding of the real-time evolution of the system of interest, and the final simulator state |ψ(t) is measured, yielding information about the dynamics of the original system. From Ref. [8].
We can look to theoretical quantum computer science to better understand the power of quantum computers in quantum simulation. The natural problem to solve on a quantum computer is the time evolution of a quantum system given some initial state, This problem is representative for the complexity class BQP, i.e. it is of polynomial cost on a quantum computer and believed to offer a clear separation from the classical case (an exponential quantum speedup has been rigorously proven only in the query complexity setting [9]). However, it is necessary to prepare the initial state, which may be difficult by itself. In particular, preparing a low-energy state may be challenging, which naturally leads to considering other important problems, Thermal averages: Â = Tr Â e −βĤ Tr e −βĤ Ground state determination lies in complexity class QMA [4], a class of problems not known to be efficiently solvable (in general) on a quantum computer. This also means that thermal averages cannot in general be computed efficiently on a quantum computer, since in the limit of zero temperature, this problem reduces to ground-state determination. Although it is the case that there are many physical ground state and thermal problems that are not so hard to solve in practice (as demonstrated by efficient classical simulation of many problems) and similarly many initial states of interest in quantum dynamics that are easy to prepare, the above tells us that in a rigorous sense, we do not have a straightforward theoretical guarantee that a quantum computational advantage can be naturally achieved for the central problems in quantum simulation. Given the limits to the guidance that can be provided by rigorous computational complexity results, it is clear that to understand quantum advantage in chemistry, condensed matter physics, and quantum materials science, we must be guided by actual empirical data in the form of numerical and theoretical experiments with quantum algorithms and quantum devices on simulation problems of interest. This requires making progress on both theoretical and practical questions of quantum simulations, ranging from the basic algorithms and choices of encoding and representation to issues of circuit compilation, readout, and mapping to specialized hardware. A central purpose of this review is to provide a perspective on what the relevant chemical and materials problems are today; to give a snapshot of the limitations of classical methods for these problems; and in these contexts to understand the strengths, weaknesses, and bottlenecks of existing ideas for quantum algorithms, and where they need to be improved, both in terms of fundamental theoretical aspects as well as practical methods of implementation.

B. Current quantum architectures
The idea of using a quantum mechanical device to perform a computation was first considered in earnest by Richard P. Feynman in a famous lecture in 1982 [7]. His suggestion was to build a lattice of spins with tunable interactions. He conjectured that by appropriately tuning these interactions, such a system could be made to imitate the behavior of any other (bosonic) quantum system with the same dimensionality, and thus could serve as a way to compute the properties of some other system that one would like to study. This idea, which is often referred to as analog quantum computation, is still very much alive today and embodied in the field of cold atomic gases and related techniques, which has made great progress in simulating complicated physics of strongly correlated systems in a controlled environment. A general schematic of the idea is shown in 1.
However, such systems come with some limitations: for example, it appears unlikely that a system of spins can naturally emulate a system of fermions, or that a two-dimensional array of spins could emulate a three-dimensional material; in other words, a given analog quantum computer can only simulate a certain set of physical systems. Furthermore, it is itself an experiment with limited accuracy; while it may be easier to perform than experiments on the original system of interest, the accuracy nonetheless remains limited. A more general approach, which is the focus of this review, is digital quantum computation. Here, very much in analogy to the classical computers we are used to, one considers a system of quantum registers -qubits -that are controlled through some set of instructions -the quantum gate operations. One can show that with a finite (and indeed relatively small) set of such gate operations, one can in principle generate arbitrary quantum evolution to arbitrary precision! In other words, every problem in the complexity class BQP can be mapped into a quantum program, i.e. the quantum circuit composed of such gates.
A third model of quantum computation is adiabatic quantum computation. Here, the computation is encoded into a time-dependent Hamiltonian, and the system is evolved slowly to track the instantaneous ground state of this Hamiltonian. It can be shown that this model is equivalent to circuitbased digital quantum computation [10], but it is usually considered to be less practical. However, a restricted version of it, adiabatic quantum optimization [11], has gained some popularity. Here, a classical optimization problem is encoded into a quantum Hamiltonian, to which one adds some noncommuting terms to endow the system with non-trivial quantum dynamics. If one then slowly turns off the quantum terms, the optimal solution to the classical problem should be obtained. In practice, however, one may have to go impractically slowly for this to be true; the detailed analysis of this approach is a complex problem that is not covered in this review.

C. Building a circuit-based digital quantum computer
Returning to circuit-based quantum simulation, what is the status of the field today? The natural enemy of quantum computation is decoherence of the qubit, i.e. the tendency of the stored quantum state to decay into a classical state. After decades of research, a number of qubit technologies [12,13], for example superconducting qubits [14] and ion traps [15], have reached the point where small devices of a few dozen qubits can be sufficiently isolated from decoherence to execute non-trivial quantum algorithms of a few tens to hundreds of gates. This has been termed the noisy intermediate-scale quantum (NISQ) era [16]. With these devices, the first demonstrations of quantum supremacy appear to be possible. For this, an artificial but well-defined problem that is intractable classically is solved on a quantum computer [17]. However, these problems do not generally have practical relevance, and a more important and challenging question is when a scientifically or commercially relevant problem can be solved on a quantum computer.
To address large-scale problems, it is necessary to correct errors that occur on physical qubits during the computation. While the no-cloning theorem prevents error correction by simple redundancy, it can be shown that quantum error correction is possible nonetheless by encoding a single qubit into an entangled state of many qubits [19], to use suitable measurement to detect errors occurring in the system and to apply FIG. 2. Examples of error correction schemes. (a) Quantum circuit implementing the bit-flip code. The state of a single logical qubit, |ψ = α|0 + β|1 is encoded in the state of three physical qubits as α|000 + β|111 . The system is then sent through a channel E bit = E ⊗3 , where E flips a qubit with probability p (X error). Measuring the Z operator on the last two qubits permits to determine whether a qubit has been flipped by the channel, and to correct such an error without corrupting the transmitted state. (b) Schematic of the surface code for quantum error correction and quantum fault-tolerance. Top: a two-dimensional array implementation of the surface code, data qubits are open circles and measurement qubits are filled circles, with measure-Z qubits colored green (dark) and measure-X qubits colored orange (light). Middle: sequence of operations (left), and quantum circuit (right) for one surface code cycle for a measure-Z qubit, to detect sign flip errors. Bottom: similar, but for measure-X qubits. Adapted from Ref. [18]. suitable unitary transformations correcting such errors without disturbing the information encoded in the system (see 2). This leads to the very important distinction between physical and logical qubits. The latter are error-corrected and encoded in the state of many physical qubits. Quantum algorithms are performed on the logical qubits, and the error-correction scheme translates the operations on logical qubits into physical operations. This incurs significant overhead: depending on the error rate of the underlying physical qubits and the target error rate of the logical qubits, one may need hundreds or even thousands of physical qubits to realize a single logical qubit. Therefore, when evaluating the capabilities of some qubit platform with respect to an algorithm, one must be careful to include the cost of encoding logical qubits into physical qubits to achieve the required error rates.
FIG. 3. Single-reference (left) and multi-reference (right) wavefunctions. The former is qualitatively described by a single Slater determinant, the latter by a linear combination of a potentially large number of Slater determinants. Often such determinants correspond to different configurations of electrons in an "active space" of orbitals.

II. SIMULATION CHALLENGES IN MOLECULAR AND MATERIALS SCIENCE
In this section we describe a representative, but certainly non-exhaustive, set of scientific problems of relevance to quantum simulation.

A. Quantum chemistry
Quantum chemistry is concerned with determining the lowlying eigenstates of the electronic Hamiltonian of a molecule. The eigenstates are determined for fixed sets of nuclear positions, i.e. within the Born-Oppenheimer approximation. Determining the main features of the resulting potential energy surface, i.e. the electronic energy as a function of nuclear positions, its minima and saddle points, is key to understanding chemical reactivity, product distributions, and reaction rates.
There exists a wide range of quantum chemical methods with different accuracy and speed tradeoffs, ranging from density functional methods that routinely treat thousands of atoms or more [20], to high-level many-electron wavefunction methods, such as coupled cluster theory, that can attain chemical accuracy of 1 kcal/mol and better, on systems of tens of atoms [21,22]. However, most methods in quantum chemistry are most accurate for problems where there is a dominant electronic configuration, a subset of the quantum chemistry problem known as the single-reference problem. Singlereference quantum chemistry is found in the ground-states of many simple molecules (e.g. hydrocarbons), but in many molecular excited states, at stretched bond geometries, and in transition metal chemistry, multiple electronic configurations can come into play, which is referred to as multi-reference quantum chemistry. Despite much work (and progress) in extending quantum chemistry to multi-reference situations, the attainable accuracy in molecules with more than a few atoms is significantly lower than in the single-reference case. Some examples of multi-reference quantum chemical problems in-clude: • The chemistry of enzyme active sites. Such active sites can involve multiple coupled transition metals, famous examples being the 4 manganese ions in the oxygen evolving complex [25], or the 8 transition metals in the iron-sulfur clusters of nitrogenase [26], shown in 4. They present some of the most complicated multireference quantum chemistry problems in the biological world. Combined theoretical and experimental studies, primarily at the level of density functional theory, have proven successful in unravelling many structural and electronic features of such enzyme active sites [27][28][29]. However, a detailed understanding of the interplay between spin-coupling and delocalization between metals, which requires true multi-reference quantum chemistry and is needed to interpret aspects of experimental spectroscopy, remains in its infancy [23,24,[30][31][32].
• Transition metal nanocatalysts and surface catalysts.
Similarly to enzyme active sites, simulating the mechanism of action of synthetic heterogeneous catalysts remains a major challenge. While density functional theory has been widely employed, predictions of even basic quantities such as the adsorption energy of small molecules are unreliable [33,34]. While not all such problems are expected to be multi-reference in character, even the single-reference modeling of such chemistry, at a level significantly beyond density functional theory, is currently challenging or impossible. In addition, multi-reference effects are expected to play a role in certain catalysts, such as transition metal oxides, or at intermediate geometries in reaction pathways [33,[35][36][37].
• Light harvesting and the vision process. The photochemistry of conjugated organic molecules is the means by which Nature interacts with light. Some prominent examples of such natural conjugated systems include the carotenoid and chlorophyll pigments in the light-harvesting complex of plants [38,39], as well as the rhodopsin system associated with vision [40,41]. While describing the interaction with light is not purely a question of electronic structure, as it involves the quantum dynamical evolution of quantum states, the quantum chemical questions revolve around the potential energy surfaces of the ground-and excited-states, and the influence of the environment on the spectrum [42]. These questions are currently challenging due to the size of the systems involved as well as the varying degree of single-and multi-reference character in many of the conjugated excited states [43][44][45].
The basic quantum simulation problem is the ground-state (or low-energy eigenstate) problem for the electronic Hamiltonians, and the basic metric is whether ground-state or lowenergy eigenstate quantum algorithms yield more accurate energies (for comparable computational resources) than the best classical algorithms, for the problem sizes of interest. Proofof-principle demonstrations could be carried out in simplified FIG. 4. Left: Iron-sulfur clusters associated with different oxidation states (P N , Syn (a synthetic model of P N ), P 1+ and P ox ) of the P cluster of the nitrogenase enzyme, from [23]. Right: Mn4Ca core of the oxygen evolving complex of photosystem II, from Ref. [24]. models of the above problems (e.g. in small active spaces of orbitals). However, to make real progress, one should also consider more quantitative models, which requires treating a large number of electrons in a large number of orbitals; at minimum, tens of electrons in hundreds of orbitals. This poses new challenges for ground-state algorithms, and raises questions of how best to represent and encode the resulting Hamiltonians and states. In addition, there are many aspects of the chemical problems beyond the modeling of the electronic wavefunctions, for example, to treat environmental, solvent, and dynamical effects [46]. These will require interfacing the quantum simulation with other classical simulation methods. Finally, although the above examples have focused on multireference and strongly correlated quantum chemistry, weakly correlated chemistry itself becomes hard to model classically when the number of degrees of freedom is very large [47]. These may also be interesting to target with quantum algorithms when sufficiently large quantum machines are available.

B. Quantum molecular spectroscopy
High-resolution gas-phase rovibrational spectroscopy provides an extremely precise experimental probe of molecular structure [48]. Such spectroscopy is important not only for the fundamental understanding of small molecules and the quantum control of atomic and molecular states, but also to provide insight into the basic chemical processes and species involved, for example, in atmospheric chemistry [49] and in astrochemistry [50]. In larger molecules, with more than a few atoms, even the low-energy rovibrational spectrum contains many peaks which cannot be interpreted without theoretical simulation [51]. The theoretical goal is to compute the eigenstates of the nuclear Schrödinger equation [52]. However, unlike the electronic structure problem, there are several challenges even in setting up the best form of the nuclear Schrödinger equation to solve. The first is that the nuclear Hamiltonian (in particular, the nuclear-nuclear interactions) are not known a priori because the interactions are mediated by the electrons. This nuclear potential energy term must instead be determined from quantum chemical calculations at fixed nuclear geometries and then fitted to an appropriate functional form; this requires a large number of high accuracy quantum chemistry calculations. The second is that nuclear ro-vibrational motion is often far from harmonic and not well approximated by simple mean-field theories, unlike many electronic structure problems. Thus there is a need for a proper choice of curvilinear nuclear coordinates that decreases coupling in the nuclear potential energy (e.g. in a harmonic system, normal modes are such a choice of coordinates) while retaining a simple form for the kinetic energy operator, and which also exposes the symmetry of the molecular system.
Once the nuclear Schrödinger equation has been properly formulated, one then faces the problem of representing the eigenstates. While methods such as diffusion Monte Carlo have made progress on vibrational ground-states [53], spectroscopy involves transitions to excited states. In this setting, tensor factorization [54,55] and other approaches [56,57] have been explored to approximate the rovibrational wavefunctions [58]. However, the high dimensionality and spectral congestion, requiring resolution between peaks on the scale of 1 wavenumber, proves extremely challenging [55,59]. Some famous examples include: • Spectra of floppy molecules. Floppy molecules are by their nature very anharmonic and thus far from a simple vibrational description. CH + 5 is a prototypical floppy molecule (see 5): the five hydrogen atoms move around the central carbon and the molecule has almost no structure in the traditional sense [55,[61][62][63][64].
• Hydrogen bonded clusters. Another vibrational problem with large anharmonicity is found in hydrogen bonded clusters, such as in the spectroscopy of water molecules and protonated water clusters [65,66]. The hydrogen bond network is fluid and even small clusters can transition between many different minima on the potential energy surface [67]. Resolving the peaks and tunnelling splittings is important for interpreting water spectra in the atmosphere, as well as in understanding reaction mechanisms in water. Further, the spectroscopy of molecules with intermolecular hydrogen bonds, such as the malonaldehyde molecule, has also posed long-standing challenges for the field [68][69][70][71].
From a quantum algorithms perspective, although there are similarities with the quantum chemistry problem (in particular one is interested in low-energy eigenstates), there are significant differences. One important difference is that the Hamiltonian is no longer of simple two-particle form due to the effective nuclear-nuclear interaction and typically includes important three-and four-mode terms. Also, one is often interested in an order of magnitude more states (e.g. hundreds of excited states) than in the electronic structure problem. All these features are sufficiently distinct from the usual quantum chemical scenarios that quantum algorithms are likely to require additional innovation to be useful in the nuclear problem. Some steps in this direction have recently appeared [48,72,73]. One simplification is that many nuclei are distinguishable avoiding the need to consider indistinguishable particles. The lack of a good mean-field starting point together with the various technical complications means that one can find relatively small systems (in terms of the Hilbert space size) where classical methods already have trouble; for example, a minimal quantum model of the CH + 5 molecule can be formulated as a 12 dimensional problem with 10 basis functions per mode [64].

C. Chemical quantum dynamics
Chemical quantum dynamics is another important target for quantum simulation [74]. This field is concerned with modeling time-dependent electronic and nuclear quantum effects in molecules (as distinct from computing the time-independent electronic or nuclear eigenstates in quantum chemistry and quantum molecular spectroscopy). Quantum molecular dynamics is primarily concerned with the nuclear motion and describes the rates of chemical processes as well as the dynamical interaction of molecules with light, as involved in spectroscopy and quantum control. However, with the devel-opment of short X-ray pulses direct experimental access to electron dynamics in molecules is now also available.
Currently these dynamical simulations are challenging. For example, nuclear motion is poorly described by mean-field theory and the classical limit is often a better starting point, but offers no zeroth order description of quantum effects. Thus classical simulations of quantum dynamics either invoke methods based on the classical limit that scale to large systems but which are difficult to systematically improve (such as approximate path integral methods [75], [76]), or methods which model the wave-function dynamics or the path integral accurately for a small number of degrees of freedom, but which are not scalable due to dimensionality or the dynamical sign problem [77,78].
A subfield of quantum molecular dynamics, but one of important chemical interest, is the description of non-adiabatic quantum effects [74,79]. At nuclear configurations where different electronic surfaces approach each other, the Born-Oppenheimer approximation can break down and the quantum behavior of the nuclei, coupled indirectly via the electrons, is enhanced. The associated quantum non-adiabatic effects govern non-radiative energy relaxation via the crossing of electronic surfaces (so-called conical intersections) and are thus central in describing energy transfer. The faithful description of non-adiabatic quantum effects requires the simultaneous treatment of quantum electrons and quantum nuclei. The complexity of this problem together with the often large system sizes where non-adiabatic effects are of interest means that current classical methods rely on simple heuristic approximations, such as surface hopping [80], for which a rigorous quantum formulation is lacking. Examples of relevant chemical problems in the area of chemical quantum dynamics include: • Proton coupled electron transfer (PCET) [81]. PCET is known to be an important mechanism in catalysis and energy storage: electrons are transferred at lower overpotentials when thermodynamically coupled to proton transfer. Examples range from homogeneous catalysts [82][83][84][85] and heterogeneous electrocatalysts [86,87] to enzymes that perform PCET, such as soybean lipoxygenase (that catalyzes the oxidation of unsaturated fatty acids) [88], photosystem II (in the tyrosine oxidation step) [89,90] and the redox-leveling mechanism of catalytic water oxidation [87]. While semiclassical predictions are often good enough for describing electron transfer, the quantum nature of molecular dynamics is paramount with PCET because of the quantum nature of the proton [88,91] as evidenced in kinetic isotope effects (referring to the ratio of the proton to deuteron reaction rates) which can be very large [88]. Thus, classical mechanics is not sufficient, and quantum simulations of PCET would be extremely helpful for making predictions in large (especially biological) systems.
• Vibrational dynamics in complex environments. For many systems of interest, vibrational spectroscopy is the key tool available for characterization. There is overlap with the problems in section II B but the ques-tions here focus on larger scale systems and condensed phase problems, where the line-shapes as well as frequencies are important, and the system size limits the use of fully quantum formalisms. Indeed, from the librations of water to the high frequency motion of C-H bonds, much of physical chemistry uses nuclear vibrational frequencies to characterize complex systems.
Nevertheless, because of computational limitations, the standard approach today for modeling vibrational dynamics in large complex environments is to invoke a very old flavor of theory: Kubo theory [92]. In this case, one focuses on energy gaps and uses a semiclassical expansion of the line-shape. However, even for the energy gaps one relies on diagonalization of the quantum subsystem and this limits the problems that are accessible; systems with many interacting quantum states -for instance, very large H/J molecular aggregates [93] -are difficult or out of reach.
• Plasmonic chemistry. A recent development is the possibility of using metal particles (with large cross sections) as a tool to absorb light and, with the resulting plasmonic excitations, initiate "plasmonic chemistry". Already, there are a few examples in the literature of hot plasmonic chemistry, including the dissociation of hydrogen molecules [94]. To better understand this arena, however, one of the key questions is: how do we characterize plasmonic excitations? While classical descriptions of plasmons are easy to obtain, quantum descriptions are necessary to model quantum processes, e.g. electron transfer. Covering the disparate length scales of plasmonic excitations and the chemical process remains extremely difficult.
While quantum dynamics is in principle an ideal simulation problem for a quantum computer, the quantum simulation of quantum molecular dynamics entails several practical challenges. Some of this may be viewed as an issue of representation. As in the problem of quantum molecular spectroscopy, the nuclear quantum Hamiltonian contains complicated interactions which must be tabulated or calculated on the fly [95]. In addition to this, the dynamical quantum state involves near continuum degrees of freedom, posing a challenge for the standard discretizations of Hilbert space considered in quantum algorithms. Finally, typical spectroscopic observables may be accessible to relatively simple treatments, not requiring the full fidelity of the quantum wave-function evolution. These technical issues mean that implementing quantum molecular dynamics with a quantum advantage is likely to remain challenging in practice, despite the favorable theoretical complexity on a quantum device.

D. Correlated electronic structure in materials
The goal of electronic structure calculations in materials is to determine their low-energy properties. There is some overlap in methods and ideas between the materials elec-FIG. 6. Qualitative phase diagram of cuprate high-temperature superconductors, which has challenged theory and simulation in condensed matter for decades. Figure adapted from Ref. [96]. tronic structure problem and the problem of quantum chemistry. When electron-electron interactions are weak, the low energy material properties can normally be described by computing the band structure using density-functional theory and one of the many popular density functionals, such as the local density approximation (LDA) or generalized gradient approximations (GGA). However, in some materials, commonly referred to as strongly correlated, the electron-electron interactions fundamentally alter the behavior and such an effective non-interacting description is no longer appropriate. A paradigmatic example are Mott insulators, which appear as conductors in band structure theory but in fact become insulating due to electron-electron interactions. While the mechanism behind the insulating behavior of Mott insulators is wellunderstood, for many other phenomena in strongly correlated systems the underlying microscopic mechanism is not fully known, let alone a quantitative and predictive theory of the associated physics with material specificity. Some famous examples of such problems include: • Originally discovered in 1986, high-temperature superconductivity has eluded a complete theoretical explanation to date. The experimental phase diagram, which is sketched in Fig. 6, is accurately characterized experimentally across several materials. While the general properties of the superconducting phase itself are relatively well characterized, the mechanism driving superconductivity is not yet fully elucidated. Also, two nearby regimes, the pseudogap and strange metal phase, continue to puzzle theorists [96]. In both cases, their nature as well as their precise relation to the superconducting phase are not understood. The strange metal phase (also known as non-Fermi liquid) exhibits behavior inconsistent with a simple weakly interacting metal even at high energies, and has motivated a whole area of research on exotic metallic systems [97]. Meanwhile, the pseudogap phase exhibits several competing ordering tendencies [98], which are extremely challenging to resolve in numerical methods because most methods naturally favor a particular ordering pattern, thus making it challenging to disentangle physical effects from method biases.
• The non-Fermi liquid behavior exemplified in the strange-metal phase of cuprates appears also in other classes of materials, such as heavy fermion compounds and fermionic systems near criticality [99,100]. Often, these systems are amenable to classical simulation only at special points where quantum Monte Carlo methods do not suffer from the infamous sign problem [101].
• Since quantum fluctuations are enhanced in lower dimensions, two dimensional systems have long been of central interest in strongly correlated physics. Many material systems realize effectively two-dimensional physics, including two-dimensional electron gases in semiconductor heterostructures, where the integer and fractional quantum Hall effect were first discovered [102], layered materials (including cuprate high temperature superconductors), graphene [103] and transition-metal dichalcogenides (TMDs). More recently, it has been found that so-called Moiré materials exhibit rich phase diagrams due to strong interactions, including exotic superconductivity [104] and exotic topologically non-trivial phases. A paradigmatic example is twisted bilayer graphene [105], which consists of two graphene layers that are slightly twisted with respect to each other. This leads to a Moiré pattern with a very large unit cell, which effectively quenches the kinetic energy (i.e., leads to almost flat bands) and drastically enhances the effect of Coulomb interaction.
• Frustrated spin systems have long been an important topic especially for numerical simulations in condensed-matter physics. These systems potentially realize a host of high non-trivial phases, in particular topological and gapless spin liquids [106]. They have historically been the testbed for computational methods such as tensor networks and variational methods. As such, they appear as good test cases also for quantum simulations. Furthermore, recent developments in particular in materials with strong spin-orbit coupling have opened the door on a variety of new materials that may exhibit exotic topological phases, and in particular realize a non-Abelian spin liquid [107,108].
While many methods have been developed to accurately include electron-electron interactions, their scope generally remains limited. For example, tensor network methods have revolutionized the study of one-and to a limited extent two-dimensional effective models for magnets and itinerant fermions. However, these methods so far have not been successfully applied to more realistic models, and in particular in three dimensions. On the other hand, quantum embedding methods such as the dynamical mean-field theory (DMFT) and its many cousins can capture interaction effects in threedimensional systems, including multi-band systems. However, they require some approximations to the correlations of the state; for example, in its simplest form, DMFT disregards momentum dependence of the electron self-energy. While many improved variants of these embedding methods exist, their accuracy is often difficult to control, and so far they have not been applied to realistic models without further approximations. Finally, quantum Monte Carlo methods have been extremely successful for bosonic systems and unfrustrated spin systems, but the sign problem hinders their application to frustrated or fermionic systems (away from special points) without other uncontrolled approximations.
From the perspective of quantum algorithms, the materials electronic structure problem is both simpler and more difficult than the quantum chemistry problem. Some ways in which it is simpler include the fact that often very simple Hamiltonians describe the main physics, as well as the potential presence of translational invariance. A major way in which it is more complicated is the fact that one needs to treat systems approaching the thermodynamic limit, which involves a very large number of degrees of freedom. This not only increases the number of qubits required but also heavily impacts the circuit depth of algorithms, such as state preparation. The thermodynamic limit can also lead to small energy scales for excitations and energy differences between competing phases. For these reasons, it remains to be understood whether the quantum algorithms of relevance to quantum chemistry are those of relevance to materials electronic structure.

E. Dynamical quantum effects in materials
Many experiments on condensed-matter systems do not probe the equilibrium properties of the system, but rather dynamical properties. For example, the main workhorse of mesoscopic quantum physics is electron transport, i.e. the response of the system when it is coupled to electron reservoirs and a voltage is applied [110]. Likewise, material properties are often probed through scattering experiments, such as neutron scattering or angle-resolved photoemission spectroscopy (ARPES) [111], which probe dynamical properties such as the structure factor or spectral function. Going beyond spectral properties, the non-equilibrium real-time dynamics of quantum systems have increasingly come into focus, both because of experiments that can probe quantum dynamics at atomic scales and because of fundamental interest in equilibration of quantum systems. Experimental setups that can probe ultra-fast dynamics in materials include, for example, FIG. 8. Cold atomic gases are controllable quantum systems that can be arranged in a lattice (a). They can be tuned to provide analog simulators of real crystal potentials (b). From Ref. [117].
free-electron lasers [112,113] as well as other pulsed laser systems. These allow the application of experimental techniques, such as pump-probe spectroscopy [109,114], to provide novel insights into the behavior of correlated quantum systems.
On the other hand, cold atomic gases [115] provide a highly controllable environment that allows systematic exploration of quantum dynamics even in the strongly interacting regime [116], see 8. A key advantage is that one can engineer the evolution of the system to closely follow a target model; this approach is also referred to as analog quantum simulation (see also Sec. I B). However, classical simulation still plays a crucial role in establishing the accuracy of cold atom setups.
From a conceptual point of view, a central question has become the connection between statistical mechanics and the dynamics of closed quantum systems. The general goal is to put quantum statistical mechanics on a solid conceptual foundation as well as understanding the cases where it does not apply, such as many-body localized systems [118][119][120].
Numerically simulating all these systems has been a severe challenge. While many approaches to quantum dynamics exist, none are universally applicable. At short times, tensor-network methods can accurately describe the dynamics. However, in general the computational cost grows exponentially with the desired simulation time, thus severely limiting the timescales that can be resolved [121]. Conversely, at long times, the system is often effectively described by classical dynamics controlled by conservation laws. However, the interesting strongly correlated behavior is generally exhibited at intermediate times, inaccessible to the established classical methods. For example, quantum Monte Carlo methods can scale exponentially in time even for unfrustrated systems. Non-equilibrium dynamical mean-field theory [122] has emerged as a powerful method especially for systems in high dimensions, but requires uncontrolled approximations (both in the setup of the method and its numerical solution). Finally, in the regime of weak interactions, time-dependent density functional theory can be used, but likewise implies uncontrolled approximations [123].

A. Overview of algorithms
In Section I A we described three quantum problems that lie at the heart of chemistry and materials physics: the problem of quantum dynamics, representative of computational tasks that can be efficiently tackled by a quantum computer, aside from the initial state preparation; quantum ground-and low-energy state determination, central to quantum chemistry, condensed phase electronic structure, and quantum molecular spectroscopy; and thermal averages or quantum statistical mechanics, to describe finite-temperature chemistry and physics. In next sections, we survey the current status and theoretical and practical challenges to implement quantum algorithms for these problems, and some future directions.
Regardless of the problem we are studying, the first step in a quantum simulation is to choose a representation for the Hamiltonian and the states. Thus, in III B we first examine the possibilities for different qubit representations, and open questions in that area.
Quantum ground-state algorithms fall into several classes. Quantum phase estimation is a direct route to (near exact) eigenstate determination, but has been challenging to implement in the near-term era. A complementary technique is to prepare the exact ground-state via a prescribed "exact" evolution path, either in real time (adiabatic state preparation), or in imaginary time (quantum imaginary time evolution). On the other hand, variational quantum algorithms provide the possibility to introduce approximations with adjustable circuit depth, which are determined via a non-linear optimization, usually implemented in a hybrid-quantum classical loop. The variety of algorithms currently explored for quantum ground-(and low-lying excited state) determination are discussed in III C.
While polynomial time algorithms for quantum time evolution have been known for some time, algorithms with optimal asymptotic complexity (linear in time) as well as favorable scaling with error and low prefactors, remain an area of active research. Also, much work remains to be done to optimize Hamiltonian simulation algorithms for the Hamiltonians of interest in chemistry and materials science. Quantum time evolution is a fundamental building block in many quantum algorithms, such as phase estimation. The current status of quantum time evolution algorithms is summarized in III H.
How best to simulate thermal states in chemical and materials science applications remains an open question. A wide variety of techniques have been discussed, ranging from eigenstate thermalization, to state preparation methods, to hybrid quantum-classical algorithms, though few have been implemented. The current status of thermal state methods, the prospects for implementing them, and other open questions are discussed in III K.
Many quantum algorithms involve interfacing with classical data and classical algorithms. This can be to take advantage of classical optimization strategies, as in variational quantum algorithms. Another reason is to enable a multi-level/multi-scale representation of the problem. Quantum embedding provides a framework for such multi-scale quantum/classical hybrids, with the quantum representation of a subsystem coupled either to a classical environment, or another quantum representation via the exchange of classical data. We discuss the current status of hybrid quantumclassical algorithms and quantum embedding in particular in III L.
Finally, an important consideration when developing improved quantum algorithms for real chemical and materials science problems is to establish benchmark systems and results, from the best available classical simulation data. The possibilities and prospects for such benchmarks are discussed in IV.

B. Qubit representation of many-body systems
Many-body systems in chemistry and materials physics are composed of interacting electrons and atomic nuclei. An exact quantum mechanical treatment involves continuous variables such as the particles' positions and momenta. To simulate such systems on a digital computer (either quantum or classical), the infinite-dimensional Hilbert space of a many-body system has to be truncated.
The most direct route is to define a finite set of basis functions and then to project the exact many-body Hamiltonian onto the chosen basis. The resulting discretized system is then expressed in terms of qubits. Depending on the problem, the Hamiltonian of interest may be different, e.g. in electronic structure it is the electronic Hamiltonian, while in molecular vibrational problems, it is the nuclear Hamiltonian. Alternatively, one can write down a simple form of the Hamiltonian a priori that contains the main interactions (a model Hamiltonian) with adjustable parameters. This latter approach is particularly popular in condensed matter applications. Finally, depending on the particles involved it may also be necessary to account for their fermionic or bosonic nature, in which case a suitable encoding of the statistics is required.
A choice of a good representation is important as it may affect the simulation cost dramatically. In this section we briefly summarize known methods for the qubit representation of many-body systems, discuss their relative merits, and outline important directions for future research.

Ab initio electronic structure qubit representations
The main objective of electronic structure in chemistry and physics to understand the low-energy properties of the electronic structure Hamiltonian that describes a system of interacting electrons moving in the potential created by atomic nu-clei [124,125], Here K is the number of electrons, r i is the position operator of the i-th electron, ∆ i is the corresponding Laplacian, and V (r) is the electric potential created by atomic nuclei at a point r. The termĤ 1 includes the kinetic and the potential energy of non-interacting electrons whileĤ 2 represents the Coulomb repulsion. Here we ignore relativistic effects and employ the standard Born-Oppenheimer approximation to solve the electronic Hamiltonian for fixed nuclei positions. Each electron is described quantum mechanically by its position r i ∈ R 3 and spin ω i ∈ {↑, ↓}. Accordingly, a quantum state of K electrons can be specified by a wave function The wave function must obey Fermi statistics, that is, Ψ(x 1 , . . . , x K ) must be antisymmetric under exchanging any pair of coordinates x i and x j .
The first step of any quantum electronic structure simulation algorithm is to approximate the electronic HamiltonianĤ with a simpler simulator Hamiltonian that describes a system of interacting qubits. This is usually achieved by truncating the Hilbert space of a single electron to a finite set of basis functions ψ 1 , . . . , ψ N known as (spin) orbitals. For example, each orbital could be a linear combination of atom-centered Gaussian functions with a fixed spin orientation.
Electronic structure simulation algorithms based on the first quantization method [126-128] describe a system of K electrons using the Configuration Interaction (CI) space (in classical simulations, this would be called the Full Configuration Interaction space). This is a linear space spanned by all possible Slater determinant states that can be formed by distributing K electrons over N orbitals. The CI space has dimension N K and can be identified with the anti-symmetric subspace of (C N ) ⊗K .
The projection of the full electronic HamiltonianĤ onto the CI space has the form Here |p ≡ |ψ p are the chosen spin-orbitals. The coefficients t p,q and u pqrs are known as one-and two-electron integrals. For example, Likewise, u pqrs is the matrix element of the Coulomb interaction operator 1/|r 1 − r 2 | between anti-symmetrized versions of the states ψ p ⊗ ψ q and ψ r ⊗ ψ s . Each copy of the singleelectron Hilbert space C N is then encoded by a register of log 2 N qubits. This requires n = K log 2 N qubits in total. The CI Hamiltonian H includes multi-qubit interactions among subsets of 2 log 2 N qubits. The full Hilbert space of n qubits contains many unphysical states that do not originate from the CI space. Such states have to be removed from simulation by enforcing the anti-symmetry condition. This can be achieved by adding suitable energy penalty terms to the CI Hamiltonian [129].
An important parameter that affects the runtime of quantum simulation algorithms is the sparsity of the simulator Hamiltonian. A Hamiltonian H is said to be d-sparse if the matrix of H in the standard n-qubit basis has at most d non-zero elements in each row (each column). For example, the runtime of simulation algorithms based on quantum signal processing [130] scales linearly with the sparsity d. The CI Hamiltonian H has sparsity d ∼ (KN ) 2 . Thus the first-quantization method is well-suited for high-precision simulation of small molecules when the number of electrons K = O(1) is fixed and the number of orbitals N is a large parameter. As one approaches the continuum limit N → ∞, the number of qubits grows only logarithmically with N while the sparsity of H scales as d ∼ N 2 .
The second quantization approach often results in a simpler simulator Hamiltonian and requires fewer qubits, especially in the case when the filling fraction K/N is not small. This method is particularly well suited for quantum simulation algorithms [131] and has been experimentally demonstrated for small molecules [132]. Given a set of N orbitals ψ 1 , . . . , ψ N , the second-quantized simulator Hamiltonian is whereĉ † p andĉ p are the creation and annihilation operators for the orbital ψ p . The Hamiltonian H acts on the Fock space spanned by 2 N basis vectors |n 1 , n 2 , . . . , n N , where n p ∈ {0, 1} is the occupation number of the orbital ψ p . The advantage of the second quantization method is that the Fermi statistics is automatically enforced at the operator level. However, the number of electrons can now take arbitrary values between 0 and N . The simulation has to be restricted to the subspace with exactly K occupied orbitals. The secondquantized Hamiltonian H can be written in terms of qubits using one of the fermion-to-qubit mappings discussed in Section III B 3.
Within the above outline, there are several active areas of research. For example, one may wonder if the redundancy of the space in the first quantized representation can be reduced or completely avoided. Other questions include the choice of basis functions and fermion-to-qubit mappings. These are discussed in the next sections.

Electronic basis functions
The discretization of the electronic Hilbert space for a quantum simulation requires balancing two concerns. We need to represent the state with as few qubits, but also, retain maximal Hamiltonian sparsity. These requirements only partially align with those of classical many-particle quantum simulations. In the classical setting, a compact state is crucial due to the exponential Hilbert space, while Hamiltonian sparsity is less so; the choice of basis functions has historically been made so that matrix elements of the Hamiltonian (the one-and two-electron "integrals") can be analytically evaluated [124,133].
There are two families of basis functions in wide use in quantum chemistry and quantum materials science: atomic orbital Gaussian bases and plane waves. Gaussian bases are most commonly employed in molecular simulations due to their compactness, while plane waves are most often used in crystalline materials simulation, due to their intrinsic periodicity and ease of regularizing the long-range contributions of the Coulomb operator (which are conditionally convergent in an infinite system).
In Gaussian bases, linear combinations of Gaussian functions (referred to as simply Gaussian basis functions) are placed at the nuclear positions. As they are placed where the ground-state electron density is highest, they give a compact representation of the wavefunction for bound states, but the Hamiltonian is not sparse, with O(N 4 ) second-quantized matrix elements. In a quantum algorithm, this leads to high gate counts even for simple quantum primitives such as a Trotter step.
Plane-waves offer greater simplicity as the accuracy of the basis is controlled by a single parameter, the kinetic energy cutoff. While the number of plane waves needed to reach a desired accuracy is larger than the number of Gaussian states, the Hamiltonian contains fewer terms (O(N 3 )) due to momentum conservation. To reduce the number of required plane waves, it is essential to employ pseudopotentials to remove the sharp nuclear cusp [135,136]. Furthermore, the asymptotic basis convergence of Gaussian and pseudopotential plane wave calculations is the same: the feature governing the rate of convergence is the wavefunction discontinuity or electron-electron cusp due to the singularity of the Coulomb interaction (see 9). In classical simulations, so-called explicit correlation methods can be used to remove the slow convergence due to the singularity [137,138]. How to use such techniques with quantum computers has yet to be explored.
The need to expose more sparsity in the Hamiltonian while retaining a reasonably compact wavefunctions is an active area of research in both classical and quantum algorithms. Recent ideas have included new types of basis function that return to a more grid-like real-space basis [139][140][141][142][143] where the Coulomb operator and thus Hamiltonian has only a quadratic number of terms, as well as factorizations of the Coulomb operator itself [144,145]. The best choice of basis for a quantum simulation remains very much an open question. FIG. 9. Extrapolation of the Hartree-Fock (left) and correlation energy (right) for chains of 10 hydrogen atoms, using cc-pVxZ Gaussian bases, and various methods from Ref. [134]. While the mean-field energy converges exponentially fast in the basis size (roughly given by x 3 ), the correlation energy converges as x −3 (inversely proportional to the number of basis functions) due to the electron-electron cusp.

Fermion-to-qubit mappings
Since the basic units of a quantum computer are qubits rather than fermions, any quantum simulation algorithm of fermions (e.g. for electronic structure) employs a suitable encoding of fermionic degrees of freedom into qubits. For example, the standard Jordan-Wigner mapping (sketched in 10) identifies each Fermi mode (orbital) with a qubit such that the empty and the occupied states are mapped to the qubit basis states |0 and |1 respectively. More generally, the Fock basis vector |n 1 , . . . , n N is mapped to a qubit basis vector |x 1 ⊗ · · · ⊗ |x N , where each bit x j stores a suitable partial sum (modulo two) of the occupation numbers n 1 , . . . , n N . The Jordan-Wigner mapping corresponds to x j = n j for all j. This is not quite satisfactory since single-mode creation/annihilation operators contain Jordan-Wigner strings (−1) n1+...+nj which may be very non-local in terms of qubits. On the other hand, updating the qubit state x upon application of a single creation/annihilation operator requires a single bit flip (see 11). More efficient fermion-toqubit mappings balance the cost of computing Jordan-Wigner strings and the bit-flip cost of updating the qubit state. For example, the binary tree encoding proposed by Kitaev and one of the authors [146] maps any fermionic single-mode operator (e.g.ĉ p orĉ † p ) to a qubit operator acting non-trivially on roughly log N or less qubits. Generalizations of this encoding were studied in [147][148][149]. As a consequence, the secondquantized Hamiltonian Eq. (7) expressed in terms of qubits becomes a linear combination of Pauli terms with weight at most O(log N ). This is important in the context of VQE-type quantum simulations (see III E) since Pauli operators with an extensive weight (of order N ) cannot be measured reliably in the absence of error correction.
A natural question is whether the number of qubits required to express a Fermi system can be reduced by exploiting symmetries such as the particle number conservation or the point group symmetries of molecules. For example, zero temperature simulations often target only one symmetry sector containing the ground state. This motivates the study of symmetry-adapted fermion-to-qubit mappings. The goal here is to reduce the number of qubits required for the simulation without compromising the simple structure of the resulting qubit Hamiltonian (such as sparsity). The simplest case of Z 2 symmetries is now well understood and the corresponding symmetry adapted mappings are routinely used in experiments [129,132,150,151]. The U (1) symmetry underlying particle number conservation was considered in [129,152,153].
Dimension counting shows that a system of N Fermi modes with exactly K particles can be mapped to roughly log 2 N K FIG. 11. Difference between Jordan-Wigner and parity encoding. The former encodes occupation numbers xi on qubit states, the latter parities pi = j<i xi mod 2. Other strategies, such as the binarytree encoding, balance non-locality of occupation numbers and parities to achieve more efficient encodings.
qubits. However, it remains an open question whether this mapping can be chosen such that the resulting qubit Hamiltonian admits a sparse representation in some of the commonly used operator bases (such as the basis of Pauli operators) to enable applications in VQE. Mappings adapted to point group symmetries have been recently considered in [154]. It is also of great interest to explore fermion-to-qubit mappings adapted to approximate and/or emergent symmetries.
Alternatively, one may artificially introduce symmetries either to the original Fermi system or its encoded qubit version with the goal of simplifying the resulting qubit Hamiltonian. This usually requires redundant degrees of freedom such as auxiliary qubits or Fermi modes [146,[155][156][157]. In the case of lattice fermionic Hamiltonians such as the 2D Fermi Hubbard model or more general models defined on bounded degree graphs, such symmetry-adapted mappings produce a local qubit Hamiltonian composed of Pauli operators of constant weight independent of N .

Model and non-electronic problems
While much of the above has focused on ab initio quantum chemistry and electronic structure in quantum simulations, the diverse questions of chemistry and materials physics discussed in section II raise additional issues of representation. For example, model Hamiltonians avoid the problem associated with a choice of basis by restricting the Hamiltonian to a predetermined simple form on a lattice. The lattice structure permits specialized techniques, such as specific fermionic encodings. Developing specialized representations for model problems is of particular importance in simulating condensed matter systems.
Other kinds of non-electronic simulations may involve different requirements on the basis than electronic problems. For example, in quantum reactive scattering processes, there is little a priori information on the positions of the particles; instead various grid representations, often in non-Cartesian coordinate systems, are used [158,159]. Alternatively, the particles of interest may be bosons which engender new encoding considerations. Relatively little attention has been paid to these questions so far in quantum algorithms.

C. Quantum algorithms for ground and excited states
There are many approaches to obtaining ground states or excited states on a quantum computer. State preparation procedures attempt to construct a circuit to prepare a state with as large as possible overlap with the desired eigenstate. One set of such procedures, which we review in Sec. III D and which includes adiabatic state preparation and quantum imaginary time evolution, uses a prescribed evolution path. An alternative strategy is based on variational methods (often called variational quantum eigensolvers) and is reviewed in Sec. III E. Here, the preparation circuit itself is defined via the optimization of the energy with respect to parameters of the circuit.
Given some state with reasonably large overlap with the desired state, one can perform quantum phase estimation, which simultaneously projects the state onto an eigenstate of the Hamiltonian and obtains an estimate for the energy of this eigenstate. The error depends inversely on the simulation time. The probability of successfully projecting onto the desired state is given by the square overlap of the input state and the desired state, and it is thus necessary to use some other method (such as the state preparation procedures above) to prepare an input state with sufficient overlap with the desired state.
The various approaches come with different strengths and weaknesses. While phase estimation allows the deviation of the final state from an exact eigenstate (although not necessarily the desired eigenstate) to be systematically reduced, it can require deep circuits with many controlled gates that are challenging for devices with limited coherence and without error correction. Variational methods or quantum imaginary time evolution replace such circuits by a large number of potentially shorter simulations, which is expected to be easier to implement on near-term machines. However, if one does not measure the energy by phase estimation, but instead by expressing the Hamiltonian as a sum of multi-qubit Pauli operators and measuring the terms individually, the state preparation and measurements must be repeated many times, with the error converging only as the square root of the number of repetitions. Variational methods are also limited by the variational form and ability to solve the associated optimization problem, which may by itself represent a difficult classical optimization. Finally, adiabatic state preparation and quantum imaginary time evolution become inefficient for certain Hamiltonians.
State preparation methods are first discussed in III D and III E (variational state preparation is discussed separately due to the large number of different types of ansatz). Considerations for excited states are discussed in III F. Phase estimation is reviewed in III G. D. Preparing ground states along a prescribed path

Adiabatic state preparation
One general route to prepare the ground-state of a physical system on a quantum device is through adiabatic state preparation. This relies on the well-known adiabatic theorem [160,161], which states that a system evolved under a time-dependent Hamiltonian will remain in the instantaneous ground state as long as the evolution is sufficiently slow and the spectrum of the Hamiltonian remains gapped. To make use of this, one chooses a Hamiltonian pathĤ(λ), 0 ≤ λ ≤ 1, such that the ground state ofĤ(0) is easily prepared, whilê H(1) is the Hamiltonian whose ground state one wants to obtain. The system is then evolved under the time-dependent Schrödinger equation, In the limit T → ∞ and if the spectrum ofĤ t/T is gapped for all t, the final state |Φ T is the exact ground state ofĤ. Away from the adiabatic limit T → ∞, corrections arise that depend on the instantaneous gap ofĤ(λ) and the total time T . Furthermore, if degeneracies occur along the path, a different time-dependent Hamiltonian path (although with the same endpoints) must be chosen. Thus while this approach is very general, its practical applicability is limited by the requirement of having to take the limit of large T , and to choose a path without degeneracies (see 12). The limit of large T may require deep circuits that may not be practical in near-term quantum machines. Analyzing the errors (and optimizing the path, for example by choosing an improved f (s)) is also challenging due to the dependence on the unknown spectrum ofĤ t/T . Some of these questions have been studied in the more general context of adiabatic quantum computation [162]. However, while this is one of the first state preparation methods discussed for chemical systems [131], more heuristic work in this area for problems of interest to chemistry and physics is needed [163].

Quantum imaginary-time evolution
In classical simulations, one popular approach to prepare (nearly exact) ground-states is imaginary-time evolution, which expresses the ground-state as the long-time limit of the imaginary-time Schrödinger equation, Imaginary time-evolution underlies the family of projector quantum Monte Carlo methods in classical algorithms [125]. To perform imaginary time evolution on a quantum computer, it is necessary to implement the (repeated) action of the short-imaginary-time propagator e −∆τĤ on a state. Given a Hamiltonian that can be decomposed into geometrically local terms,Ĥ = mĥ m , and a state |Ψ with finite correlation length C, the action of e −∆τĥi can be generated by a unitarŷ U = e iÂ acting on O(C) qubits surrounding those acted on byĥ i , i.e.
where the coefficients of the Pauli strings inÂ can be determined from local measurements of the qubits aroundĥ i . This is the idea behind the quantum imaginary time evolution (QITE) algorithm [164]. Like adiabatic state preparation, quantum imaginary time evolution can in principle prepare exact states without the need for variational optimization. Also, the total length of imaginary time propagation to achieve a given error is determined by the spectrum ofĤ and the initial overlap, rather than by the spectrum ofĤ t/T along the adiabatic path. However, the method becomes inefficient in terms of the number of measurements and complexity of the operatorÂ if the domain C grows to be large along the imaginary time evolution path. In these cases, QITE can be used as a heuristic for approximate ground-state preparation, analogous to using adiabatic state preparation for fixed evolution time. While initial estimates in a limited set of problems show QITE to be resource efficient compared to variational methods due to the lack of an optimization loop [164], a better numerical understanding of its performance and cost across different problems, as well as the accuracy of inexact QITE in different settings, remains to be developed.

E. Variational state preparation and variational quantum eigensolver
A class of state preparation methods that have been argued to be particularly amenable to near-term machines is variational state preparation [165][166][167]. Here, similar to classical variational approaches, one chooses a class of ansatz states for the ground state of the Hamiltonian of interest. Generally speaking, such an ansatz consists of some initial state and a unitary circuit parametrized by some set of classical variational parameters. Applying this circuit to the initial state yields a guess for the ground state, whose energy is then evaluated. This yields an upper bound to the true ground state energy. One then varies the variational parameters to lower the energy of the ansatz state.
In choosing the class of ansatz states, one pursues several goals: on the one hand, it is crucial that the class contains an accurate approximation to the true ground state of the system. On the other hand, one desires a class of circuits that are easily executed on the available quantum computer, i.e. for a given set of available gates, connectivity of the qubits, etc. Finally, it is important for the classical optimization over the variational parameters to be well-behaved, so as to be able to find lowenergy minima. While we cannot list all possible ansatz states below, we provide a representative sample.

Unitary coupled cluster
An early example of a particular class of ansatz states that has been suggested for applications in quantum chemistry is the unitary coupled-cluster (uCC) ansatz [166,168,169] Here, d denotes the maximum order of excitations in the uCC wavefunction (for example d = 1, 2, 3 for singles, doubles and triples respectively),ĉ † a k . . .ĉ † a1 (ĉ i k . . .ĉ i1 ) are creation (destruction) operators relative to orbitals unoccupied (occupied) in the Hartree-Fock state, and t is a rank-2k tensor, antisymmetric in the a k . . . a 1 and i k . . . i 1 indices. This choice of ansatz is motivated by the success of mean-field theory, which suggests that the density of excitations in the true wavefunction should be small relative to the mean-field state. Standard coupled cluster theory -written as eT |Ψ HF -is widely used in classical quantum chemistry but is challenging to implement on a quantum device, whereas the reverse is true for the unitary variant. Understanding the theoretical and numerical differences between standard and unitary coupled cluster is an active area of research [170,171]. The variational quantum eigensolver algorithm applied to the unitary coupled-cluster Ansatz is depicted in 13.

Hardware-efficient ansatz
The unitary coupled cluster ansatz involves non-local gate operations and is expensive to implement on near-term devices with limited qubit connectivity. An alternative variational approach, pursued e.g. in Ref. [132] and termed "hardware-efficient" there, is to tailor the ansatz specifically to the underlying hardware characteristics. The circuits considered in Ref. [132], sketched in 14, consist of alternating layers of arbitrary single-qubit gates and an entangling gate that relies on the intrinsic drift Hamiltonian of the system. While this drift Hamiltonian and thus the entangling gate is not known precisely, for the variational approach it is sufficient to know that the gate is reproducible. The variational parameters are only the rotations in the layer of single-qubit gates. While it is not guaranteed that such an ansatz contains a good approximation to the state of interest, it is an example of an adaption of a method to NISQ devices [132]. An application to BeH 2 is seen in 14.

Adapt-VQE ansatz
In the adapt-VQE scheme, a collection of operatorsÂ i (operator pool) is chosen in advance, and the ground state is approximated by |Ψ adapt-VQE = e θnÂn . . . e θ1Â1 |Ψ HF , Given a current parameter configuration θ 1 . . . θ n , the commutator of the Hamiltonian with each operator in the pool is measured to obtain the gradient of the energy with respect to the parameters θ. Repeating this multiple times and averaging over the obtained samples gives the gradient of the expectation value of the Hamiltonian with respect to the coefficient of each operator. The ansatz is improved by adding the operatorÂ i with the largest gradient to the left end of the ansatz with a new variational parameter, thereby increasing n. The operation is repeated until convergence of the energy [173]. Numerical simulations, for example for short hydrogen chains, show that adapt-VQE can improve over the unitary coupled cluster approach in terms of the accuracy reached for a given circuit depth.

Tensor networks
Tensor networks are a class of variational states which construct the global wavefunction amplitude from tensors associated with local degrees of freedom. They specify a class of quantum states that can be represented by an amount of classical information proportional to the system size. There are two main families of tensor networks: those based on matrix product states and tree tensor network states (MPS, TTNS) [174,175] (also known in the numerical multi-linear algebra community as the tensor-train decomposition [176,177] and the tree-structured hierarchical Tucker representation [178,179]) and their higher dimensional analogs, projected entangled pair states (PEPS) [180]; and those based on the multi-scale entanglement renormalization ansatz (MERA) [181]. Because of their success in representing low-energy states in classical simulations, they are a natural class of variational wavefunctions to try to prepare in a quantum algorithm for low-energy states. These tensor networks are schematically depicted in 15. There are many analogies between tensor network algorithms and quantum circuits. This analogy can be exploited to develop an efficient preparation mechanism for these states on a quantum computer. By recognizing that the tensors in an MPS or a MERA can be associated with a block of unitaries (with the bonds between tensors playing the role of circuit lines in a quantum circuit) it is possible to prepare an MPS or MERA state on a quantum computer [182][183][184]. Because the dimensions of the associated tensor grow exponentially with the depth of the quantum circuit associated with it, it is possible to prepare certain tensor networks with large bond dimension on a quantum computer that presumably cannot be efficiently simulated classically; an example of this is the socalled deep MERA [185].
There are many open questions in the area of tensor networks and quantum computing. For example, preparing PEPS on a quantum computer appears to be much less straightforward than preparing a matrix product state or a MERA [186]. Similarly, although "deep" tensor network states can only be efficiently simulated on a quantum computer, their additional representational power over classically efficient tensor networks for problems of physical or chemical interest is poorly understood.

Other considerations
Besides the choice of ansatz state, the computational challenges of variational methods and VQE are twofold: • Potentially, a very large number of measurements must be performed to accurately estimate the energy. Indeed, the scaling is quadratically worse than when using quantum phase estimation. We will discuss this point further in V.
• The optimization of the variational parameters may be very difficult, in particular if the energy exhibits a very non-trivial dependence on classical parameters with many local minima, and if gradient information is not easily available. For some discussion of optimization algorithms in this context, see Refs. [188,189].
Some key advantages of the VQE approach are that it can often be carried out with a large number of independent, short quantum simulations. This is more suitable to NISQ machines than the long coherent circuits required for approaches based on quantum phase estimation, which has been demonstrated in several experiments [132,150]. Furthermore, the approach is more resilient against certain types of errors. For example, as mentioned already above, it is generally not necessary to know exactly what circuit is executed for some variational parameters as long as it is reproducible; therefore, systematic coherent tuning errors of the qubits (for example systematic deviations between the desired and the actually applied single-qubit rotations) do not adversely affect the results. In addition to studying the robustness of VQE against errors, it has become a very active field to develop techniques that mitigate such physical errors. Such approaches promise to reduce the impact of errors on near-term machines before error correction becomes available. For work in this direction, see Refs. [188,[190][191][192][193][194][195].  [187] and Ref. [185].

F. Excited states
While much of the above discussion of state preparation and variational algorithms has focused on ground-states, most of the same methods can also be used with minor extensions for excited states. For example, adiabatic state preparation can be used to prepare an excited state, so long as it is connected to the initial state without a vanishing gap.
In the area of variational methods, it is often useful to choose the excited state ansatz to be related to that of the ground-state, since at low-energies much of the physics is the same. This is widely used in classical simulations and essentially the same ideas have been ported to the quantum algorithm setting. For example, in the quantum subspace expansion (QSE) [196], the excited state is made via the ansatz |Ψ = α c αÊα |Ψ , where {E α } is a set of "excitation" operators and |Ψ is the ground-state constructed in VQE. In QSE, one needs to measure all the subspace matrix elements Ψ|Ê † αÊβ |Ψ , Ψ|Ê † αĤÊβ |Ψ , thus the number of measurements grows quadratically with the subspace. In the quantum Lanczos method, the QITE algorithm is used to construct the subspace {e −λĤ |Ψ , e −2λĤ |Ψ , . . .} and the special structure of this space means that all subspace matrix elements can be constructed with a number of measurements that grows only linearly with the size of the subspace [164]. Alternatively, one can fix the coefficients c α and reoptimize the quantum circuit in the variational method; this is the basis of the multi-state VQE method; other similar ideas have also been proposed. Connections between quantum subspaces and error correction have been explored in [197].
The above methods compute total energies of excited states, which have to be subtracted from the ground state energy to give the excitation energies of the system. A method to directly access excitation energies is desirable. One route to achieve this goal is provided by the equation-of-motion (EOM) approach, also widely used in classical simulations [198] and recently extended to quantum computing [199]. In the EOM approach, excitation energies are obtained as where Ψ is an approximation to the ground state (such as the VQE ansatz) andÔ n is an excitation operator expanded on a suitable basis. The variational problem of finding the stationary points of ∆E n leads to a generalized eigenvalue equation, the solutions of which are the excited-state energies.
FIG. 16. Quantum circuit for iterative QPE. The first two singlequbit gates bring the ancilla to the state 1 √ 2 |0 + e iθ |1 . The controlled-U (t) gate, with U (t) = e −itĤ , brings the register into the state (15). The last Hadamard gate is needed to measure in the X basis.

G. Phase estimation
Quantum Phase Estimation (QPE) is a crucial step in many quantum algorithms. In the context of quantum simulation, QPE enables high-precision measurements of the ground and excited energy levels. This is achieved by preparing a trial initial state |ψ(0) that has a non-negligible overlap with the relevant eigenvector of the target HamiltonianĤ and applying a quantum circuit that creates a superposition of time evolved states |ψ(t) = e −iĤt |ψ(0) over a suitable range of the evolution times t. In the simplest case, known as the iterative QPE [4,131] and sketched in 16, the final state is a superposition of the initial state itself and a single time-evolved state, Here one ancillary control qubit has been added that determines whether each gate in the quantum circuit realizing e −iĤt is turned on (control is 1) or off (control is 0). The extra phase shift θ coordinates interference between the two computational branches such that useful information can be read out with high confidence. Finally, the control qubit is measured in the so-called X-basis, |± = (|0 ± |1 )/ √ 2 and the measurement outcome b ∈ {+1, −1} is recorded.
The iterative QPE works by performing many runs of the above subroutine with a suitable choice of parameters t, θ in each run and performing a classical post-processing of the observed measurement outcomes. The ancillary control qubit stays alive only over the duration of a single run since its state is destroyed by the measurement. However, the remaining qubits that comprise the simulated system stay alive over the entire duration of the QPE algorithm. More specifically, let E α and |ψ α be the eigenvalues and eigenvectors ofĤ such thatĤ = α E α |ψ α ψ α |. The trial state can be ex-panded in the eigenbasis ofĤ as |ψ(0) = α c α |ψ α . Then the joint probability distribution describing measurement outcomes b 1 , . . . , b N = ±1 observed in N runs of QPE has the form (16) Here θ i and t i are the phase shift and the evolution time used in the i-th run and b i = ±1 is the observed measurement outcome. This has the same effect as picking an eigenvector ofĤ at random with the probability |c α | 2 and then running QPE on the initial state |ψ(0) = |ψ α . Accordingly, QPE aims at estimating a random eigenvalue E α sampled from the probability distribution |c α | 2 . We note that a suitably formalized version of this problem with the trial state |00 . . . 0 and a local Hamil-tonianĤ composed of few-qubit interactions is known to be BQP-complete [200]. In that sense, QPE captures the full computational power of quantum computers and any quantum algorithm can be expressed as a special case of QPE. A common application of QPE is the task of estimating the smallest eigenvalue E 0 = min α E α . This requires a trial state |ψ(0) that has a non-negligible overlap with the true ground state ofĤ to ensure that the minimum of a few randomly sampled eigenvalues E α coincides with E 0 . For example, in the context of molecular simulations, |ψ(0) is often chosen as the Hartree-Fock approximation to the ground state. Such a state is easy to prepare as it can be chosen to correspond to a standard basis vector, see III B 3.
The problem of obtaining a good estimate of the eigenvalue E α based on the measured outcomes b 1 , . . . , b N is an active research area, see [201][202][203][204]. Assuming that the trial state has a constant overlap with the ground state ofĤ, the smallest eigenvalue E 0 can be estimated using QPE within a given error using N ∼ log (1/ ) runs such that each run evolves the system over time at most O(1/ ). As discussed in III H, the evolution operator e −iĤt (as well as its controlled version) can be approximated by a quantum circuit of size scaling almost linearly in t (neglecting logarithmic corrections). Thus QPE can achieve an approximation error at the computational cost roughly 1/ even if the trial state has only a modest overlap with the ground state. This should be contrasted with VQE algorithms that have cost at least 1/ 2 due to sampling errors and where the trial state must be a very good approximation to the true ground state, see III E. On the other hand, QPE is much more demanding in terms of the required circuit depth and the gate fidelity. It is expected that quantum error correction will be required to implement QPE in a useful way (e.g. to outperform VQE in ground-state determination).
QPE also has a single-run (non-iterative) version where the time evolution of the simulated system is controlled by a multi-qubit register and the X-basis measurement of the control qubit is replaced by the Fourier basis measurement [205]. The iterative version of QPE has the clear advantage of requiring fewer qubits. It also trades quantum operations required to realize the Fourier basis measurement for classical postprocessing, thereby reducing the overall quantum resource cost of the simulation.
Since QPE is used ubiquitously in a variety of quantum applications, it is crucial to optimize its performance. Below we list some open problems that are being actively investigated; see Ref. [204] for a recent review.
• Given limitations of near-term quantum devices, of particular interest are tradeoffs between the depth of the QPE circuit and its spectral-resolution power as well as its sensitivity to noise. It was shown [204] that that the computational cost of QPE interpolates between 1/ and 1/ 2 as the depth (measured by the number of runs N per iteration of the algorithm) is reduced from O(1/ ) to O(1). A particular version of QPE with a tunable depth that interpolates between the standard iterative QPE and VQE was proposed in [206].
• Several methods have been proposed for mitigating experimental errors for VQE-type simulations [192,194,195]. Such methods enable reliable estimation of expected values of observables on a given trial state without introducing any overhead in terms of extra qubits or quantum gates. Generalizing such error mitigation methods to QPE is a challenging open problem since QPE performs a non-trivial postprocessing of the measurement outcomes that goes beyond computing mean values.
• Classical post-processing methods that enable simultaneous estimation of multiple eigenvalues are highly desirable [204].
• Finally, a natural question is whether the time evolution operator e −iĤt in QPE can be replaced by some other functions ofĤ that are easier to implement [207,208].
H. Quantum algorithms for time evolution

Hamiltonian simulation problem
It was recognized early on [7,209] that a quantum computer can be programmed to efficiently simulate the unitary time evolution of almost any physically realistic quantum system. The time evolution of a quantum system initialized in a given state |ψ(0) is governed by the Schrödinger equation whereĤ is the system's Hamiltonian. Since any fermionic or spin system can be mapped to qubits, see III B, below we assume thatĤ describes a system of n qubits. By integrating Eq. (17) for a time-independent Hamiltonian one obtains the time-evolved state A quantum algorithm for Hamiltonian simulation takes as input a description ofĤ, the evolution time t, and outputs a quantum circuitÛ that approximates the time evolution operator e −itĤ within a specified precision , that is, More generally, the circuitÛ may use some ancillary qubits initialized in the |0 state. The simulation cost is usually quantified by the runtime of the algorithm (the gate count ofÛ ) and the total number of qubits. Applying the circuitÛ to the initial state |ψ(0) provides an -approximation to the time-evolved state |ψ(t) . The final stateÛ |ψ(0) can now be measured to access dynamical properties of the system such as timedependent correlation functions. The time evolution circuitÛ is usually invoked as a subroutine in a larger enveloping algorithm. For example, the quantum phase estimation method employs a controlled version ofÛ to measure the phase accumulated during the time evolution, see III G. The enveloping algorithm is also responsible for preparing the initial state |ψ(0) . While practical applications are concerned with specific Hamiltonian instances, quantum simulation algorithms apply to general classes of Hamiltonians satisfying mild technical conditions that enable a quantum algorithm to access the Hamiltonian efficiently. For example, a Hamiltonian can be specified as a linear combination of elementary interaction terms denotedV 1 , . . . ,V L such that Here α i are real coefficients and V i is the operator norm (the maximum magnitude eigenvalue). In the case of local Hamiltonians [4,209], each termV i acts non-trivially only a few qubits. This includes an important special case of lattice Hamiltonians where qubits are located at sites of a regular lattice and the interactionsV i couple small subsets of nearest-neighbor qubits. Molecular electronic Hamiltonians mapped to qubits assume the form Eq. (20), whereV i are tensor products of Pauli operators which may have a superconstant weight. This situation is captured by the Linear Combination of Unitaries (LCU) model [210]. It assumes that each termV i is a black-box unitary operator that can be implemented at a unit cost by querying an oracle (more precisely, one needs a "select-V " oracle implementing a controlled version ofV 1 , . . . ,V m ). The LCU model also assumes an oracle access to the coefficients α i , see [210] for details. Alternatively, a quantum algorithm can access the Hamiltonian through a subroutine that computes its matrix elements. A HamiltonianĤ is said to be d-sparse [211] if the matrix ofĤ in the standard n-qubit basis has at most d non-zero entries in a single row or column. The sparse Hamiltonian model assumes that positions and values of the nonzero entries can be accessed by querying suitable oracles [212]. Most physically realistic quantum systems can be mapped to either local, LCU, or sparse qubit Hamiltonians such that the corresponding oracles are realized by a short quantum circuit. As described below, the runtime of quantum simulation algorithms is controlled by a dimensionless parameter T proportional to the product of the evolution time t and a suitable norm of the Hamiltonian. One can view T as an effective evolution time.
A formal definition of T for various Hamiltonian models is as follows.
Here H max denotes the maximum magnitude of a matrix element.
It is strongly believed that the Hamiltonian simulation problem is hard for classical computers. For example, Ref. [213] showed that any problem solvable on a quantum computer can be reduced to solving an instance of a suitably formalized Hamiltonian simulation problem with a local Hamiltonian. Technically speaking, the problem is BQP-complete [4]. All known classical methods capable of simulating general quantum systems of the above form require resources (time and memory) exponential in n. On the other hand, Feynman's original insight [7] was that a quantum computer should be capable of simulating many-body quantum dynamics efficiently, such that the simulation runtime grows only polynomially with the system size n and the evolution time T . This intuition was confirmed by Lloyd who gave the first quantum algorithm for simulating local Hamiltonians [209]. The algorithm exploits the Trotter-Suzuki product formula By choosing a sufficiently small Trotter step δ one can approximate the evolution operator e −iHt by a product of fewqubit operators describing evolution under individual interaction terms. Each few-qubit operator can be easily implemented by a short quantum circuit. The runtime of Lloyd's algorithm scales as [209,214] where T is the effective evolution time for local Hamiltonians, see Eq. (21). Importantly, the runtime scales polynomially with all relevant parameters. Lloyd's algorithm was a breakthrough result demonstrating that quantum computers can indeed provide an exponential speedup over the best known classical algorithms for the task of simulating time evolution of quantum systems. However, it was quickly realized that the runtime of Lloyd's algorithm is unlikely to be optimal. Indeed, since any physical system simulates its own dynamics in a real time, one should expect that a universal quantum simulator can attain a runtime scaling only linearly with t. Moreover, for any realistic Hamiltonian composed of short-range interactions on a regular lattice, one should expect that the simulation runtime is linear in the space-time volume nt. Clearly, the scaling Eq. (23) falls far behind these expectations.

I. Algorithmic tools
The last decade has witnessed several improvements in the runtime scaling based on development of new algorithmic tools for Hamiltonian simulation. Most notably, a breakthrough work by Berry et al. [210,212,215] achieved an exponential speedup over Lloyd's algorithm with respect to the precision . A powerful algorithmic tool introduced in [212] is the so-called LCU lemma [212]. It shows how to construct a quantum circuit that implements an operatorÛ = M i=1 β iÛi , where β i are complex coefficients andÛ i are black-box unitary operators. Assuming that U is close to a unitary operator, the lemma shows that U can be well approximated by a quantum circuit of size roughly M i |β i | using roughly i |β i | queries to the oracle implementinĝ U 1 , . . . ,Û M (and their inverses). The simulation algorithm of Ref. [210] works by splitting the evolution into small intervals of length τ and using the truncated Taylor series approximation e −iτĤ ≈ K m=0 (−iτĤ) m /m! ≡Û τ . Accordingly, e −itĤ = (e −iτĤ ) t/τ ≈ (Û τ ) t/τ . Substituting the LCU decomposition of H into the Taylor series one obtains an LCU decomposition ofÛ τ . For a suitable choice of the truncation order K, the truncated seriesÛ τ is close to a unitary operator. ThusÛ τ can be well approximated by a quantum circuit using the LCU lemma. The runtime of this simulation algorithm, measured by the number of queries to the Hamiltonian oracles, scales as [210] t run = T log (T / ) log log (T / ) .
Here T is the effective evolution time for the LCU Hamiltonian model, see Eq. (21). This constitutes a square-root speedup with respect to T and an exponential speedup with respect to the precision compared with Lloyd's algorithm. An important algorithmic tool proposed by Childs [216] is converting a Hamiltonian into a quantum walk. The latter is a unitary operator W that resembles the evolution operator e −itĤ with a unit evolution time t. For a suitable normalization ofĤ, the quantum walk operatorŴ has eigenvalues e ±i arcsin(Eα) , where E α are eigenvalues ofĤ. The corresponding eigenvectors ofŴ are simply related to those ofĤ. Unlike the true evolution operator, the quantum walkŴ can be easily implemented using only a few queries to the oracles describing the HamiltonianĤ, e.g. using the LCU or the sparse models. To correct the discrepancy betweenŴ and the true evolution operator, Low and Chuang [130] proposed the Quantum Signal Processing (QSP) method. One can view QSP as a compiling algorithm that takes as input a black-box unitary operator W , a function f : C → C, and outputs a quantum circuit that realizes f (Ŵ ). Here it is understood that f (Ŵ ) has the same eigenvectors asŴ while each eigenvalue z is mapped to f (z). The circuit realizing f (Ŵ ) is expressed using controlled-Ŵ gates and single-qubit gates on the control qubit. Remarkably, it can be shown that the Low and Chuang algorithm is optimal for the sparse Hamiltonian model [130]. Its runtime, measured by the number of queries to the Hamiltonian oracles, scales as t run = T + log(1/ ) log log(1/ ) (25) where T is the effective evolution time for the sparse model, see Eq. (21). This scaling is optimal in the sense that it matches previously known lower bounds [212,215]. We note that simulation methods based on the sparse and LCU Hamiltonian models have been recently unified using a powerful framework known as qubitization [217]. It provides a general recipe for converting a Hamiltonian into a quantum walk using yet another oracular representation of a Hamiltonian known as a block encoding [217]. Algorithms based on the quantum walk (such as the QSP) or truncated Taylor series may not be the best choice for nearterm applications since they require many ancillary qubits.
In contrast, the original Lloyd algorithm [209] and its generalizations based on higher order product formulas [218] require only as many qubits as needed to express the Hamiltonian. In addition, such algorithms are well-suited for simulating lattice Hamiltonians where qubits are located at sites of a regular D-dimensional grid and each elementary interaction V i couples a few qubits located nearby. Lattice Hamiltonians contain L = O(n) elementary interactions. Each few-qubit operator that appears in a product formula approximating e −itĤ can be expressed using a few gates that couple nearest-neighbor qubits. The corresponding quantum circuit can be easily implemented on a device whose qubit connectivity graph is a D-dimensional grid.
Quite recently, Childs and Su [219] revisited simulation algorithms based on product formulas and demonstrated that their performance is better than what one could expect from naive error bounds. More precisely, an order-p product formula approximates the evolution operator e −iHt under a HamiltonianĤ =Â+B with a simpler operator that involves time evolutions under Hamiltonians A and B such that the approximation error scales as t p+1 in the limit t → 0. Childs and Su [219] showed that the gate complexity of simulating a lattice Hamiltonian using order-p product formulas scales as (nT ) 1+1/p −1/p which shaves off a factor of n from the best previously known bound. Here T is the effective evolution time for the local Hamiltonian model.
A major breakthrough in simulating lattice Hamiltonians has been recently achieved by Haah et al. [221]. This work introduced a new class of product formulas based on the divideand-conquer strategy and showed that lattice Hamiltonians can be simulated with gate complexityÕ(nT ), whereÕ hides logarithmic corrections.
This result confirms the physical intuition that the cost of simulating lattice Hamiltonians scales linearly with the spacetime volume. The algorithm of Ref. [221], sketched in 17, approximates the full evolution operator by dividing the lattice into small (overlapping) subsystems comprising O(log n) qubits each and simulating time evolution of the individual subsystems.
The errors introduced by truncating the Hamiltonian near the boundaries are canceled by alternating between forward FIG. 17. Trotter-like time evolution (top) and divide-and-conquer strategy for lattice systems (bottom). From Ref. [220] and Ref. [221] respectively. and backward time evolutions. The error analysis is based on a skillful application of the Lieb-Robinson bound [222,223] that controls how fast information can propagate across the system during the time evolution.
Finally, we note that while product formulas achieve a better scaling with system size than LCU or QSP methods by exploiting commutativity of Hamiltonian terms, they suffer from worse scaling with the simulation time and the error tolerance. A recent approach of multiproduct formulas [224] combines the best features of both. From Trotter methods, it inherits the simplicity, low-space requirements of its circuits, and a good scaling with system size. From LCU, it inherits the optimal scaling with time and error, up to logarithmic factors. Essentially, the work [224] shows how a certain type of highorder product formula can be implemented with a polynomial gate cost in the order (scaling as p 2 ). In contrast, the standard Trotter-Suzuki formulas scale exponentially with the order (scaling as 5 p ).

J. Open problems
The existing simulation methods such as QSP are optimal in terms of the query complexity. However their runtime may or may not be optimal if one accounts for the cost of implementing the Hamiltonian oracles by a quantum circuit. Indeed, in the case of lattice Hamiltonians, system-specific simulation strategies that are not based on oracular models are known to achieve better runtime [221]. Improved system-specific simulation methods are also available for quantum chemistry systems. For example Ref. [141] achieved a quadratic reduction in the number of interaction terms present in molecular Hamiltonians by treating the kinetic energy and the potential/Coulomb energy operators using two different sets of basis functions -the plane wave basis and its dual. The two bases are related by the fermionic version of the Fourier transform which admits a simple quantum circuit [141]. This simplification was shown to reduce the depth of simulation circuits (parallel runtime) based on the Trotter-Suzuki and LCU decompositions [141]. One may anticipate that further improvements can be made by exploiting system-specific information. It is an interesting open question whether molecular Hamiltonians can be simulated in depth scaling polylogarithmically with the size of the electronic basis. More recent advanced sparse Hamiltonian simulation algorithms can exploit prior knowledge of other Hamiltonian norms to get better performance. For instance, the algorithm of Ref. [225] has an effective time T = d H max H 1 , which is a tighter bound than the one displayed in Eq. (21). This can be further improved to T = t √ d H 1→2 which is a tighter bound than both, see [226] for details. However, these scaling improvements come with a larger constant factor that has not been thoroughly characterized.
In many situations the Hamiltonian to be simulated can be written as H = H 0 + V , where the norm of H 0 is much larger than the norm of V , while the time evolution generated by H 0 can be "fast-forwarded" by calculating the evolution operator analytically. For example, H 0 could represent the kinetic en-ergy in the momentum basis or, alternatively, potential energy in the position basis. Hamiltonian simulation in the interaction picture [227] is an algorithmic tool proposed to take advantage of such situations. It allows one to pay a logarithmic cost (instead of the usual linear cost) with respect to the norm of H 0 . This tool has been used to simulate chemistry in the plane wave basis with O(N 2 ) gates in second quantization, where N is the number of plane waves [227]. In the firstquantization, the gate cost scales as O(N 1/3 K 8/3 ), where K is the number of electrons [228].
Likewise, better bounds on the runtime may be obtained by exploiting the structure of the initial state ψ(0). In many applications, ψ(0) is the Hartree-Fock approximation to the true ground state. As such, the energy of ψ(0) tends to be small compared with the full energy scale ofĤ and one may expect that the time evolution is confined to the low-energy subspace ofĤ. How to develop state-specific simulation methods and improved bounds on the runtime is an intriguing open problem posed in [229].
An interesting alternative to quantum simulations proposed by Poulin et al. [207] is using the quantum walk operator W described above directly in the quantum phase estimation method to estimate eigenvalues ofĤ. This circumvents errors introduced by the Trotter or LCU decomposition reducing the total gate count. It remains to be seen whether other physically relevant quantities such as time dependent correlation functions can be extracted from the quantum walk operator sidestepping the unitary evolution which is relatively expensive compared with the quantum walk.
While the asymptotic runtime scaling is of great theoretical interest, practical applications are mostly concerned with specific Hamiltonian instances and constant approximation error (e.g. four digits of precision). To assess the practicality of quantum algorithms for specific problem instances and compare their runtime with that of state-of-the-art classical algorithms, one has to examine compiling methods that transform a high-level description of a quantum algorithm into a quantum circuit that can be run on particular hardware. This motivates development of compilers tailored to quantum simulation circuits. For example, Hastings et al. [230] examined Trotter-Suzuki type simulation of chemical Hamiltonians mapped to qubits using the Jordan-Wigner transformation. An improved compiling method was proposed reducing the runtime by a factor O(N ), where N is the number of orbitals. Childs et al. [231] and Campbell [214] recently proposed a randomized compiler tailored to Trotter-Suzuki simulations. It improves the asymptotic runtime scaling Eq. (23) to L 5/2 T 3/2 −1/2 and achieves almost 10 5 speedup for simulation of small molecule quantum chemistry Hamiltonians [214,231]. Different types of compilers may be needed for NISQ devices and fault-tolerant quantum machines [232]. Indeed, in the context of NISQ devices, arbitrary single-qubit gates are cheap and the simulation cost is dominated by the number of two-qubit gates. In contrast, the cost of implementing an error-corrected quantum circuit is usually dominated by the number of non-Clifford gates such as the T-gate or the CCZ-gate. Since large-scale quantum simulations are expected to require error correction, developing compiling al-gorithms minimizing the T-gate count is vital. For example, Low et al. [233] recently achieved a square-root reduction in the T-count for a state preparation subroutine employed by the LCU simulation method. Compiling algorithms that minimize the number of generic single-qubit rotations (which are expensive to implement fault-tolerantly) were investigated by Poulin et al. [207]. This work considered implementation of the quantum walk operator W associated with a lattice Hamiltonian that contains only a few distinct parameters (e.g. translation invariant models).
Estimating resources required to solve practically important problems is linked with the study of space-time tradeoffs in quantum simulation [234][235][236]. For example, in certain situations the circuit depth (parallel runtime) can be reduced at the cost of introducing ancillary qubits [237] and/or using intermediate measurements and feedback [238]. Of particular interest are methods for reducing the size of quantum simulation circuits using "dirty" ancilla -qubits whose initial state is unknown and which must be restored to their original form upon the completion of the algorithm [233,237,239], see also [240]. For example, if a Hamiltonian simulation circuit is invoked as a subroutine from a larger enveloping algorithm, the role of dirty ancillas may be played by data qubits borrowed from different parts of the algorithm.
The Hamiltonian simulation problem has a purely classical version. It deals with the ODEṗ = −∂H/∂q,q = ∂H/∂p, where p, q are canonical coordinates of a classical Hamiltonian system. Numerical algorithms for integrating classical Hamiltonian dynamics known as symplectic integrators have a long history [241,242] and are widely applied in simulations of molecular dynamics [243]. One may ask whether new advances in quantum or classical Hamiltonian simulation algorithms can be made by an exchange of ideas between these two fields.
Finally, the problem of simulating Hamiltonians in the presence of noise, without active error correction in the NISQ era, is an open problem. The primary question in this context is whether noise tolerant methods, analogous to variational algorithms in ground-state energy optimization, can also be developed for Hamiltonian simulation.

K. Finite-temperature algorithms
How a quantum computer can be used to simulate experiments on quantum systems in thermal equilibrium is an important problem in the field of quantum simulation. Early algorithms for the simulation of Gibbs states [244][245][246] were based on the idea of coupling the system of interest to a set of ancillary qubits and letting the system and bath together evolve under a joint Hamiltonian, thus mimicking the physical process of thermalization. The main disadvantages of thermalization-based methods are the presence of additional ancillary qubits defining the bath states, and the need to perform time evolution underĤ for a thermalization time t that could be very long.
More recent proposals have focused on ways to generate finite-temperature observables without long system-bath ther-malization times and/or large ancillae bath representations. For example, Ref. [247] showed how to realize the imaginary time evolution operator e −βĤ/2 using Hamiltonian simulation techniques. Applying a suitable version of the Hubbard-Stratonovich transformation the authors of Ref. [247] obtained a representation e −βĤ/2 = α c α U α , where c α are real coefficients and U α are unitary operators describing the time evolution under a HamiltonianĤ 1/2 . Although the square-rootĤ 1/2 is generally not easily available, Ref. [247] showed how to realize it for Hamiltonians composed of fewqubit positive semidefinite terms using an ancillary system. Applying state-of-the-art methods to simulate time evolution underĤ 1/2 and the LCU Lemma (see Section III I) to realize the desired linear combination of unitaries, Ref. [247] obtained a quantum algorithm for preparing the thermal Gibbs state with gate complexity β 1/2 2 n/2 Z −1/2 , where Z is the quantum partition function. Here we ignored a prefactor scaling poly-logarithmically with β and the inverse error tolerance. A closely related but slightly less efficient algorithm was discussed in Ref. [248].
Several more heuristic quantum algorithms for finitetemperature simulations have been proposed recently. However, most of these algorithms are challenging to analyze mathematically and generally do not have performance guarantees. For example, variational ansatz states for the Gibbs state that can be prepared with simple circuits have been proposed to bypass possibly long thermalization times. One example is the product spectrum ansatz [249] (PSA), where a shallow unitary circuit applied to a product thermal state is chosen to minimize the free energy of the system. A different avenue is to sample from the Gibbs state rather than generate it explicitly on the quantum computer. For example, quantum Metropolis sampling [250] samples from the Gibbs state in an analog of classical Metropolis sampling, using phase estimation on a random unitary applied to the physical qubits to "propose" moves, and an iterative amplification procedure to implement the "rejection". Much like the classical Metropolis algorithm, the fixed point of this procedure samples the Gibbs state. Alternatively, the quantum minimally entangled typical thermal state (METTS) algorithm [164] samples from the Gibbs state using imaginary time evolution applied to pure states, implemented via the quantum imaginary time evolution algorithm. One strength of quantum METTS is that it uses only the physical qubits of the system and potentially shallow circuits, thus making it feasible even in the NISQ era, where it has been demonstrated on quantum hardware for small spin systems. An application to the Heisenberg model is shown in 18.
Another alternative is to work within the microcanonical ensemble. This is the basis of the minimal effective Gibbs ansatz (MEGA) [251], which attempts to generate pure states within the energy window corresponding to a microcanonical ensemble (for example, using phase estimation). The basic challenge is to ensure that the energy window is chosen according to the desired temperature. The MEGA method estimates the temperature from asymptotic properties of the ratio of the greater and lesser Green's functions, which in principle FIG. 18. Application of the quantum METTS algorithm to the 6site 1D Heisenberg model. Energy is shown as a function of inverse temperature β, using unitaries acting on D = 2 . . . 5 spins. From [164].
can be measured on the quantum device using the techniques in V.
While there are many different proposed techniques for estimating observables of thermal states which all appear quite plausible on theoretical grounds, little is known about their heuristic performance, and almost none have been tested on real devices. In this sense, thermal state algorithms lag greatly behind those for ground-states for problems of interest in chemistry and materials science. To identify the best way forward, heuristic benchmarking for systems of relevance in physics and chemistry will be of major importance. The problem of benchmarking is further discussed in IV.
L. Hybrid quantum-classical methods

Quantum embedding methods
Embedding algorithms use a divide-and-conquer strategy to break a large quantum simulation into smaller pieces that are more amenable to simulation. The properties of the original model and the reduced models are related to each other in an (ideally) self-consistent fashion. These methods are popular both in condensed matter physics to study correlated electronic materials, where they reduce the problem of solving a bulk fermionic lattice model to that of studying a simpler Anderson-like impurity model, as well as in molecular applications, to reduce the computational scaling of methods in the simulation of complicated molecules.
In classical simulations, there are many flavors of quantum embedding. These can be grouped roughly by the choice of variable used to communicate between regions; dynamical mean-field theory (DMFT) works with the Green's function and self-energy [252][253][254]; density matrix embedding theory with the one-body density matrix [254,255]; density functional embedding via the electron density [254,256], and other methods, such as QM/MM, ONIOM, and fragment MO methods communicate via the electrostatic potential [257].
There has been growing interest in quantum embedding methods in the quantum information community, with the quantum computer playing the role of the quantum me- FIG. 19. Illustration of hybrid quantum-classical embedding for a quantum chemistry system. Fragmentation of the original problem is performed on a classical computing device, and the more complex task of simulating each subproblem is handled by the quantum algorithm. From Ref. [258]. chanical solver for the fragment/impurity problem (see 19). Ref. [220] suggested that a small quantum computer with a few hundred qubits could potentially speed up material simulations based on the DMFT method, and proposed a quantum algorithm for computing the Green's function of a quantum impurity model. Kreula et al. [259] subsequently proposed a proof-of-principle demonstration of this algorithm. Similarly, Rubin [260] and Yamazaki et al [258] have explored the potential of DMET in conjunction with a quantum computer for both condensed-phase lattice models as well as for large molecular calculations. Ground states of quantum impurity models and their structural properties in the context of quantum algorithms have been analyzed in Ref. [261].

Other hybrid quantum-classical algorithms
Beyond standard quantum embedding, there are many other possibilites for hybrid quantum-classical algorithms. Variational quantum eigensolvers have previously been discussed for eigenstate problems, see III E. One can also classically postprocess a quantum simulation of an eigenstate to improve it, in an analog of post-Hartree-Fock and post-completeactive-space methods, as has been explored in [262]. In quantum molecular dynamics, it is natural to use a quantum computer to propagate the wavefunction or density matrix subject to motion of the nuclei, as is done today in classical Born-Oppenheimer dynamics. Similarly, the use of quantum optimizers and quantum annealers, for example to assist classical conformational search [263], can also be viewed as types of quantum-classical hybrids.

M. Open questions
It is clear that one will rely on hybrid quantum-classical algorithms for many years to come, and there remain many open questions. One is how to best adapt quantum algorithms within existing quantum-classical frameworks. For example, Green's function embedding methods are generally formulated in terms of actions rather than Hamiltonians; unfolding into a bath representation, consuming additional qubits, is currently required. More compact representations of the retarded interactions suited for quantum simulation should be explored (see for example, Ref. [264] for a related proposal to generate effective long-range interactions). Similarly one should explore the best way to evaluate Green's functions or density matrices, minimizing the coherence time and number of measurements. It is also possible that new kinds of embedding frameworks should be considered. For example, quantumquantum embedding algorithms within the circuit model of quantum computation have been proposed to simulate largescale "clustered" quantum circuits on a small quantum computer [265,266]. These are circuits that can be divided into small clusters such that there are only a few entangling gates connecting different clusters. Another promising class of embedding algorithms known as holographic quantum simulators was also recently proposed in [267,268]. Such algorithms enable the simulation of 2D lattice models on a 1D quantum computer by converting one spatial dimension into time. Whether it is useful to incorporate a classical component into such quantum-quantum frameworks clearly needs to be explored.
Another question is how best to implement the feedback between the quantum and classical parts of the algorithm. For example, as already discussed in the variational quantum eigensolver, such optimizations require the evaluation of approximate gradients on the quantum device, and noisy optimization with limited gradient information on the classical device. Improving both aspects is clearly needed.

IV. BENCHMARK SYSTEMS
Researchers working on quantum simulation algorithms would greatly benefit from having access to well-defined benchmarks. Such benchmarks help the community by defining common conventions (e.g. choosing specific bases) and curating the best results.
There are two types of benchmarks to develop.
• the first is a benchmark that allows quantum algorithms (possibly on different hardware) to be compared against each other. For example, such problems could include a test suite of molecular Hamiltonian simulation problems for some specific choice of electron basis functions, fermion-to-qubit mapping, evolution time, and the desired approximation error.
• The second is curated data from the best classical methods for specific problems and well-defined Hamiltonians. Wherever possible, the data should not only include ground-state energies, but also excited-states and other observables, and if exact results are not available, an estimate of the precision should be given.
In the near-term era, suitable candidates for benchmarking may include molecular or material science Hamiltonians that can be expressed with about 50 or fewer qubits. Some promising candidates discussed in the literature include lattice spin Hamiltonians [269], and models of correlated electrons such as the 2D Fermi-Hubbard model [270], the uniform electron gas (jellium) [141], and the Haldane pseudo-potential Hamiltonian that models FQHE systems [271][272][273]. While some aspects of these models are easy to solve classically for systems of moderate size, others remain difficult, providing room for quantum advantage.
Naturally, there is a wide range of molecular or materials problems that could be chosen, a small number of which are highlighted in section II. One relevant factor to check when constructing a Hamiltonian benchmark problem is to verify that it is indeed difficult to simulate classically [274]. Ideally, there should be a way to maximally tune the "complexity" for classical simulation, which then defines a natural setting for demonstrating quantum supremacy. Commonly, a way to tune the Hamiltonian (in model Hamiltonians) is to change the parameters of the Hamiltonian directly. In more realistic settings, one may change the size or geometry of the system or the chemical identity of the atoms.

A. Equal-time measurements
It is of course an essential part of any quantum computation to perform a measurement on the final state and thus read out the result of the computation. Conventionally, this is achieved through projective measurements of individual qubits in the computational basis. More complicated operators can be measured through standard techniques, for example by first applying a unitary rotation to the state that maps the operator onto Pauli-Z, or by using an ancillary qubit.
For very complicated operators, however, this can become quite resource-intensive. Consider for example measuring the expectation value of the Hamiltonian, which is required e.g. in the VQE (see III E). This can be done by writing the Hamiltonian as a sum of products of Pauli operators and measuring each one individually. Each measurement must be repeated a sufficient number of times to collect accurate statistics. Since for many applications, the number of operators grows quite quickly with the number of qubits (for example as N 4 in typical quantum chemistry applications) and the state may have to be prepared anew after a projective measurement, it is important to organize the terms in such a way that the number of operations to achieve a desired accuracy is minimized. Some work in this direction appears in Sec. V of Ref. [132], as well as more recently in Refs. [275,276]. Nonetheless, the need to converge a large number of measurements to high precision presents a practical problem in many applications, particularly in many hybrid quantum-classical algorithms. FIG. 20. Quantum circuit to measure Ψ|Â † i (t)Bj|Ψ . HereÂi, Bj are unitaries,Â † i (t) =Û (t)Â † iÛ † (t) is the time evolution ofÂ † i underÛ (t) = e −itĤ , and 2σ− = 2|0 1| = X + iY . Adapted from Ref. [277].

B. Dynamical properties and Green's functions
Much of the experimentally relevant information about a system, for example as obtained in scattering experiments such as optical or X-ray spectroscopy or neutron scattering, is encoded in dynamical properties. Access to these properties allows for a more direct comparison between theoretical predictions and experiments, thus allowing to infer microscopic information that is difficult to extract solely from the experiment.
A convenient way to capture this information is via the single-or few-particle Green's functions, which can be simply related to time-correlation functions of observables, such as the dynamic structure factors and dipole-dipole correlation functions. For example, the particle and hole components of the single-particle Green's function in real time (t ≥ 0) are given by G p αβ (t) = ψ|ĉ α (t)ĉ † β (0)|ψ , G h αβ (t) = ψ|ĉ α (t) †ĉ β (0)|ψ , where α, β can be spin, orbital or site indices, |ψ is the quantum state of interest (for example the ground state), and c ( †) α (t) = e itĤĉ ( †) α e −itĤ . These can be measured by decomposing the fermion creation and annihilation operators into unitary combinations and applying standard techniques, see e.g. Ref. [220,277] (a quantum circuit for measuring products of unitaries, one of which time-evolved, is shown in 20. However, this procedure is expensive as it requires separate calculations for each time t, the particle/hole, real/imaginary and potentially orbital/spatial components of the Green's function, and finally a potentially large number of repetitions for each to achieve some desired target accuracy. While some initial work has shown how to improve this slightly by avoiding having to re-prepare the initial state every time (which could be prohibitive), further improvements and developments of alternatives, would be very valuable.
It should be noted that the behavior of the Green's function is very constrained both at short times (due to sum rules) and at long times (where the decay is governed by the longest time-scale of the system). Therefore, of primary interest is the regime of intermediate times where classical methods are most difficult to apply. However, in these other regimes, the additional structure could potentially be employed to reduce the computational effort also on a quantum computer.

VI. CONCLUSIONS
In this review we surveyed some of the quantum problems of chemistry and materials science that pose a challenge for classical computation, and the wide range of quantum algorithms that can potentially target such problems. The field of quantum algorithms for physical simulation is growing at an exponential pace, thus this review only provides a snapshot of progress. However, there are many shared ideas and components and we have tried to bring out these common features, as well as highlight the many open problems that remain to be explored.
Although very few experiments have yet been carried out in actual quantum hardware, it will soon become possible to test many proposals in real quantum computational platforms. While it is not at all certain then which techniques will perform best in practice, what is clear is that such experiments will greatly enhance our understanding of the heuristics of quantum algorithms. Indeed, the widespread availability of classical computers led to a new wave of algorithmic innovation, and the establishment of computational chemistry, materials science, and physics as scientific fields. Should quantum hardware continue to develop at the current rate, we can thus look forward also to a new wave of developments in quantum algorithms, and the emergence of new computational disciplines dedicated to the study of chemistry, materials science, and physics on quantum computers.

ACKNOWLEDGEMENT
This review was adapted and condensed from a report produced for a US National Science Foundation workshop held between January 22-24 2019 in Alexandria VA. Funding for the workshop and for the production of the initial report was provided by the US National Science Foundation via award no. CHE 1909531.