# Fault-Tolerant Quantum Simulations of Chemistry in First Quantization

Yuan Su  corresponding author<EMAIL_ADDRESS>Google Quantum AI, Venice,
CA 90291, United States  Institute for Quantum Information and Matter,
Caltech, Pasadena, CA 91125, United States  Dominic W. Berry  corresponding
author<EMAIL_ADDRESS>Department of Physics and Astronomy, Macquarie
University, Sydney, NSW 2109, Australia  Nathan Wiebe  corresponding author:
<EMAIL_ADDRESS>Department of Computer Science, University of Toronto,
ON M5S 2E4, Canada  Pacific Northwest National Laboratory, Richland, WA 99354,
United States  Google Quantum AI, Venice, CA 90291, United States  Nicholas
Rubin  Google Quantum AI, Venice, CA 90291, United States  Ryan Babbush
corresponding author<EMAIL_ADDRESS>Google Quantum AI, Venice, CA
90291, United States

###### Abstract

Quantum simulations of chemistry in first quantization offer some important
advantages over approaches in second quantization including faster convergence
to the continuum limit and the opportunity for practical simulations outside
of the Born-Oppenheimer approximation. However, since all prior work on
quantum simulation of chemistry in first quantization has been limited to
asymptotic analysis, it has been impossible to directly compare the resources
required for these approaches to those required for the more commonly studied
algorithms in second quantization. Here, we compile, optimize and analyze the
finite resources required to implement two first quantized quantum algorithms
for chemistry from Babbush _et al._ Babbush  _et al._ (2019) that realize
block encodings for the qubitization and interaction picture frameworks of Low
_et al._ Low and Chuang (2019); Low and Wiebe (2018). The two algorithms we
study enable simulation with gate complexities of $\widetilde{\cal
O}(\eta^{8/3}N^{1/3}t+\eta^{4/3}N^{2/3}t)$ and $\widetilde{\cal
O}(\eta^{8/3}N^{1/3}t)$ where $\eta$ is the number of electrons, $N$ is the
number of plane wave basis functions, and $t$ is the duration of time-
evolution ($t$ is linearly inverse to target precision when the goal is to
estimate energies). In addition to providing the first explicit circuits and
constant factors for any first quantized simulation, and then introducing
improvements which reduce circuit complexity by about a thousandfold over
naive implementations for modest sized systems, we also describe new
algorithms that asymptotically achieve the same scaling in a real space
representation. Finally, we assess the resources required to simulate various
molecules and materials and conclude that the qubitized algorithm will often
be more practical than the interaction picture algorithm. We demonstrate that
our qubitization algorithm often requires much less surface code spacetime
volume for simulating millions of plane waves than the best second quantized
algorithms require for simulating hundreds of Gaussian orbitals.

###### Contents

  1. I Introduction
    1. I.1 Background on representing molecular Hamiltonians in a first quantized plane wave basis
    2. I.2 Overview of the Hamiltonian simulation frameworks deployed in this work
  2. II The qubitization based algorithm
    1. II.1 Circuit implementation of qubitization
    2. II.2 Block-encoding the kinetic operator
    3. II.3 State preparation for $U$ and $V$
    4. II.4 The $\mathrm{SEL}$ operations for $U$ and $V$
    5. II.5 Total cost of constructing the qubitization operator
  3. III The interaction picture based algorithm
    1. III.1 Implementation choices for the interaction picture algorithm
    2. III.2 Exponentiating the kinetic operator
    3. III.3 Implementing $\mathrm{SEL}$ and $\mathrm{PREP}$
    4. III.4 Total cost of the interaction picture simulation
  4. IV Results and discussion
    1. IV.1 The resources required to simulate real systems
    2. IV.2 Conclusion and outlook
  5. A Variable and function names
  6. B Background on Galerkin representations
    1. B.1 Galerkin basis representations in first quantization
    2. B.2 Galerkin basis representations in second quantization
  7. C Qubit costings
    1. C.1 Qubit costings for qubitization approach
    2. C.2 Qubit costings for interaction picture approach
  8. D Selection between $T$, $U$, and $V$
  9. E Detailed costing of preparation of $k$ state
  10. F The modified phase gradient state
  11. G Complexity of binary squaring
  12. H Bound on time-discretization error
  13. I Interaction picture algorithm for simulating generic Hamiltonians
  14. J Bound on the success probability for preparing equal superposition
  15. K Simulating first quantized real space grid representations

##  I Introduction

The last several years have seen monumental developments in quantum
algorithms, resulting in significant progress towards the goal of showing that
a modest sized quantum computer can provide a decisive advantage for
scientifically relevant problems. The expectation that quantum simulation
might enable an exponential speedup for physical problems stretches back to
Feynman’s original proposal for quantum computers Feynman (1982); Lloyd (1996)
and was first concretely proposed for solving the the central problem of
quantum chemistry by Aspuru-Guzik _et al_. Aspuru-Guzik  _et al._ (2005).
Since then, progress has accelerated on multiple fronts: improvements in
general purpose simulation frameworks Berry  _et al._ (2015); Wiebe and
Granade (2015); Low and Chuang (2019); Low and Wiebe (2018); Poulin  _et al._
(2017); Haah  _et al._ (2018); Campbell (2019); Childs  _et al._ (2021); Motta
_et al._ (2020a); Şahinoğlu and Somma (2020); Tran  _et al._ (2021), more
efficient algorithms for quantum chemistry Whitfield  _et al._ (2011); Wecker
_et al._ (2014); Hastings  _et al._ (2015); Poulin  _et al._ (2015); Babbush
_et al._ (2016, 2018a); Berry  _et al._ (2018, 2019); Kivlichan  _et al._
(2018, 2020); McArdle  _et al._ (2020); von Burg  _et al._ (2020); Lee  _et
al._ (2021); Su  _et al._ (2021), more advantageous-to-quantum-simulate
representations of molecular Hamiltonians Babbush  _et al._ (2018b); Seeley
_et al._ (2012); Bravyi  _et al._ (2017); Setia and Whitfield (2018); Motta
_et al._ (2021); Bauman  _et al._ (2019); Takeshita  _et al._ (2020); McClean
_et al._ (2020); Motta  _et al._ (2020b) and an improved understanding of the
problems that genuinely require a quantum computer versus those that can be
classically computed Reiher  _et al._ (2017); Li  _et al._ (2019).

However, a problem lurks within this discussion that has not received as much
attention as it deserves. Specifically, most existing methods for simulating
quantum chemistry (especially those leveraging simple basis functions like
plane waves) become impractical as we scale to the continuum limit. This is
because the second quantized simulations that are deployed in the vast
majority of quantum computing algorithms for chemistry require a number of
qubits that scales linearly with the number of spin-orbital basis functions.
Dating back to early work by Kassal _et al._ Kassal  _et al._ (2008), first
quantized simulations of chemistry have been proposed as a method to
ameliorate this problem. The central idea behind such simulations is to track
the basis state that each particle is in rather than storing the occupancy of
each basis state. This has two significant benefits: the number of qubits
needed to represent the state scales logarithmically with the number of basis
functions (although linearly with the number of particles) and also, such
simulations are easy to adapt to cases where the entanglement between the
electronic and nuclear subsystems is non-negligible. Despite the promise of
these methods, no prior work investigating first quantized approaches for
chemistry has managed to rigorously estimate all of the constant factors
associated with the realization of these algorithms111In the first paper to
discuss error-correcting quantum algorithms for chemistry, Jones _et al._ Cody
Jones  _et al._ (2012) estimate the number of Toffoli gates required to
implement subroutines (Trotter steps) of the Kassal _et al._ Kassal  _et al._
(2008) algorithm. However, their work stops short of the error analysis
required to predict the constant factors associated with performing fixed
precision quantum simulation. For example, they do not estimate how many
Trotter steps would be required to achieve target precision.

Our aim in this paper is to address this shortcoming by providing a full
costing of the two leading quantum simulation algorithms for chemistry in a
first quantized representation, both first described by Babbush _et al._
Babbush  _et al._ (2019). These algorithms comprise asymptotically efficient
schemes for realizing the oracles required by the qubitization and interaction
picture frameworks of Low _et al._ Low and Chuang (2019); Low and Wiebe
(2018). Here, we compile and optimize these algorithms within an error-
correctable gate set. We develop new techniques to reduce the complexity,
resulting in about three orders of magnitude improvement over a naive
implementation for even modest sized systems. Our detailed cost analysis
reveals that despite the improved scaling of the interaction picture based
method, the qubitization based method often proves to be more practical. We
conduct numerical experiments to characterize the overhead of both approaches
and find that they are highly competitive with second quantized approaches,
often requiring significantly fewer resources to reach comparable accuracy.

The layout of this paper is as follows. In the remainder of this section, we
provide the formalism and context needed to understand the two most efficient
schemes for simulating the molecular Hamiltonian in first quantization.
Section II develops and analyzes our optimized implementation of the
qubitization based algorithm in first quantization. Section III contains a
similar analysis for the interaction picture based algorithm. Finally, Section
IV numerically compares these algorithms (to each other and to other methods
in the literature), estimating the resources required to simulate real systems
in a first quantized basis.

###  I.1 Background on representing molecular Hamiltonians in a first
quantized plane wave basis

The accurate simulation of quantum systems was Feynman’s original vision for
quantum computing Feynman (1982). There is perhaps no more natural application
of this idea than the simulation of atoms and molecules, governed by the
interactions of electrons and nuclei. Such systems give rise to the properties
of most matter, determining everything from the rates of chemical reactions to
the conductivity of materials. In the non-relativistic case, the dynamics of
these systems are governed by the Coulomb Hamiltonian:

$\displaystyle
H=\underbrace{-\sum_{i=1}^{\eta}\frac{\nabla^{2}_{i}}{2}}_{T}\underbrace{-\sum_{i=1}^{\eta}\sum_{\ell=1}^{L}\frac{\zeta_{\ell}}{\left\|R_{\ell}-r_{i}\right\|}}_{U}\underbrace{+\sum_{i\neq
j=1}^{\eta}\frac{1}{2\left\|r_{i}-r_{j}\right\|}}_{V}\underbrace{-\sum_{\ell=1}^{L}\frac{\nabla^{2}_{\ell}}{2m_{\ell}}}_{T_{\rm
nuclei}}\underbrace{+\sum_{\ell\neq\kappa=1}^{L}\frac{\zeta_{\ell}\zeta_{\kappa}}{2\left\|R_{\ell}-R_{\kappa}\right\|}}_{V_{\textrm{nuclei}}}\,.$
(1)

Here we have used atomic units where $\hbar$ as well as the mass and charge of
the electron are unity, and $\|\cdot\|$ denotes the two-norm. In the above
expression $i,j$ represent electronic degrees of freedom and $\ell,\kappa$
represent nuclear degrees of freedom; thus, $r_{i}$ represent the positions of
electrons whereas $R_{\ell}$ represent the positions of nuclei, $m_{\ell}$ the
atomic masses of nuclei, and $\zeta_{\ell}$ the atomic numbers of nuclei.
Throughout this work we will use $\eta$ to denote the number of electrons in
our simulation and $L$ to denote the number of nuclei in our simulation. We
keep a catalogue of the different symbols used throughout this paper in
Appendix A.

The methods we analyze and develop in this work are useful for simulations of
the dynamics of (1) (treating both nuclear and electronic degrees of freedom
as fully quantum) with only minor modifications. But in order to simplify and
reduce our analysis, we will primarily focus on simulating (1) under the Born-
Oppenheimer approximation: i.e., $H(R)=T+U(R)+V+C(R)$ where $C(R)$ is a
constant given by the Coulomb repulsion of the point charges of “classical”
nuclei with locations $R_{\ell}$. Thus, the Born-Oppenheimer approximation
assumes that one can decouple nuclear and electronic degrees of freedom,
treating the former as essentially “classical”Polanyi (1972); Born and
Oppenheimer (1927). This is appropriate for many systems near room temperature
because the electronic degrees of freedom are thousands of times faster than
the nuclear degrees of freedom (due to the disparity in mass between electrons
and nuclei) and hence, electrons often relax nearly instantaneously to their
quantum ground states whereas the nuclei are often heavy enough to be modelled
as classical objects responding to forces arising from the electronic
interactions. However, the Born-Oppenheimer approximation is known to break
down in some circumstances, including at low temperatures or when chemical
bonds involve low atomic mass nuclei (e.g., due to the tunneling of hydrogen
atoms); hence, methods for performing simulations of (1), referred to as “non-
Born-Oppenheimer dynamics” are also of great interest Yarkony (2012); Butler
(1998); Curchod and Martínez (2018); Zhu and Yarkony (2016).

The frequently studied “electronic structure problem” is to solve for the
ground state energies of electrons interacting in the external potential of
nuclei under the Born-Oppenheimer approximation (i.e., to solve for the ground
state energies of $H(R)$). As a function of the nuclear coordinates $R$, these
energies define an energy surface and those energy surfaces can be extremely
helpful for understanding the mechanisms and dynamics of chemical reactions as
well as material properties. For example, energy surfaces are commonly sought
after so that one may simulate the “classical” dynamics of nuclei moving on
these “quantum” energy surfaces by integrating Newton’s equations of motion
for the nuclei. When the trajectories of these classical nuclei responding to
quantum electronic energy surfaces are coupled to a finite temperature bath so
that the system may sample the canonical or grand canonical ensemble (rather
than just Hamiltonian dynamics), these simulations are referred to as “ _ab
initio_ molecular dynamics” Car and Parrinello (1985); Marx and Hutter (2009);
Hutter (2012). While the electronic energy surfaces (referred to in the
molecular dynamics community as “force-fields”) are often determined
empirically or by crude calculations that introduce large errors into the end
result, molecular dynamics simulations are among the most common calculations
in scientific computing and can provide invaluable insights into the
thermodynamic properties of real materials and chemical reactions.

In order to simulate molecular systems on a computer one must discretize them
in some fashion. For reasons discussed in more depth in Appendix B, most
molecular modeling employs a Galerkin representation which involves projecting
the system onto some well behaved set of basis functions $\\{\phi_{p}(r)\\}$.
The matrix elements of the Hamiltonian operators are then given by the
following integrals over these basis functions:

$\displaystyle T_{pq}^{(m)}$
$\displaystyle=\int\\!dr\,\phi_{p}^{*}\left(r\right)\left(-\frac{\nabla^{2}}{2m}\right)\phi_{q}\left(r\right)$
(2) $\displaystyle U_{pq}$
$\displaystyle=\sum_{\ell=1}^{L}\int\\!dr\,\phi_{p}^{*}\left(r\right)\left(\frac{\zeta_{\ell}}{\left\|r-R_{\ell}\right\|}\right)\phi_{q}\left(r\right)$
(3) $\displaystyle V_{pqrs}^{(\alpha,\beta)}$
$\displaystyle=\int\\!dr_{1}\,dr_{2}\,\phi_{p}^{*}\left(r_{1}\right)\phi_{q}^{*}\left(r_{2}\right)\left(\frac{\alpha\beta}{\left\|r_{1}-r_{2}\right\|}\right)\phi_{r}(r_{2})\phi_{s}\left(r_{1}\right)\,.$
(4)

We note that simple grid representations (as opposed to discrete-variable-
representation-like grids, e.g., those discussed for quantum simulation in
Babbush  _et al._ (2018b)) are incompatible with Galerkin discretizations. For
example, if one associates a grid with compact basis functions having disjoint
support (e.g., delta functions or step functions) then evaluation of the above
integrals would give both potential operators and kinetic operators that are
simultaneously diagonal. For more discussion on real space representations see
Appendix K.

A key consideration for Galerkin discretizations is the compactness of the
basis set for the states of interest (often the ground state); i.e., how many
basis functions are required to converge the discretization of states of
interest to within some threshold $\epsilon$ of the continuum limit.
Throughout this paper we will use $N$ to refer to the number of basis
functions used to discretize the system. For most reasonable choices of basis
sets, the asymptotic behavior is that $\epsilon={\cal O}(1/N)$ Grüneis  _et
al._ (2013); Hättig  _et al._ (2011); Harl and Kresse (2008); Shepherd  _et
al._ (2012); Helgaker  _et al._ (1998); Klopper (1995); Halkier  _et al._
(1998). However, the constant factors in this scaling can differ considerably
depending on the context of the simulation and the choice of basis. The
primary factor limiting the convergence of these discretizations is the
resolution of “cusps” that are known to be features of eigenstates of
molecular Hamiltonians Kato (1957). Wavefunction cusps appear at all points in
space where particles overlap. Away from these cusps, the wavefunction is
generally smooth and easier to converge. When working within the Born-
Oppenheimer approximation, the nuclei have fixed position and there is a clear
cusp in the electronic charge density at the location of those nuclei. The
electronic density tends to be sharply peaked at these nuclei and then falls
off exponentially away from the nuclei. To capture this structure, one of the
most common types of basis sets used for electronic structure simulations are
nuclei centered Gaussian orbitals. Often, one uses these “primitive” orbitals
to construct even more compact numerical orbitals; e.g., by using the
primitive orbitals to discretize the diagonalization of an approximate single-
particle description of the physics. When the single-particle description is a
mean-field model of molecule (usually Hartree-Fock) these numerically
optimized orbitals are referred to as “molecular orbitals”.

Aside from Gaussian orbitals, the other most commonly used class of basis
functions are plane waves. Plane waves are the type of basis function that we
will focus on in this paper. Plane wave basis functions are eigenstates of the
linear momentum operator (and thus, also the operator $T$), expressed as

$\phi_{p}\left(r\right)=\sqrt{\frac{1}{\Omega}}e^{-\mathrm{i}\,k_{p}\cdot
r}\,,$ (5)

where $r$ is a position vector in real space, $\Omega$ is the computational
cell volume and $k_{p}$ is a reciprocal lattice vector in three dimensions. In
this paper we will focus on methods to perform simulations with plane waves
defined over a cubic reciprocal lattice so that

$k_{p}=\frac{2\pi p}{\Omega^{1/3}}\qquad\qquad p\in G\qquad\qquad
G=\left[-\frac{N^{1/3}-1}{2},\frac{N^{1/3}-1}{2}\right]^{3}\subset\mathbb{Z}^{3}\,.$
(6)

We note that while this reciprocal lattice is appropriate for simulating non-
periodic systems or those with cubic periodicity, to simulate systems with
other crystal symmetries one may need to generalize the above equation to the
case of non-orthogonal Bravais vectors (we leave this to future work) Martin
(2004).

Plane waves are a proper basis (unlike a grid) but, like a grid, provide an
unbiased discretization of space that is especially well suited to
representing dynamics. Due to their regularity, plane waves also give very
systematic convergence to the continuum limit, allowing one to accurately
estimate properties in the continuum by extrapolating from a series of
calculations with progressively more plane waves (the extrapolation results in
higher accuracy predictions than even the largest simulation used in the
extrapolation). But from a quantum algorithms perspective, perhaps the most
useful property of plane waves is that the integrals defining their Galerkin
representation have a convenient closed form (arising from the Fourier
transform of the Coulomb kernel) Martin (2004):

$T_{pq}^{(m)}=\delta_{pq}\frac{\left\|k_{p}\right\|^{2}}{2m}\qquad\qquad
U_{pq}=\frac{4\pi}{\Omega}\sum_{\ell=1}^{L}\zeta_{\ell}\frac{e^{ik_{q-p}\cdot
R_{\ell}}}{\left\lVert k_{p-q}\right\rVert^{2}}\qquad\qquad
V_{pqrs}^{(\alpha,\beta)}=\delta_{p-s,r-q}\frac{4\pi\,\alpha\beta}{\Omega\left\|k_{\nu}\right\|^{2}}\,,$
(7)

where $\nu=p-s=r-q\neq 0$ in the last expression. As we will see, this closed
form enables some especially efficient quantum algorithms.

The structure of the plane wave basis can be used to give quantum algorithms
with lower asymptotic gate complexity compared to algorithms with Gaussian
orbitals. One reason for this is that most methods for simulating second
quantized Hamiltonians tend to have some cost that scales with the number of
terms in the Hamiltonian (or specifically, the number of terms in the two-body
operator since that is the difficult part). A second quantized plane wave
Hamiltonian has ${\cal O}(N^{3})$ terms in the two-body operator (less than
the ${\cal O}(N^{4})$ terms of a molecular orbital representation by a factor
of $N$ due to conservation of linear momentum) but only ${\cal O}(N^{2})$
terms in the dual basis obtained by Fourier transforming the basis. This
argument was first made in the context of second quantized quantum simulations
in Babbush  _et al._ (2018b). Today, the most efficient second quantized
algorithm for simulating plane wave electronic structure Hamiltonians with $N$
plane waves and $\eta$ electrons has Toffoli complexity
$(N^{5/3}/\eta^{2/3}+N^{4/3}\eta^{2/3})N^{o(1)}$ with space complexity ${\cal
O}(N\log N)$ Su  _et al._ (2021), where this complexity assumes that particle
density is held fixed as the system size grows. We compare the scaling of all
prior quantum algorithms for plane wave electronic structure in more detail in
Table 1. By contrast, the best second quantized algorithms for simulating
arbitrary basis (e.g., molecular orbital) electronic structure Hamiltonians
with $N$ orbitals have Toffoli complexity that is roughly $\widetilde{\cal
O}(N^{3})$ with space complexity ${\cal O}(N\log N)$ Lee  _et al._
(2021)222Here and throughout the paper we use $\widetilde{\cal O}(\cdot)$ to
indicate an asymptotic upper bound suppressing polylogarithmic factors and
$o(1)$ to represent a positive number that approaches zero as some parameter
grows.. Thus, the best plane wave algorithms have better asymptotic scaling
than the best arbitrary basis algorithms. However, while the basis set
discretization error of both plane waves and Gaussians is asymptotically
refined as $\epsilon={\cal O}(1/N)$, in practice one often requires far fewer
Gaussian orbitals than plane waves in order to achieve “chemically accurate”
electronic structure simulations333Conventionally, “chemical accuracy” is
defined as having precision in the energy that is better than 1 kcal/mol, or
equivalently, 0.0016 Hartree. An error in the energy on the order of this
quantity corresponds to an error in the predicted rates of chemical reactions
by roughly an order of magnitude at room temperature Helgaker  _et al._
(2000).. For some systems, one could easily require thousands of times more
plane waves to reach chemical accuracy.

Plane waves are especially well suited to treat regular solid-state systems
(like crystals) since the basis is naturally periodic. Furthermore, there are
some special systems, such as the uniform electron gas, for which plane waves
are essentially the most compact basis. Thus, several papers Babbush  _et al._
(2018a); Kivlichan  _et al._ (2020); McArdle  _et al._ (2021) suggest that
using a plane wave basis in second quantization one can encounter interesting
and challenging instances of electronic structure that require fewer resources
to solve than any second quantized Gaussian orbital representation. However,
for most chemical calculations on a quantum computer, a second quantized plane
wave representation is not realistic to use and this is why the vast majority
of work on quantum computing quantum chemistry focuses on using Gaussian
orbitals. Even though second quantized plane wave representations have
asymptotically lower gate complexity than second quantized Gaussian basis
algorithms, the space complexity required for accurate calculations using
second quantized plane waves is daunting. For example, if a calculation
requires one-hundred Gaussian orbitals but one-hundred thousand plane waves
that means in second quantization the former would require hundreds of logical
system qubits (reasonable) whereas the latter would require hundreds of
thousands of logical system qubits (an extravagant cost). Thus, for chemistry
calculations in second quantization, Gaussian orbitals are more practical than
plane waves.

However, in first quantization the story is rather different. Because the
space complexity of first quantized representations is only ${\cal O}(\eta\log
N)$ one can easily simulate systems with a very large number of plane waves
without the number of qubits becoming too extravagant. As we will later
discuss, the gate complexity of first quantized plane wave algorithms can also
have sublinear scaling in $N$; for example, the algorithm of Babbush  _et al._
(2019) has gate complexity $\widetilde{\cal O}(\eta^{8/3}N^{1/3})$. Sublinear
scaling in $N$ would be impossible in second quantization as each of $N$
qubits must be acted upon at least once. In principle one can also perform
simulations of Gaussian orbitals in first quantization but due to the lack of
structure in the integrals this would likely give a relatively inefficient
gate complexity. The work of Babbush  _et al._ (2018c) does something that is
very close to this; in that work, the authors simulate a fixed particle number
manifold of the second quantized molecular orbital Hamiltonian (this is known
in chemistry as the “configuration interaction” representation). Their
approach has space complexity ${\cal O}(\eta\log N)$ but it is technically a
second quantized approach in a fixed particle manifold rather than a first
quantized approach because the symmetries are still encoded by how the
Hamiltonian operator acts on the state rather than the algorithm requiring
anti-symmetrized initial states. However, the gate complexity of that approach
is $\widetilde{\cal O}(\eta^{2}N^{3})$ which is still relatively high in $N$,
meaning that one cannot refine the basis to extremely large sizes the way that
one can with first quantized plane wave approaches.

Year | Reference | Primary innovation | Logical qubits | Toffoli/T complexity
---|---|---|---|---
2017 | Babbush _et al_. Babbush  _et al._ (2018b) | Using plane waves with Trotter | $N+{\cal O}(\log N)$ | $\widetilde{\cal O}(\eta^{2}N^{17/6}\sqrt{1+\eta\Omega^{1/3}/N^{1/3}}/(\Omega^{5/6}\epsilon^{3/2}))$
2017 | Babbush _et al_. Babbush  _et al._ (2018b) | Using plane waves with LCU | $N+{\cal O}(\log N)$ | $\widetilde{\cal O}((N^{4}/\Omega^{1/3}+N^{11/3}/\Omega^{2/3})/\epsilon)$
2018 | Babbush _et al_. Babbush  _et al._ (2018a) | Linear scaling quantum walks | $N+{\cal O}(\log N)$ | $\widetilde{\cal O}((N^{10/3}/\Omega^{1/3}+N^{8/3}/\Omega^{2/3})/\epsilon)$
2018 | Low _et al_. Low and Wiebe (2018) | Using the interaction pic. | ${\cal O}(N\log N)$ | $\widetilde{\cal O}(N^{8/3}/(\Omega^{2/3}\epsilon))$
2018 | Babbush _et al_. Babbush  _et al._ (2019) | First quantized qubitization | ${\cal O}(\eta\log N)$ | $\widetilde{\cal O}((\eta^{3}N^{1/3}/\Omega^{1/3}+\eta^{2}N^{2/3}/\Omega^{2/3})/\epsilon)$
2018 | Babbush _et al_. Babbush  _et al._ (2019) | First quantized interaction pic. | ${\cal O}(\eta\log N)$ | $\widetilde{\cal O}(\eta^{3}N^{1/3}/(\Omega^{1/3}\epsilon))$
2019 | Kivlichan _et al_. Kivlichan  _et al._ (2020) | Better Trotter steps | $N+{\cal O}(\log N)$ | $\widetilde{\cal O}(N^{3}/(\Omega^{2/3}\epsilon^{3/2}))$
2019 | Childs _et al_. Childs  _et al._ (2021) | Tighter Trotter bounds | ${\cal O}(N\log N)$ | $N^{7/3+o(1)}/(\Omega^{1/3}\epsilon^{1+o(1)})$
2020 | Su _et al_. Su  _et al._ (2021) | Yet tighter Trotter bounds | ${\cal O}(N\log N)$ | $(\eta/\Omega^{1/3}+N^{1/3}/\Omega^{2/3})N^{4/3+o(1)}/\epsilon^{1+o(1)}$
Table 1:  Best quantum algorithms for phase estimating chemistry in a plane
wave basis. $N$ is number of basis functions, $\eta<N$ is number of electrons,
$\Omega$ is the computational cell volume, and $\epsilon$ is target precision.
In some of these papers it is assumed that $\Omega\propto N$ or that
$\Omega\propto\eta$, but here we report the complexities without any such
assumptions. Unlike other entries here, the scaling of the Kivlichan _et al._
Kivlichan  _et al._ (2020) algorithm is empirically observed for specific
systems (with some systems scaling somewhat better and some systems scaling
somewhat worse) as opposed to rigorous upper-bounds. The ${\cal O}(N\log N)$
space complexity of Low and Wiebe (2018); Childs  _et al._ (2021); Su  _et
al._ (2021) results from a Fourier transform based method of computing the
potential operator discussed in Low and Wiebe (2018). If those papers were
instead to use the Trotter steps introduced in Kivlichan  _et al._ (2018) or
Kivlichan  _et al._ (2020), the algorithms would have $N+{\cal O}(\log N)$
space complexity but gate complexity that is worse by approximately a factor
of $N$. Algorithms with ${\cal O}(N)$ space complexity are likely impractical
for any real problems in chemistry when using plane waves due to the large
number of plane waves required to reach chemical accuracy; thus, ${\cal
O}(N\log N)$ space complexity is even less viable. We note that Appendix K
introduces new algorithms defined in real space that match the scaling of
those from Babbush _et al._ Babbush  _et al._ (2019).

We have thus argued that plane waves are especially well suited for use within
first quantized quantum simulations whereas Gaussian orbitals are especially
well suited for use within second quantized quantum simulations. The relative
merits of these two representations is a fairly subtle issue because they have
different properties. One requires more plane waves than Gaussian orbitals to
converge a calculation but the first quantized plane wave algorithm
asymptotically scales better for essentially all values of $\eta$ and $N$.
That is because $\eta<N$, so the $\widetilde{\cal O}(\eta^{8/3}N^{1/3})$
Toffoli complexity of the best first quantized plane algorithm Babbush  _et
al._ (2019) reduces to the roughly $\widetilde{\cal O}(N^{3})$ Toffoli
complexity of the best second quantized arbitrary basis algorithm Lee  _et
al._ (2021) in the worst case that $\eta=\Theta(N)$. There is also a
substantial difference in how the algorithms approach the continuum limit. For
example, the Toffoli complexity to achieve basis error $\epsilon$ scales as
$\widetilde{\cal O}(1/\epsilon^{3})$ with $\widetilde{\cal O}(1/\epsilon)$
space complexity for the arbitrary basis algorithm but the Toffoli complexity
is $\widetilde{\cal O}(1/\epsilon^{1/3})$ with space complexity ${\cal
O}(\log(1/\epsilon))$ for the plane wave algorithm. This enormous difference –
a polynomial scaling difference by a power of nine – suggests a strong
preference for the first quantized plane wave algorithm when targeting high
precision444The first quantized plane wave algorithm of Babbush  _et al._
(2019) has a Clifford complexity matching the Toffoli complexity whereas the
second quantized arbitrary basis algorithm of Lee  _et al._ (2021) has
Clifford complexity that is worse than the Toffoli complexity by a factor of
$\widetilde{\cal O}(N)$. Hence, in terms of overall gate complexity the
scaling to the continuum limit is actually $\widetilde{\cal
O}(1/\epsilon^{4})$ versus $\widetilde{\cal O}(1/\epsilon^{1/3})$ – a
difference by a power of twelve!. Furthermore, the plane wave representation
is the clear preference when performing simulations beyond the Born-
Oppenheimer approximation. The reason for this is because the advantage of
Gaussian basis sets lies in their ability to resolve cusps in the wavefunction
centered on the nuclei. But outside of the Born-Oppenheimer approximation the
nuclei are treated explicitly and thus one cannot center Gaussians on the
nuclei because the nuclei do not have a fixed position that is known in
advance. Thus, for simulations of nuclear quantum dynamics plane waves are
going to be the better choice.

By plugging the plane wave integrals given in (7) into the definitions of the
Born-Oppenheimer and non-Born-Oppenheimer first quantized Galerkin
discretizations discussed in Appendix B, we obtain

$\displaystyle H_{\rm BO}$
$\displaystyle=T+U+V+\frac{1}{2}\sum_{\ell\neq\kappa=1}^{L}\frac{\zeta_{\ell}\zeta_{\kappa}}{\left\|R_{\ell}-R_{\kappa}\right\|}$
(8) $\displaystyle H_{\rm non-BO}$ $\displaystyle=T+T_{\rm nuc}+U_{\rm non-
BO}+V+V_{\rm nuc}$ (9) $\displaystyle T$
$\displaystyle=\sum_{i=1}^{\eta}\sum_{p\in
G}\frac{\left\|k_{p}\right\|^{2}}{2}\lvert p\rangle\\!\langle
p\rvert_{i}\qquad\qquad T_{\rm nuc}=\sum_{\ell=\eta+1}^{L+\eta}\sum_{p\in
G}\frac{\left\|k_{p}\right\|^{2}}{2\,m_{\ell}}\lvert p\rangle\\!\langle
p\rvert_{\ell}$ (10) $\displaystyle U$
$\displaystyle=-\frac{4\pi}{\Omega}\sum_{\ell=1}^{L}\sum_{i=1}^{\eta}\sum_{\begin{subarray}{c}p,q\in
G\\\ p\neq q\end{subarray}}\bigg{(}\zeta_{\ell}\frac{e^{ik_{q-p}\cdot
R_{\ell}}}{\left\lVert k_{p-q}\right\rVert^{2}}\bigg{)}\lvert
p\rangle\\!\langle q\rvert_{i},$ (11) $\displaystyle U_{\rm non-BO}$
$\displaystyle=-\frac{4\pi}{\Omega}\sum_{i=1}^{\eta}\sum_{\ell=\eta+1}^{L+\eta}\sum_{p,q\in
G}\sum_{\begin{subarray}{c}\nu\in G_{0}\\\ (p+\nu)\in G\\\ (q-\nu)\in
G\end{subarray}}\frac{1}{\left\|k_{\nu}\right\|^{2}}\lvert
p+\nu\rangle\\!\langle p\rvert_{i}\lvert q-\nu\rangle\\!\langle
q\rvert_{\ell}$ (12) $\displaystyle V$
$\displaystyle=\frac{2\pi}{\Omega}\sum_{i\neq j=1}^{\eta}\sum_{p,q\in
G}\sum_{\begin{subarray}{c}\nu\in G_{0}\\\ (p+\nu)\in G\\\ (q-\nu)\in
G\end{subarray}}\frac{1}{\left\|k_{\nu}\right\|^{2}}\lvert
p+\nu\rangle\\!\langle p\rvert_{i}\lvert q-\nu\rangle\\!\langle q\rvert_{j}$
(13) $\displaystyle V_{\rm nuc}$
$\displaystyle=\frac{2\pi}{\Omega}\sum_{\ell\neq\kappa=\eta+1}^{L+\eta}\sum_{p,q\in
G}\sum_{\begin{subarray}{c}\nu\in G_{0}\\\ (p+\nu)\in G\\\ (q-\nu)\in
G\end{subarray}}\frac{\zeta_{\ell}\zeta_{\kappa}}{\left\|k_{\nu}\right\|^{2}}\lvert
p+\nu\rangle\\!\langle p\rvert_{\ell}\lvert q-\nu\rangle\\!\langle
q\rvert_{\kappa}\,,$ (14)

where
$G_{0}=[-N^{1/3},N^{1/3}]^{3}\subset\mathbb{Z}^{3}\backslash\\{(0,0,0)\\}$ is
the set of the difference of frequencies from $G$ excluding the singular zero
mode, and $\lvert p\rangle\\!\langle q\rvert_{j}$ is a shorthand for
$I_{1}\otimes\cdots\otimes\lvert p\rangle\\!\langle
q\rvert_{j}\otimes\cdots\otimes I_{\eta}$ where $I_{i}$ is the identity
operator on the $\log N$ bit register associated with electron $i$. We will
store our wavefunction by having a computational basis that encodes
configurations of the electrons in $N$ basis functions such that a
configuration is specified as $\lvert\phi_{1}\phi_{2}\cdots\phi_{\eta}\rangle$
where each $\phi$ encodes the index of an occupied basis function.

The representation of the Hamiltonian we show discretized above is in “first
quantization”. The distinction between “first” and “second” quantization is
fundamentally about how one deals with the symmetries of electrons and nuclei.
Electrons (and some nuclei) are fermions. Some nuclei are bosons. Whereas
identical fermions must be antisymmetric with respect to particle label
exchange, identical bosons must be symmetric with respect to particle label
exchange. For example, if our state corresponds to a valid first quantized
wavefunction for identical particles then the system must have the property
that

$\lvert\psi\rangle=\sum_{\phi_{p}\in\\{\phi\\}}a_{\phi_{1}\cdots\phi_{\eta}}\lvert\phi_{1}\cdots\phi_{i}\cdots\phi_{j}\cdots\phi_{\eta}\rangle=\left(-1\right)^{\pi}\sum_{\phi_{p}\in\\{\phi\\}}a_{\phi_{1}\cdots\phi_{\eta}}\lvert\phi_{1}\cdots\phi_{j}\cdots\phi_{i}\cdots\phi_{\eta}\rangle$
(15)

where $\pi=1$ if fermions and $\pi=0$ if bosons and
$a_{\phi_{1}\cdots\phi_{\eta}}$ are the amplitudes of the associated
computational basis state. No special symmetries need to exist between
distinguishable particles. The work of Berry  _et al._ (2018) shows that one
can symmetrize (or antisymmetrize) any set of $\eta$ registers which each
index $N$ basis functions with a quantum algorithm that has gate complexity
${\cal O}(\eta\log\eta\log N)$. Once symmetrized (or antisymmetrized) any
simulation under the Hamiltonian will maintain that property as a consequence
of (1) commuting with the particle permutation operator between any identical
particles. Since gate complexity of ${\cal O}(\eta\log\eta\log N)$ is
significantly less than the cost of time-evolving these systems, it can be
regarded as a negligible additive cost to the overall simulation.

Finally, one might notice that the spin degree of freedom does not appear in
the Hamiltonians written above, nor does it appear in (1). A spin is
associated with each of the $\eta$ registers but because there are not
magnetic field terms in (1) that would explicitly interact with these spins,
it is not necessary to include the spin label and so we avoid doing so for
simplicity. However, because particles of different spins are distinguishable,
this should be accounted for during the antisymmetrization of the initial
wavefunction.

###  I.2 Overview of the Hamiltonian simulation frameworks deployed in this
work

In the prior section we described the full non-Born-Oppenheimer Hamiltonian in
a plane wave basis and suggested that plane waves would be an ideal
representation in which to perform non-adiabatic dynamics. However, the
algorithms we analyze in detail in this paper focus on the Born-Oppenheimer
Hamiltonian. Still, we emphasize that the non-Born-Oppenheimer application
helps to motivate our work since the algorithms used to simulate the Born-
Oppenheimer Hamiltonian apply to the non-Born-Oppenheimer Hamiltonian with
very minor modifications; the only additional terms in the non-Born-
Oppenheimer Hamiltonian are $U_{\rm non-BO}$ (identical to $V$ up to different
charges) and $T_{\rm nuc}$ (identical to $T$ up to different masses).

Because we will be simulating the Born-Oppenheimer Hamiltonian, which is
typically of interest in electronic structure, our focus will be preparing
eigenstates rather than effecting dynamics. Specifically, we will use the
phase estimation algorithm to sample in the eigenbasis of the Hamiltonian
Kitaev (1995); Abrams and Lloyd (1999). Among several possible uses of this
technique is the refinement of molecular ground states from initial guess
states with non-vanishing support on the exact ground state Aspuru-Guzik  _et
al._ (2005). The requirement for such algorithms is that we are able to
synthesize a unitary circuit encoding eigenvalues of the Born-Oppenheimer
Hamiltonian as a simple function of their eigenphases. In this paper, we will
study two different algorithms for this task based on the qubitization Low and
Chuang (2019) and interaction picture simulation Low and Wiebe (2018)
algorithms of Low _et al._ , with subroutines adapted for simulating first
quantized chemistry by Babbush _et al._ Babbush  _et al._ (2019). We will now
outline these approaches at a high level and establish conventions useful for
later (more technical) sections.

Both of the algorithms we will explore are based (at least in part) on block
encoding the Hamiltonian using qubitization operators Low and Chuang (2019).
This method begins from the observation that we express our target Hamiltonian
as a linear combination of unitaries (often referred to here as “LCU”) Childs
and Wiebe (2012), i.e.

$H=\sum_{\ell}\alpha_{\ell}H_{\ell},$ (16)

where $H_{\ell}$ are unitary operators and $\alpha_{\ell}$ are positive
coefficients (by convention we choose to incorporate any phases into the
$H_{\ell}$). Using this notation we then define the quantum operators
$\mathrm{PREP}_{H}$ and $\mathrm{SEL}_{H}$ as

$\mathrm{PREP}_{H}\lvert
0\rangle=\sum_{\ell}\sqrt{\frac{\alpha_{\ell}}{\lambda}}\lvert\ell\rangle,\qquad\mathrm{SEL}_{H}\lvert\ell\rangle=\lvert\ell\rangle\otimes
H_{\ell},\qquad\lambda=\sum_{\ell}\alpha_{\ell}\,$ (17)

and we notice that $\langle
0\rvert\mathrm{PREP}_{H}^{\dagger}\cdot\mathrm{SEL}_{H}\cdot\mathrm{PREP}_{H}\lvert
0\rangle=H/\lambda$ Low and Wiebe (2018). Because of this property, if we can
realize $\mathrm{PREP}_{H}$ and $\mathrm{SEL}_{H}$ as quantum circuits then by
calling them in a certain sequence it is possible to “block encode” the
Hamiltonian in a particular subspace of a qubitization operator that we define
as

$Q=\left(2\lvert 0\rangle\\!\langle
0\rvert-I\right)\cdot\mathrm{PREP}_{H}^{\dagger}\cdot\mathrm{SEL}_{H}\cdot\mathrm{PREP}_{H}\,.$
(18)

What this means is that $Q$ acts as a quantum walk Szegedy (2004) with
eigenvalues $e^{\pm i\arccos(E_{k}/\lambda)}$ Low and Chuang (2019), related
to the energy $E_{k}$ of the target Hamiltonian $H$. The work of Berry  _et
al._ (2018); Poulin  _et al._ (2017) showed that we can use these quantum
walks to project to an eigenstate of the Hamiltonian and estimate its
eigenvalue by performing phase estimation on $Q$ and applying the cosine
function on the resulting phase. To estimate with accuracy $\epsilon$, we need
$\mathcal{O}(\lambda/\epsilon)$ uses of controlled-$Q$ operations.
Alternatively, the resource requirements of performing phase estimation can be
reduced by using the multiplexed qubitization operation Babbush  _et al._
(2018a):

$\lvert 0\rangle\\!\langle 0\rvert\otimes Q^{\dagger}+\lvert
1\rangle\\!\langle 1\rvert\otimes Q\,,$ (19)

instead of the controlled-$Q$ operation. This amplifies the difference between
the eigenvalues by a factor of 2, reducing the complexity by a factor of 2.
Provided $\mathrm{SEL}_{H}$ is self-inverse, this control can be achieved
simply by making the reflection controlled by a qubit from the control
register for phase estimation, as per Figure 2 of Lee  _et al._ (2021). This
control only adds a cost of one Toffoli per step.

We note that in addition to being useful for phase estimation, one can also
realize time-evolution under the Born-Oppenheimer Hamiltonian for duration $t$
by invoking these quantum walks $Q$ within the quantum signal processing
framework of Low _et al._ Low and Chuang (2017) a number of times scaling as

${\cal O}\left(\lambda
t+\frac{\log(1/\epsilon)}{\log\log(1/\epsilon)}\right)\,.$ (20)

Thus, while our explicit focus will be on phase estimation of the Born-
Oppenheimer Hamiltonian, the same algorithms can be used to perform time-
evolution under the non-Born-Oppenheimer Hamiltonian, with trivial
modifications and marginal additional cost. The similarity of the algorithms
enabling these disparate applications results from our choice to work in a
first quantized representation with basis functions that are independent of
nuclear coordinates; if we were instead to work in second quantization or use
Gaussian orbitals, there would be a substantial gap between the algorithms
required for these two applications. Getting back to the implementation of the
qubitization algorithm, in order to express the Hamiltonian of (8) as a linear
combination of unitaries we rewrite the $T$, $U$ and $V$ operators as follows:

$\displaystyle T$
$\displaystyle=\frac{2\pi^{2}}{\Omega^{2/3}}\sum_{j=1}^{\eta}\sum_{p\in
G}\sum_{w\in\\{x,y,z\\}}\sum_{r=0}^{n_{p}-2}\sum_{s=0}^{n_{p}-2}2^{r+s}p_{w,r}p_{w,s}\lvert
p\rangle\\!\langle p\rvert_{j},$ (21) $\displaystyle U$
$\displaystyle=\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{2\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\sum_{j=1}^{\eta}\sum_{b\in\\{0,1\\}}\left(-e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}\sum_{q\in G}(-1)^{b[(q-\nu)\notin G]}\lvert q-\nu\rangle\\!\langle
q\rvert_{j}\right),$ $\displaystyle V$ $\displaystyle=\sum_{\nu\in
G_{0}}\frac{\pi}{\Omega\left\lVert k_{\nu}\right\rVert^{2}}\sum_{i\neq
j=1}^{\eta}\sum_{b\in\\{0,1\\}}\left(\sum_{p,q\in G}(-1)^{b([(p+\nu)\notin
G]\lor[(q-\nu)\notin G])}\lvert p+\nu\rangle\langle p\rvert_{i}\cdot\lvert
q-\nu\rangle\\!\langle q\rvert_{j}\right)\,,$

where

$n_{p}=\left\lceil\log\left(N^{1/3}+1\right)\right\rceil\,$ (22)

represents the number of qubits required to store a signed binary
representation of one component of the momentum of a single electron. We can
see that the potential operators $U$ and $V$ are now expressed as linear
combinations of unitaries. To express the $T$ operator as a linear combination
of unitaries we note that the quantity $p_{w,r}$ is bit $r$ of component $w$
of $p$, and the sum of these products of bits over $r$ and $s$ gives
$\|p\|^{2}$; we apply the identity

$p_{w,r}p_{w,s}=\frac{1-(-1)^{p_{w,r}p_{w,s}}}{2},$ (23)

which gives

$T=\frac{\pi^{2}}{\Omega^{2/3}}\sum_{j=1}^{\eta}\sum_{w\in\\{x,y,z\\}}\sum_{r=0}^{n_{p}-2}\sum_{s=0}^{n_{p}-2}2^{r+s}\sum_{b\in\\{0,1\\}}\left(\sum_{p\in
G}(-1)^{b(p_{w,r}p_{w,s}\oplus 1)}\lvert p\rangle\\!\langle
p\rvert_{j}\right).$ (24)

With the Hamiltonian components expressed as a linear combination of unitaries
we can now discuss the normalizing parameter $\lambda$ appearing in (17),
which will have significant ramifications for the complexity of both
algorithms studied in this work. The parameter $\lambda$ is an induced one-
norm of the Hamiltonian expressed as a linear combination of unitaries. For
the entire Hamiltonian we have $\lambda=\lambda_{T}+\lambda_{U}+\lambda_{V}$
where

$\displaystyle\lambda_{T}$
$\displaystyle=\frac{6\eta\pi^{2}}{\Omega^{2/3}}\left(2^{n_{p}-1}\\!-1\right)^{2}={\cal
O}\left(\frac{\eta
N^{2/3}}{\Omega^{2/3}}\right),\qquad\qquad\lambda_{U}=\frac{\eta\sum_{\ell}\zeta_{\ell}}{\pi\Omega^{1/3}}\lambda_{\nu}={\cal
O}\left(\frac{\eta^{2}N^{1/3}}{\Omega^{1/3}}\right)$ (25)
$\displaystyle\lambda_{V}$
$\displaystyle=\frac{\eta(\eta-1)}{2\pi\Omega^{1/3}}\lambda_{\nu}={\cal
O}\left(\frac{\eta^{2}N^{1/3}}{\Omega^{1/3}}\right),\qquad\lambda_{\nu}=\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert^{2}}\leq\int_{0}^{2\pi}\\!\\!\\!\\!\int_{0}^{\pi}\\!\\!\\!\int_{0}^{N^{1/3}}\\!\\!\\!\\!\\!\\!{\rm
d}r\,{\rm d}\phi\,{\rm
d}\theta\,\left(\frac{1}{r^{2}}\right)r^{2}\sin\theta={\cal
O}\left(N^{1/3}\right)\,,$

where in the asymptotic scaling for $\lambda_{U}$ we have used the fact that
when including all electrons explicitly in a charge neutral simulation,
$\sum_{\ell}\zeta_{\ell}=\eta$. Further, we note that it is often appropriate
to take $\Omega\propto\eta$ (e.g., whenever growing condensed phase systems
towards their thermodynamic limit), which simplifies the asymptotic
complexities to $\lambda_{T}={\cal O}(\eta^{1/3}N^{2/3})$ and
$\lambda_{U}\propto\lambda_{V}={\cal O}(\eta^{5/3}N^{1/3})$.

The main challenge of realizing an algorithm like qubitization with low
complexity is coming up with efficient circuits for the operators defined in
(17). Describing such algorithms asymptotically for first quantized chemistry
is one of the main contributions of the work of Babbush _et al._ Babbush  _et
al._ (2019) and this work will loosely follow the methods outlined therein.
However, the fact that this paper is an order of magnitude longer than the
former is a testament to the fact that bounding asymptotic complexity is
significantly less work than coming up with, improving, and compiling, a
practical implementation of that algorithm with low constant factors in a
fault-tolerant cost model. Thus, the main technical contribution of Section II
is to describe a detailed algorithm for realizing the qubitization operators
corresponding to the first quantized plane wave electronic structure
Hamiltonian. We introduce various new techniques, including a strategy for
controlled swapping the momentum registers (Section II.2 and Section II.4) and
a proposal for performing the kinetic term when the state preparation for the
potential fails (Appendix D), which lead to a significant reduction in the
algorithmic complexity. Ultimately, we show that the operator $Q$ can be
realized with a gate complexity of $\widetilde{\cal O}(\eta)$ and that one can
implement phase estimation to estimate eigenvalues of (8) to within precision
$\epsilon$ with a leading order Toffoli complexity given by Theorem 4.

By combining the $\widetilde{\cal O}(\eta)$ gate complexity of $Q$ with the
${\cal O}(\lambda/\epsilon)$ times that $Q$ must be invoked within phase
estimation, we see that the asymptotic complexity of the qubitized algorithm
for phase estimation is

$\widetilde{\cal
O}\left(\frac{\left(\lambda_{T}+\lambda_{U}+\lambda_{V}\right)\eta}{\epsilon}\right)=\widetilde{\cal
O}\left(\frac{\eta^{2}N^{2/3}}{\epsilon\,\Omega^{2/3}}+\frac{\eta^{3}N^{1/3}}{\epsilon\,\Omega^{1/3}}\right)=\widetilde{\cal
O}\left(\frac{\eta^{4/3}N^{2/3}+\eta^{8/3}N^{1/3}}{\epsilon}\right)\,$ (26)

as reported for the qubitized algorithm complexity in Babbush  _et al._
(2019). While the scaling of this algorithm in terms of $\eta$ is worse than
the scaling in terms of $N$, for practical applications we may be more
concerned about the scaling with $N$. The reason is because in order to
accurately solve the electronic structure problem in a plane wave basis it
will be necessary to use a very large number of plane waves. For example, some
enticing applications may require roughly $\eta=100$ and $N=10^{6}$. While
such simulations would be impossible within second quantization because that
would (in this case) require at least $N=10^{6}$ logical qubits, within first
quantization the space complexity required for the system register is only
$\eta\log N\approx 2\times 10^{3}$ qubits. Here and throughout the manuscript
we use the convention that $\log$ indicates a logarithm base 2.

The vast majority of this space complexity arises from representing the
simulated system, which has basis states of the form

$\lvert p_{1}\rangle\cdots\lvert p_{j}\rangle\cdots\lvert p_{\eta}\rangle.$
(27)

Representing coordinates of the orbitals explicitly, we have

$\lvert p_{1x},p_{1y},p_{1z}\rangle\cdots\lvert
p_{jx},p_{jy},p_{jz}\rangle\cdots\lvert p_{\eta x},p_{\eta y},p_{\eta
z}\rangle.$ (28)

As in (6), these registers hold values in the range
$[-(N^{1/3}-1)/2,(N^{1/3}-1)/2]$. Thus, we will use $n_{p}$ qubits with
$N^{1/3}=2^{n_{p}}-1$ to encode each register. Furthermore, for reasons of
algorithmic efficiency described later, we will use signed binary numbers to
represent each register of length $n_{p}$, so

$\lvert q\rangle=\lvert q_{\text{sign}},q_{n_{p}-2},\ldots,q_{0}\rangle.$ (29)

The total number of qubits in the simulated system is thus

$n_{s}=3\,\eta\,n_{p}.$ (30)

So while the space complexity is likely not a problem for performing
simulations at large $N$, the gate complexity scaling as $\widetilde{\cal
O}(N^{2/3})$ is somewhat concerning and we might wonder if it is possible to
reduce this scaling. This is the motivation for the second algorithm we will
explore (also introduced in Babbush  _et al._ (2019)), which provides an
algorithm that scales linearly in $\lambda_{U}+\lambda_{V}$ but avoids scaling
polynomially in $\lambda_{T}$ (which is where the $N^{2/3}$ scaling was
entering previously). The approach is based on the interaction picture
simulation technique introduced by Low _et al._ Low and Wiebe (2018) which
begins by supposing that we want to simulate a target Hamiltonian containing
two terms $A$ and $B$, where $\left\lVert A\right\rVert$ is significantly
larger than $\left\lVert B\right\rVert$ but time-evolution of $A$ can be fast-
forwarded. Then, we can represent the evolution in the interaction picture of
the rotating frame of $A$,

$e^{-\mathrm{i}\tau(A+B)}=e^{-\mathrm{i}\tau
A}\mathcal{T}\exp\left(-\mathrm{i}\int_{0}^{\tau}ds\
e^{\mathrm{i}sA}Be^{-\mathrm{i}sA}\right),$ (31)

where $\mathcal{T}$ is the time ordering operator. The potential benefit of
using this representation is that the larger term $A$ only appears in the
exponent of matrix exponentials and its implementation can be efficient even
when $\left\lVert A\right\rVert$ is large. For the first quantized plane wave
Hamiltonian in the regime that $N\gg\eta$, we have the norm relation
$\left\lVert T\right\rVert\geq\left\lVert U+V\right\rVert$ so we choose $A=T$
and $B=U+V$ Babbush  _et al._ (2019).

To perform quantum simulation in the interaction picture, we truncate and
discretize the Dyson series of the time-ordered exponential, obtaining

$\displaystyle e^{-\mathrm{i}(A+B)\tau}$ $\displaystyle=e^{-\mathrm{i}\tau
A}\mathcal{T}\exp\left(-\mathrm{i}\int_{0}^{\tau}ds\,e^{\mathrm{i}sA}Be^{-\mathrm{i}sA}\right)$
(32)
$\displaystyle=\sum_{k=0}^{\infty}\frac{(-\mathrm{i})^{k}}{k!}\int_{0}^{\tau}d\tau_{1}\int_{0}^{\tau}d\tau_{2}\cdots\int_{0}^{\tau}d\tau_{k}\,e^{-\mathrm{i}(\tau-\tau^{\prime}_{k})A}Be^{-\mathrm{i}(\tau^{\prime}_{k}-\tau^{\prime}_{k-1})A}B\ldots
Be^{-\mathrm{i}(\tau^{\prime}_{2}-\tau^{\prime}_{1})A}Be^{-\mathrm{i}\tau_{1}A}$
$\displaystyle\approx\sum_{k=0}^{K}\frac{(-\mathrm{i}\tau)^{k}}{M^{k}k!}\sum_{m_{1}=0}^{M-1}\sum_{m_{2}=0}^{M-1}\cdots\sum_{m_{k}=0}^{M-1}e^{-\mathrm{i}\tau(M-1/2-m^{\prime}_{k})A/M}Be^{-\mathrm{i}\tau(m^{\prime}_{k}-m^{\prime}_{k-1})A/M}B\ldots$
$\displaystyle\quad\ldots
Be^{-\mathrm{i}\tau(m^{\prime}_{2}-m^{\prime}_{1})A/M}Be^{-\mathrm{i}\tau(m^{\prime}_{1}+1/2)A/M}$

where $0\leq\tau^{\prime}_{1}\leq\ldots\leq\tau^{\prime}_{k}\leq\tau$ are
sorted times from $\tau_{1},\ldots,\tau_{k}$, and $0\leq
m^{\prime}_{1}\leq\ldots\leq m^{\prime}_{k}\leq M-1$ are sorted integers from
$m_{1},\ldots,m_{k}$. The resulting expression is only an approximation of the
ideal evolution, but the truncation and discretization error can be made
arbitrarily small by choosing $K$ and $M$ sufficiently large. These parameters
give a contribution to the complexity which scales as a polynomial in
$\log(1/\epsilon)$, which is omitted when giving complexity using
$\widetilde{O}$.

As discussed in Section III, the bottleneck requirements for realizing this
Dyson series are that one implements $e^{-\mathrm{i}T}$ and also block encodes
$U+V$ a number of times scaling as $(\lambda_{U}+\lambda_{V})/\epsilon$. For
the most part, one can reuse parts of the circuits in Section II to perform
the block encoding of $U+V$ with complexity $\widetilde{\cal O}(\eta)$.
Furthermore, it is relatively straightforward to provide circuits that fast-
forward the operator $e^{-\mathrm{i}T}$ with gate complexity $\widetilde{\cal
O}(\eta)$. Thus, this suggests that the overall asymptotic complexity will
scale as

$\widetilde{\cal
O}\left(\frac{\left(\lambda_{U}+\lambda_{V}\right)\eta}{\epsilon}\right)=\widetilde{\cal
O}\left(\frac{\eta^{3}N^{1/3}}{\epsilon\,\Omega^{1/3}}\right)=\widetilde{\cal
O}\left(\frac{\eta^{8/3}N^{1/3}}{\epsilon}\right)\,.$ (33)

However, the most challenging part of the interaction picture algorithm is to
efficiently tie together these primitives within the Dyson series simulation
framework and propagate errors from that simulation through the phase
estimation circuits. While past work has discussed quantum simulation of the
Dyson series Kieferová  _et al._ (2019); Low and Wiebe (2018) in asymptotic
terms, no prior work has provided concrete circuits for the task and worked
out the constant factors in the scaling. One of the primary contributions of
Section III is to work out those circuits in the context of chemistry,
optimize their implementations, and bound the overall Toffoli cost, as
reported in Theorem 5. In doing so, we make several technical contributions,
including a novel qubitization of the pre-amplified Dyson series (Section
III.1), an efficient strategy to update the kinetic energy (Section III.2),
and a higher-order discretization of the time integrals (Appendix H), which
greatly improves over the naive approaches. We note that we also work out the
constant factor scalings for the interaction picture simulation in a more
general context (not related to chemistry), which we report in Appendix I.

Finally, in Section IV our paper concludes with numerical comparisons between
these two algorithms assessing the resources required to deploy them to real
molecular systems. In addition to comparing the two algorithms, we also
compare to prior work analyzing the fault-tolerant viability of second
quantized Gaussian orbital approaches. Although the comparison is somewhat
nuanced, the first quantized approaches developed here appear to give lower
Toffoli complexities than any past work for realistic simulations of both
material Hamiltonians and non-periodic molecules. We note that throughout this
work we focus on estimating the constant factor Toffoli complexity and number
of ancilla required. We focus on Toffoli complexity (as opposed to the
complexity of other gates) because these are the operations that would be most
challenging to implement within most practical error-correcting codes such as
the surface code Fowler  _et al._ (2012). Furthermore, the algorithms
described here have the same asymptotic Clifford complexities as the Toffoli
complexities that we report. This makes it more likely that the implementation
would actually be bottlenecked by Toffoli distillation. This stands in
contrast to other recent chemistry algorithms compiled to fault-tolerance with
extensive use of quantum read-only memory (QROM) Babbush  _et al._ (2018a);
Low  _et al._ (2018). Approaches relying heavily on QROM such as Babbush  _et
al._ (2018a); Berry  _et al._ (2019); von Burg  _et al._ (2020); Lee  _et al._
(2021) tend to have a worse asymptotic Clifford complexity than Toffoli
complexity, thus necessitating a careful layout of how the algorithm might be
realized using lattice surgery or other surface code fault-tolerant gates, in
order to assess whether state distillation is actually the dominant
bottleneck.

There are a number of new techniques which we introduce which can be used in
many other applications. One is to note that in some cases amplitude
amplification can be eliminated from the state preparation procedure for the
linear combination of unitaries. The reason for this is that the Hamiltonian
can be composed of a number of parts, where that state preparation is only
needed for one part of the Hamiltonian. When the state preparation is regarded
as having “failed”, we can simply flag that as a part of the state which will
control the select for another part of the Hamiltonian. It is useful to
eliminate the amplitude amplification, which would otherwise triple the cost
(or more if further amplification steps are needed).

Another new technique is to eliminate the arithmetic for some parts of the
Hamiltonian. If a component of the Hamiltonian gives a multiplication of the
state by a function of the basis state, then it can be broken up into a sum of
individual bit products. The bit product only requires a single Toffoli, and
the sum is implemented via the linear combination of unitaries. In our case
the component of the Hamiltonian is the kinetic energy, and we break up the
kinetic energy into a sum of products of bits of the momenta. This
dramatically reduces the Toffoli complexity.

In the case of the interaction picture we can greatly reduce the complexity
needed for dealing with the kinetic energy (though not completely eliminate
arithmetic like we can for qubitisation). The idea is to use a kinetic energy
register. This register can be updated with relatively small complexity and
used for the phasing. Another approach we use for the interaction picture is
to use a novel qubitisation of the Dyson series. This again eliminates the
need for an amplitude amplification which would otherwise triple the
complexity. This method works generally for applications where the goal is to
measure the eigenvalue, rather than to generate the time evolution under the
Hamiltonian.

Our implementation of these algorithms allows us to report the first constant
factors for the scaling of any first quantized chemistry algorithm (or any
Dyson series based simulation), and also optimizes those constant factors more
than a thousandfold. This allows us to compare the viability of these
approaches to other methods based on second quantization, and to each other.
Ultimately, we find that the qubitization algorithm is usually more practical
and that first quantized methods often require much less surface code
spacetime volume for simulating millions of plane waves than the best second
quantized algorithms require for simulating hundreds of Gaussian orbitals.

##  II The qubitization based algorithm

We now give an overview of the qubitization algorithm and its circuit
implementation for simulating the first-quantized quantum chemistry. For this
purpose, we seek to represent the target Hamiltonian as a linear combination
of unitaries $H=\sum_{\ell=1}^{L}\alpha_{\ell}H_{\ell}$, where $H_{\ell}$ are
unitary operators and $\alpha_{\ell}$ are positive coefficients. Recall from
Section I.2 that we can construct a state preparation subroutine
$\mathrm{PREP}_{H}$ and a selection subroutine $\mathrm{SEL}_{H}$, such that
the block encoding $\langle
0\rvert\mathrm{PREP}_{H}^{\dagger}\cdot\mathrm{SEL}_{H}\cdot\mathrm{PREP}_{H}\lvert
0\rangle=H/\lambda$ holds with normalization factor
$\lambda=\sum_{\ell}\alpha_{\ell}$. Then, the operator $Q=\left(2\lvert
0\rangle\langle
0\rvert-I\right)\cdot\mathrm{PREP}_{H}^{\dagger}\cdot\mathrm{SEL}_{H}\cdot\mathrm{PREP}_{H}$
has eigenvalues $e^{\pm i\arccos(E_{k}/\lambda)}$ corresponding to the energy
$E_{k}$ of $H$. We can thus estimate the spectrum of $H$ by performing phase
estimation on $Q$ with classical postprocessing.

For the chemistry Hamiltonian $H=T+U+V$, we see from Section I.2 that the
terms $T$, $U$, and $V$ can indeed be expressed as linear combinations of
unitaries. This means that we have $\mathrm{PREP}_{T}$, $\mathrm{PREP}_{U+V}$
and $\mathrm{SEL}_{T}$, $\mathrm{SEL}_{U+V}$, such that the block encodings

$\langle
0\rvert\mathrm{PREP}_{T}^{\dagger}\cdot\mathrm{SEL}_{T}\cdot\mathrm{PREP}_{T}\lvert
0\rangle=\frac{T}{\lambda_{T}},\qquad\langle
0\rvert\mathrm{PREP}_{U+V}^{\dagger}\cdot\mathrm{SEL}_{U+V}\cdot\mathrm{PREP}_{U+V}\lvert
0\rangle=\frac{U+V}{\lambda_{U}+\lambda_{V}}$ (34)

hold with normalization factors $\lambda_{T}$ and $\lambda_{U}+\lambda_{V}$,
respectively. As a result, we find that the state preparation subroutine

$\left(\sqrt{\frac{\lambda_{T}}{\lambda_{T}+\lambda_{U}+\lambda_{V}}}\lvert
0\rangle+\sqrt{\frac{\lambda_{U}+\lambda_{V}}{\lambda_{T}+\lambda_{U}+\lambda_{V}}}\lvert
1\rangle\right)\otimes\mathrm{PREP}_{T}\lvert
0\rangle\otimes\mathrm{PREP}_{U+V}\lvert 0\rangle$ (35)

and the selection subroutine

$\lvert 0\rangle\\!\langle 0\rvert\otimes\mathrm{SEL}_{T}\otimes I+\lvert
1\rangle\\!\langle 1\rvert\otimes I\otimes\mathrm{SEL}_{U+V}$ (36)

block-encode the target Hamiltonian $H$ with normalization factor
$\lambda_{T}+\lambda_{U}+\lambda_{V}$, giving a circuit implementation of the
qubitization algorithm.

There is a subtle issue with the above naive analysis. While the kinetic
operator $T$ can be encoded with a normalization factor close to
$\lambda_{T}$, in the most efficient approach for preparation the potential
operators $U+V$ are only encoded with a normalization factor close to
$4(\lambda_{U}+\lambda_{V})$. To be more specific, we can construct a state
preparation

$\widetilde{\mathrm{PREP}}_{U+V}\lvert 0\rangle\otimes
I\approx\frac{1}{2}\lvert
0\rangle\otimes\mathrm{PREP}_{U+V}+\frac{\sqrt{3}}{2}\lvert
1\rangle\otimes\mathrm{PREP}_{U+V}^{\bot}$ (37)

that only succeeds with a probability close to $1/4$, where $\lvert 0\rangle$
flags that the preparation is successful and $\lvert 1\rangle$ indicates a
failure. If we were to use such a probabilistic subroutine in our block
encoding, we would implement qubitization with a normalization factor close to
$\lambda_{T}+4\left(\lambda_{U}+\lambda_{V}\right)$. Alternatively, we could
perform a single step of oblivious amplitude amplification to boost the
success probability close to unity, but this triples the cost of the state
preparation and can introduce overhead to the overall gate complexity.

Our new approach improves the above implementations by proceeding even when
the state preparation of $U+V$ fails. To elaborate, first consider the special
case where $\lambda_{T}=3(\lambda_{U}+\lambda_{V})$. Then, we propose to
perform the selection subroutine of $T$ following the failure of the
preparation of $U+V$, i.e.,

$I\otimes\lvert 0\rangle\\!\langle
0\rvert\otimes\mathrm{SEL}_{U+V}+\mathrm{SEL}_{T}\otimes\lvert
1\rangle\\!\langle 1\rvert\otimes I.$ (38)

This new selection operator, combined with the state preparation

$\mathrm{PREP}_{T}\lvert 0\rangle\otimes\widetilde{\mathrm{PREP}}_{U+V}\lvert
0,0\rangle,$ (39)

gives the block encoding

$\frac{3T}{4\lambda_{T}}+\frac{U+V}{4(\lambda_{U}+\lambda_{V})}=\frac{T}{4(\lambda_{U}+\lambda_{V})}+\frac{U+V}{4(\lambda_{U}+\lambda_{V})}=\frac{T+U+V}{\lambda_{T}+\lambda_{U}+\lambda_{V}}$
(40)

with the desired normalization factor $\lambda_{T}+\lambda_{U}+\lambda_{V}$.

The general case can be handled by using an ancilla qubit to adjust the
relative amplitudes between $T$ and $U+V$. If
$\lambda_{T}<3(\lambda_{U}+\lambda_{V})$, then we introduce an ancilla qubit
and only perform the selection of $T$ conditioned on both the failure of the
preparation of $U+V$ AND the ancilla qubit is in state $\lvert 0\rangle$,
i.e.,

$I\otimes I\otimes\lvert 0\rangle\\!\langle
0\rvert\otimes\mathrm{SEL}_{U+V}+\lvert 0\rangle\\!\langle
0\rvert\otimes\mathrm{SEL}_{T}\otimes\lvert 1\rangle\\!\langle 1\rvert\otimes
I+\lvert 1\rangle\\!\langle 1\rvert\otimes I\otimes\lvert 1\rangle\\!\langle
1\rvert\otimes I.$ (41)

This selection operator, combined with the state preparation

$\left(\cos\theta\lvert 0\rangle+\sin\theta\lvert
1\rangle\right)\otimes\mathrm{PREP}_{T}\lvert
0\rangle\otimes\widetilde{\mathrm{PREP}}_{U+V}\lvert 0,0\rangle,$ (42)

gives the block encoding

$\frac{U+V}{4(\lambda_{U}+\lambda_{V})}+\frac{3T\cos^{2}\theta}{4\lambda_{T}}$
(43)

up to a shifting by multiples of the identity operator. For this to be an
accurate encoding of the Hamiltonian, this expression must be proportional to
$U+V+T$. That means we need
$1/[4(\lambda_{U}+\lambda_{V})]=(3\cos^{2}\theta)/(4\lambda_{T})$; solving for
$\theta$, we find that
$\theta=\arccos\sqrt{\lambda_{T}/3(\lambda_{U}+\lambda_{V})}$, which gives the
target Hamiltonian encoded with normalization factor
$4(\lambda_{U}+\lambda_{V})$.

In the case where $\lambda_{T}>3(\lambda_{U}+\lambda_{V})$, we would still
introduce an ancilla qubit to adjust the relative amplitudes. However, we now
perform the selection of $T$ conditioned on failure of the preparation
subroutine OR the ancilla qubit being in state $\lvert 0\rangle$ (and perform
$\mathrm{SEL}_{T}$ in the remaining situations), i.e.,

$\lvert 1\rangle\\!\langle 1\rvert\otimes I\otimes\lvert 0\rangle\\!\langle
0\rvert\otimes\mathrm{SEL}_{U+V}+\lvert 1\rangle\\!\langle
1\rvert\otimes\mathrm{SEL}_{T}\otimes\lvert 1\rangle\\!\langle 1\rvert\otimes
I+\lvert 0\rangle\\!\langle 0\rvert\otimes\mathrm{SEL}_{T}\otimes I\otimes I.$
(44)

Combined with the state preparation

$\left(\cos\theta\lvert 0\rangle+\sin\theta\lvert
1\rangle\right)\otimes\mathrm{PREP}_{T}\lvert
0\rangle\otimes\widetilde{\mathrm{PREP}}_{U+V}\lvert 0,0\rangle,$ (45)

this gives the block encoding

$\frac{(U+V)\sin^{2}\theta}{4(\lambda_{U}+\lambda_{V})}+\frac{3T\sin^{2}\theta}{4\lambda_{T}}+\frac{T\cos^{2}\theta}{\lambda_{T}}.$
(46)

For this to be proportional to $U+V+T$ we need

$\frac{\sin^{2}\theta}{4(\lambda_{U}+\lambda_{V})}=\frac{3\sin^{2}\theta}{4\lambda_{T}}+\frac{\cos^{2}\theta}{\lambda_{T}}.$
(47)

Solving for $\theta$, we find that
$\theta=\arcsin(2\sqrt{(\lambda_{U}+\lambda_{V})/(\lambda_{T}+\lambda_{U}+\lambda_{V})})$,
which gives the target Hamiltonian encoded with the normalization factor
$\lambda_{T}+\lambda_{U}+\lambda_{V}$.

In the following, we sketch the high-level ideas behind the circuit
implementation of the state preparation and selection subroutines, and leave a
discussion of the details to subsequent sections, including a more detailed
description of the implementation in Section II.1. For simplicity, we focus on
the case where $\lambda_{T}<3(\lambda_{U}+\lambda_{V})$ although similar
analysis applies to the other cases as well. Note that $\lambda_{T}$ is not
significantly larger than $\lambda_{U}+\lambda_{V}$ in this parameter regime
and this is when the qubitization algorithm is expected to be advantageous.

For state preparation, we aim to prepare the quantum state

$\displaystyle(\cos\theta\lvert 0\rangle+\sin\theta\lvert
1\rangle)_{a}\lvert+\rangle_{b}\frac{1}{\sqrt{\eta}}\\!\left(\\!\sqrt{\eta-1}\lvert
0\rangle_{c}\\!\\!\sum_{i\neq j=1}^{\eta}\lvert i\rangle_{d}\lvert
j\rangle_{e}+\lvert 1\rangle_{c}\\!\sum_{j=1}^{\eta}\lvert j\rangle_{d}\lvert
j\rangle_{e}\\!\right)\\!\left(\\!\frac{1}{\sqrt{3}}\\!\sum_{w=0}^{2}\lvert
w\rangle_{f}\\!\right)\\!$
$\displaystyle\otimes\left(\\!\frac{1}{2^{n_{p}-1}-1}\\!\sum_{r,s=0}^{n_{p}-2}2^{(r+s)/2}\lvert
r\rangle_{g}\lvert
s\rangle_{h}\\!\right)\left(\sqrt{\frac{\lambda_{U}}{\lambda_{U}+\lambda_{V}}}\lvert
0\rangle+\sqrt{\frac{\lambda_{V}}{\lambda_{U}+\lambda_{V}}}\lvert
1\rangle\right)_{i}$
$\displaystyle\otimes\left(\sqrt{\frac{p_{\nu}}{\lambda_{\nu}}}\lvert
0\rangle_{j}\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert}\lvert\nu\rangle_{k}+\sqrt{1-p_{\nu}}\lvert
1\rangle_{j}\lvert\nu^{\bot}\rangle_{k}\right)\left(\frac{1}{\sqrt{\sum_{\ell}\zeta_{\ell}}}\sum_{\ell=1}^{L}\sqrt{\zeta_{\ell}}\lvert\ell\rangle\right)_{l},$
(48)

with $p_{\nu}$ the probability of successfully preparing the momentum state,
where the first ancilla register is used to adjust the relative amplitudes
between $T$ and $U+V$ (as explained above), and the $i$ register (holding
state with relative amplitudes $\sqrt{\lambda_{U}}$ and $\sqrt{\lambda_{V}}$)
is used to select between $U$ and $V$. Here, the preparation subroutine
$\mathrm{PREP}_{T}$ acts on the registers $b,f,g,h$ and outputs the
corresponding states, whereas $\widetilde{\mathrm{PREP}}_{U+V}$ acts on
registers $c,d,e,j,k,l$. Note that as some of these registers are used for
multiple terms in the Hamiltonian, the assignment of common registers such as
$\lvert\cdot\rangle_{b}$ to any particular preparation circuit is arbitrary.
Quantum circuits for these preparation subroutines are detailed in Section
II.2 and Section II.3, respectively. The normalization factor $\lambda_{\nu}$
is defined in (25), and $\lvert\nu^{\bot}\rangle$ is a state obtained in the
failure of initializing the momentum register. We will see that the
probability of successfully preparing the momentum state
$(1/\sqrt{\lambda_{\nu}})\sum_{\nu\in
G_{0}}(1/\left\lVert\nu\right\rVert)\lvert\nu\rangle$ is close to $1/4$, which
is consistent with our discussion above.

Note that several important details, including the flagging of $i\neq j$ in
register $c$ and the precision in preparing the ancilla state
$\cos\theta\lvert 0\rangle+\sin\theta\lvert 1\rangle$, are omitted in our
analysis for exposition purposes. We provide full details on how the encodings
between $T$, $U$, and $V$ are performed and what precision the ancilla state
is prepared with in Appendix D.

We now describe the key steps in the implementation of the selection
operations. For $\mathrm{SEL}_{T}$, we can compute the product
$p_{w,r}p_{w,s}$ in an ancilla register and use that to perform a controlled
$Z$ gate on $\lvert+\rangle$, followed by a final uncomputation. The net
effect of this transformation is

$\mathrm{SEL}_{T}:\lvert b\rangle_{b}\lvert j\rangle_{e}\lvert
w\rangle_{f}\lvert r\rangle_{g}\lvert s\rangle_{h}\lvert
p_{j}\rangle\mapsto(-1)^{b(p_{w,r}p_{w,s}\oplus 1)}\lvert b\rangle_{b}\lvert
j\rangle_{e}\lvert w\rangle_{f}\lvert r\rangle_{g}\lvert s\rangle_{h}\lvert
p_{j}\rangle,$ (49)

where $\lvert b\rangle$ is either the component $\lvert 0\rangle$ or $\lvert
1\rangle$ of the state $\lvert+\rangle$ and $\lvert p_{j}\rangle$ is the
$j^{\rm th}$ particle’s momentum. To implement $\mathrm{SEL}_{U+V}$, we first
add $\nu$ into the $i$th momentum register. Since this should only be done for
$V$, we perform the addition controlled by the register selecting between $U$
and $V$. Next, we subtract $\nu$ from the $j$th momentum register for both $U$
and $V$, and the subtraction does not need to be controlled by the
corresponding selection register. We now check if $(p-\nu)\notin G$ (for $U$)
or $[(p+\nu)\notin G]\lor[(q-\nu)\notin G]$ (for $V$) and use the result to
perform controlled $Z$ gate on $\lvert+\rangle$. For the case of $U$, we also
perform the phase rotation $-e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$. Overall,
this implements

$\displaystyle\mathrm{SEL}_{U}:\lvert b\rangle_{b}\lvert j\rangle_{e}\lvert
0\rangle_{i}\lvert\nu\rangle_{k}\lvert\ell\rangle_{l}\lvert
q_{j}\rangle\mapsto$ $\displaystyle\ \lvert b\rangle_{b}\lvert
j\rangle_{e}\lvert 0\rangle_{i}\lvert\nu\rangle_{k}\lvert\ell\rangle_{l}\lvert
q_{j}-\nu\rangle$ $\displaystyle\mapsto$ $\displaystyle\ (-1)^{b[(p-\nu)\notin
G]}\lvert b\rangle_{b}\lvert j\rangle_{e}\lvert
0\rangle_{i}\lvert\nu\rangle_{k}\lvert\ell\rangle_{l}\lvert q_{j}-\nu\rangle$
$\displaystyle\mapsto$ $\displaystyle\ -e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}(-1)^{b[(p-\nu)\notin G]}\lvert b\rangle_{b}\lvert j\rangle_{e}\lvert
0\rangle_{i}\lvert\nu\rangle_{k}\lvert\ell\rangle_{l}\lvert q_{j}-\nu\rangle$
(50)

and

$\displaystyle\mathrm{SEL}_{V}:\lvert b\rangle_{b}\lvert i\rangle_{d}\lvert
j\rangle_{e}\lvert 1\rangle_{i}\lvert\nu\rangle_{k}\lvert p_{i}\rangle\lvert
q_{j}\rangle\mapsto$ $\displaystyle\ \lvert b\rangle_{b}\lvert
i\rangle_{d}\lvert j\rangle_{e}\lvert 1\rangle_{i}\lvert\nu\rangle_{k}\lvert
p_{i}+\nu\rangle\lvert q_{j}-\nu\rangle$ $\displaystyle\mapsto$
$\displaystyle\ (-1)^{b\left([(p+\nu)\notin G]\lor[(q-\nu)\notin
G]\right)}\lvert b\rangle\lvert i\rangle\lvert j\rangle\lvert
1\rangle\lvert\nu\rangle\lvert p_{i}+\nu\rangle\lvert q_{j}-\nu\rangle$ (51)

for $U$ and $V$ respectively, as controlled by the register selecting them.
See Section II.2 and Section II.4 for a detailed description of the selection
operations and their circuit implementations.

We bound the error of our implementation and provide a complete costing of the
qubitization algorithm in Section II.5, with our main result summarized in
Theorem 4.

###  II.1 Circuit implementation of qubitization

Our aim in this section is to provide a high-level overview of how
qubitization can be realized to sample from the energy eigenvalues of a first
quantized Hamiltonian. The general strategy is to combine the circuits used
for state preparation and selection between $T$, $U$, and $V$ in order to
minimize the overall cost. In particular, we will be preparing superpositions
over $i$ and $j$, which are used independent of which of $T$, $U$, and $V$ is
being block encoded. Then we will need to perform operations on the momentum
given in registers $i$ and $j$, so we will use $i$ and $j$ to control the swap
of the momenta into temporary ancillae in order to perform the needed
operations. For the case of $T$ we are just computing a phase factor $\pm 1$
dependent on the value of $p$. In contrast, for $U$ we are performing
arithmetic on the momentum from register $j$, and for $V$ we are performing
arithmetic on the momenta from both registers $i$ and $j$. Therefore, the
operations $p+\nu$ and $p-\nu$ will be controlled on registers that select
between $T$, $U$, and $V$.

Our protocol for preparing the state (II) consists of the following steps.

  1. 1.

Rotate the qubit for selecting between $T$ and $U+V$.

  2. 2.

Prepare an equal superposition state over $\eta$ values of $i$ and $j$, and
check that they are not equal.

  3. 3.

Prepare superpositions over $w$, $r$ and $s$ that will be used for the $T$
part of the select operation. These preparations do not need to be controlled,
because it is only the $\mathrm{SEL}$ part that is controlled on the qubit for
$T$. The costing for this will be shown to be $2n_{p}+2b_{r}-7$ in the
subsection below for $T$.

  4. 4.

Prepare an appropriate superposition state for selecting between $U$ and $V$.
We will show that the complexity is approximately
$4\left\lceil\log\left(\eta+\lambda_{\zeta}\right)\right\rceil+2b_{r}-9$
Toffoli gates.

  5. 5.

Prepare a superposition over $\nu$ with weights proportional to $1/\|\nu\|$.
The costing for this will be described below for the $U$ and $V$ parts of the
Hamiltonian.

  6. 6.

Prepare a state with amplitudes proportional to $\sqrt{\zeta_{\ell}}$, which
may be done with complexity $\lambda_{\zeta}=\sum_{\ell}\zeta_{\ell}$.

Because we have written the state to be prepared as a tensor product of the
individual states, the overall cost of the state preparation is just the cost
of the preparation of the states on the individual subsystems.

In order to perform the $\mathrm{SEL}$ operation, the key parts are as
follows.

  1. 1.

Add $\nu$ into the momentum register at location $i$, controlled on the first
and sixth registers which select $V$.

  2. 2.

Subtract $\nu$ from the momentum register at location $j$, controlled on the
first qubit that is $\lvert 1\rangle$ for $U$ or $V$.

  3. 3.

Perform a $Z$ gate on the $\lvert+\rangle$ state if $(q-\nu)\notin G$ (for
$U$) or $[(p+\nu)\notin G]\lor[(q-\nu)\notin G]$ (for $V$).

  4. 4.

Perform a $Z$ gate on $\lvert+\rangle$ state controlled on $p_{w,r}p_{w,s}$
for the case of $T$.

  5. 5.

For the case of $U$ perform the phase rotation $-e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}$.

In addition, for our phase measurement procedure the $\mathrm{SEL}$ operation
should be self-inverse. It is easily seen that the component of the
$\mathrm{SEL}$ for $T$ is self-inverse. To account for the $\mathrm{SEL}$ for
$U$ and $V$ we can introduce an ancilla qubit that selects whether we add or
subtract $\nu$. Controlling addition versus subtraction is done with no
additional Toffoli complexity. This ancilla qubit is flipped each step, to
produce a self-inverse operation in the same way as in Eq. (8) of Berry  _et
al._ (2018).

The following discussion provides a statement of the block-encoding procedure
for the Born-Oppenheimer Hamiltonian.

######  Lemma 1.

The operator $(\mathrm{PREP}^{\dagger}\otimes
I)\mathrm{SEL}(\mathrm{PREP}\otimes I)$ forms, for
$N_{a}=10+2n_{\eta}+3n_{p}+2\lceil\log(n_{p}-1)\rceil+\lceil\log(L)\rceil$, an
$(\lambda,N_{a},0)$-block encoding of $H_{\rm BO}$, meaning that

$(\langle 0\rvert^{N_{a}}\otimes I)(\mathrm{PREP}^{\dagger}\otimes
I)\mathrm{SEL}(\mathrm{PREP}\otimes I)(\lvert 0\rangle^{N_{a}}\otimes
I)=\frac{\lvert 0\rangle\\!\langle 0\rvert^{\otimes N_{a}}\otimes H_{\rm
BO}}{\lambda}$

###### Proof.

The proof of this lemma follows by noting that the number of qubits needed to
encode the state yielded by $\mathrm{PREP}$ as given in (II). That coincides
with the ancillary space for block-encoding using the LCU lemma as in Berry
_et al._ (2015); Childs and Wiebe (2012). Further, the LCU lemma specifically
shows that if $\mathrm{PREP}^{\prime}\lvert
0\rangle^{N_{a}}=\sum_{\ell}\alpha_{\ell}\lvert j\rangle$ and
$\mathrm{SEL}^{\prime}\lvert\ell\rangle\lvert\psi\rangle=\lvert\ell\rangle
U_{\ell}\lvert\psi\rangle$ for $\alpha_{\ell}\geq 0$ then

$(\lvert 0\rangle\\!\langle 0\rvert^{\otimes N_{a}}\otimes
I)(\mathrm{PREP}^{\prime\dagger}\mathrm{SEL}^{\prime}\mathrm{PREP}^{\prime})\lvert
0\rangle^{N_{a}}\lvert\psi\rangle=\frac{\lvert
0\rangle^{N_{a}}\sum_{\ell}\alpha_{\ell}H_{\ell}\lvert\psi\rangle}{\sum_{\ell}\alpha_{\ell}}\,.$
(52)

The unitary operations $U_{\ell}$ are given implicitly in the above list for
$\mathrm{SEL}$. In this section we do not consider the complexity of
implementing the $\mathrm{SEL}$ operator and so only use the expression as a
definition of the action of the operator on input states.

Specifically, recall that up to a constant energy offset due to the nuclear
positions the Hamiltonian is $H=T+U+V$ where

$\displaystyle T$
$\displaystyle=\frac{\pi^{2}}{\Omega^{2/3}}\sum_{j=1}^{\eta}\sum_{w\in\\{x,y,z\\}}\sum_{r=0}^{n_{p}-2}\sum_{s=0}^{n_{p}-2}2^{r+s}\sum_{b\in\\{0,1\\}}\left(\sum_{p\in
G}(-1)^{b(p_{w,r}p_{w,s}\oplus 1)}\lvert p\rangle\\!\langle
p\rvert_{j}\right),$ (53) $\displaystyle U$ $\displaystyle=\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{2\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\sum_{j=1}^{\eta}\sum_{b\in\\{0,1\\}}\left(-e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}\sum_{q\in G}(-1)^{b[(q-\nu)\notin G]}\lvert q-\nu\rangle\\!\langle
q\rvert_{j}\right),$ $\displaystyle V$ $\displaystyle=\sum_{\nu\in
G_{0}}\frac{\pi}{\Omega\left\lVert k_{\nu}\right\rVert^{2}}\sum_{i\neq
j=1}^{\eta}\sum_{b\in\\{0,1\\}}\left(\sum_{p,q\in G}(-1)^{b([(p+\nu)\notin
G]\lor[(q-\nu)\notin G])}\lvert p+\nu\rangle\langle p\rvert_{i}\cdot\lvert
q-\nu\rangle\\!\langle q\rvert_{j}\right)\,,$

here $p_{w,r}$, for example, denotes the $r^{\rm th}$ bit of the $w^{\rm th}$
component of the momentum.

In this form, it is clear that the Hamiltonian operators in the LCU
decomposition satisfy for each $\ell$, $H_{\ell}\in\\{S_{T},S_{U},S_{V}\\}$
where each of these unitary operators are defined in (49), (II) and (II). Thus
we need to demonstrate two facts. First, that the operations in $\mathrm{SEL}$
selectively apply these transformations and that the coefficients of the
expansion match the square roots of those in the preceeding expression.

Consider steps $1$ and $2$ in the definition of $\mathrm{SEL}$. These steps
can be implemented using a controlled unitary adder. A unitary subtractor
subtracting the vector $\nu$ from a momentum register $\lvert q\rangle_{j}$
has the action

$\lvert\nu\rangle\lvert q\rangle_{j}\mapsto\lvert\nu\rangle\lvert
q-\nu\rangle_{j},$ (54)

This therefore has the correct action for implementing the subtraction needed
to implement $U$ and $V$.

There are two cases we need to consider. First, assume $q-\nu\in G$. In this
case the operator describing the transformation has the same action as
$\lvert\nu\rangle\\!\langle\nu\rvert\otimes\lvert p-\nu\rangle\\!\langle
p\rvert_{j}$. Now let us assume that $p-\nu\not\in G$. In this case we have an
overflow. However, since $\sum_{b\in\\{0,1\\}}(-1)^{b[(q-\nu)\not\in G]}=0$ in
(21) the action of the unitary adder need not be defined in this case and
hence an additional qubit is not needed in the register to address potential
overflow problems. (We will use that later to simplify the implementation.)
With this fact in mind, it is then straightforward to show that steps 1 and 2
block encode the $V$ and $U$ terms in the Hamiltonian.

To explain how this gives the required linear combination of unitaries, note
that the general procedure is to implement an operator of the form
$\sum_{\ell}\alpha_{\ell}H_{\ell}$, prepare an ancilla with amplitudes
proportional to $\sqrt{\alpha_{\ell}}$, then perform a controlled $H_{\ell}$
operation. The first step in $\mathrm{SEL}$ gives the addition $p+\nu$ in
momentum register $i$ required for $V$. Step 2 gives the subtraction $q-\nu$
required for $V$, the together with step 3 the complete operation needed for
$V$ is performed as described in (II). Step 2 also gives $q-\nu$ as required
for $U$, step 3 also gives the sign flip for $U$, and step 5 gives the phase
factor needed for $U$. Together they give the required operations for $U$ as
shown in (II). The correct amplitudes are obtained by the amplitudes
$1/\|\nu\|$ in (II), which gives the correct weights $1/\|\nu\|^{2}$ in the
Hamiltonian. The $\mathrm{SEL}$ operation for $T$ is obtained entirely by step
4, which gives the sign flip. The appropriate weights in the sum are obtained
by using the superpositions over $j,w,r,s$ in (II). With these definitions in
place we can formally state the following.

The block-encoding of $T$ is a little less obvious, but it ultimately follows
similarly from the definitions of the involved terms. Recalling from (25) that
$\lambda_{T}=6\eta\pi^{2}(2^{n_{p}-1}\\!-1)^{2}/\Omega^{2/3}$ and noting the
definition of $T$ given in (24) we see immediately that the theorem holds
after substituting the terms in the LCU decomposition in (24) into the $U_{j}$
in (52). ∎

For the rotation of the qubit to select between $T$ and $U+V$, we find that
error in the rotation angle of $\Delta\theta$ results in error in the energy
no more than $\lambda\Delta\theta$ (see Appendix D). As a result, if we take
the number of bits for this rotation to be $n_{T}$, then the error from this
source is bounded as

$\frac{\pi\lambda}{2^{n_{T}}}\leq\epsilon_{T}.$ (55)

The Toffoli complexity of performing this rotation is $n_{T}-3$. In some
applications it may be possible to reduce the complexity by selecting $\Omega$
appropriately. That is, the ratio $\lambda_{T}/(\lambda_{U}+\lambda_{V})$ will
depend on $\Omega$, so one can choose $\Omega$ to that the rotation may be
performed exactly with a relatively small number of bits.

In order to quantify the complexity of performing the state preparation for
selecting between $U$ and $V$, we can take account of the relative amplitudes
between the different parts. We find that

$\frac{\lambda_{U}}{\lambda_{V}}=\frac{2\lambda_{\zeta}}{\eta-1}\qquad\qquad\lambda_{\zeta}=\sum_{\ell=1}^{L}\zeta_{\ell}\,$
(56)

where we remind the reader that $\lambda_{\zeta}=\eta$ for charge neutral
systems. Even if that is not the case, we see that $\lambda_{U}/\lambda_{V}$
is simply a ratio of integers. Rather than performing a rotation, one can
instead achieve an exact state preparation as follows. First, prepare an equal
preparation over $\eta-1+2\lambda_{\zeta}$ numbers. Then, perform an
inequality test between this register and $\eta$ to give the correct relative
weightings $\eta-1$ and $2\lambda_{\zeta}$. There is a subtlety in that we
will initially prepare superpositions over all numbers $i$ and $j$ including
the case $i=j$. This effectively increases the value of $\lambda_{V}$, so the
ratio we want is $2\lambda_{\zeta}/\eta$. Then we prepare an equal
superposition over $\eta+2\lambda_{\zeta}$ numbers (rather than
$\eta-1+2\lambda_{\zeta}$).

After preparing the superposition state and applying the inequality test, use
a Toffoli between the result of the inequality test and the first qubit
selecting between $T$ and $U+V$, to give a qubit flagging that we are applying
$V$ (so will perform the addition $p+\nu$). Writing the number of qubits
needed to store $\eta+2\lambda_{\zeta}$ as $n_{\eta\zeta}$, the complexity of
preparing the equal superposition state is Lee  _et al._ (2021) (see Appendix
A.2, page 39 of version 2):

$3n_{\eta\zeta}+2b_{r}-9,$ (57)

and the cost of the inequality test is $n_{\eta\zeta}-1$. Here $b_{r}$ is a
number of bits of precision used for rotation on an ancilla. Together with the
extra Toffoli, the cost is

$4n_{\eta\zeta}+2b_{r}-9.$ (58)

The convenience of this is that it introduces little error due to the state
preparation.

The probability of success of the state preparation is usually above $0.999$
for $b_{r}=7$ (see Figure 10 of Sanders  _et al._ (2020)). The formula for the
probability of success is given in Eq. (205) of Sanders  _et al._ (2020), and
is

${\rm Ps}\left(n,b_{r}\right)=\frac{n}{2^{\lceil\log
n\rceil}}\left[\left(1+\left(2-\frac{4n}{2^{\lceil\log
n\rceil}}\right)\sin^{2}\theta(n,b_{r})\right)^{2}+\sin^{2}(2\theta(n,b_{r}))\right],$
(59)

where $n$ is the number of basis states to create an equal superposition over,
and $\theta(u,b_{r})$ is the angle of rotation. Since we are using $b_{r}$
bits for the rotation angle, the angle is of the form $2\pi n/2^{b_{r}}$ for
some $n$, and

$\theta(n,b_{r})=\frac{2\pi}{2^{b_{r}}}{\rm
round}\left(\frac{2^{b_{r}}}{2\pi}\arcsin\left(\sqrt{\frac{2^{\lceil\log
n\rceil}}{4n}}\right)\right).$ (60)

More generally, we have proven that the probability of error is no more than
$(3/2)^{2}(\pi/2^{b_{r}})^{2}$ with a slight adjustment to the formula for
$\theta$ (see Appendix J). A simpler case is that where
$\lambda_{\zeta}=\eta$, which occurs when the total charge is zero. In that
case, the ratio we want is just 2. Instead of needing to prepare an equal
superposition over $\eta+2\lambda_{\zeta}$ values, we only need a
superposition over 3. Including $n_{T}-3$ Toffolis to perform the rotation on
the first qubit, the total complexity for the preparation on the qubits
selecting $T$, $U$ and $V$ would be

$n_{T}+4n_{\eta\zeta}+2b_{r}-12.$ (61)

Next we consider the cost of the preparation of the equal superposition of $i$
and $j$ in unary. The steps needed are as follows.

  1. 1.

Prepare an equal superposition in binary, using the approach in Appendix A.2
of Lee  _et al._ (2021), with cost $3n_{\eta}+2b_{r}-9$ for each, where
$n_{\eta}=\lceil\log\eta\rceil$ and $b_{r}$ is the angle of rotation on an
ancilla qubit.

  2. 2.

Test if $i$ and $j$ are equal in binary, flagging the result in a qubit. This
qubit will need to be combined with another success qubit for applying the
potential, so the cost is $n_{\eta}$.

  3. 3.

The testing $i\neq j$ can be inverted with $n_{\eta}$ Toffolis.

  4. 4.

The inversion of the equal superposition of $i$ and $j$ can be inverted with
cost $3n_{\eta}+2b_{r}-9$ for each.

Adding all these costs together gives a total complexity

$14n_{\eta}+8b_{r}-36.$ (62)

For the failure case where $i=j$ _and_ we would otherwise be selecting $V$
(rather than $U$), we can simply combine this with the failure of the state
preparation for $\lvert\nu\rangle$, and perform $\mathrm{SEL}$ for $T$ in that
case.

For convenience of the rest of the procedure, we produce three qubits for
selecting $T$, $U$ and $V$. These are as follows.

  1. 1.

A qubit which is 0 if we are to perform $T$.

  2. 2.

A qubit which is 0 if we are performing the subtraction $q-\nu$. This will be
the case for $U$ or $V$.

  3. 3.

A qubit which is 0 for $V$ only, which controls the addition $p+\nu$.

First, there are four preparations of equal superposition states, and checking
that all of those succeeded takes 3 Toffolis. There are also another four
qubits needed to select $T,U,V$. These are the rotated qubit (register 1 in
(II)), the qubit resulting from the inequality test to select between $U$ and
$V$, the qubit flagging $i\neq j$, and the qubit flagging failure of the state
preparation for $\nu$. For simplicity of the explanation, let us call these
logic values $A,B,C,D$, respectively, with 0 corresponding to True.

To check that the state preparation for $U+V$ succeeded, we can use $(B\vee
C)\wedge D$, which takes two Toffolis (one for each logic operation). In one
scenario (OR), $T$ was selected as $A\vee\neg((B\vee C)\wedge D)$, and in the
other (AND) $T$ was selected as $A\wedge\neg((B\vee C)\wedge D)$. Either takes
one more Toffoli, but an extra Toffoli is needed to ensure that qubit is set
only for success of the preparation of the equal superposition states. That is
7 Toffolis so far.

Now, in the OR case, we set the second qubit (selecting $U$ or $V$) if
$\neg(A\vee\neg((B\vee C)\wedge D))$ and there were success of the equal
superposition preparations. That takes one more Toffoli. In the AND case we
would set that second qubit according to $(B\vee C)\wedge D$ and success of
the preparations, which is again one Toffoli. The third qubit can then be set
by taking the AND of $\neg B$ and the result in the second qubit. That is one
more Toffoli for a total of 9.

For the preparations of equal superposition states, the four needed are the
equal superposition over 3 basis states for $w$, one with
$\eta+2\lambda_{\zeta}$ basis states, and two with $\eta$. For the
superposition over three basis states required for $w$, it is highly efficient
to take $b_{r}=8$, so we use that rather than allowing arbitrary $b_{r}$
there. As a result, the block encoded state is multiplied by the probability

$P_{\rm eq}:={\rm Ps}\left(3,8\right){\rm
Ps}\left(\eta+2\lambda_{\zeta},b_{r}\right){\rm
Ps}\left(\eta,b_{r}\right)^{2}.$ (63)

This corresponds to giving an effective $\lambda$ divided by $P_{\rm eq}$.

###  II.2 Block-encoding the kinetic operator

Recall that the usual form of the kinetic operator is given by

$T=\sum_{j=1}^{\eta}\sum_{p\in G}\frac{\left\lVert
k_{p}\right\rVert^{2}}{2}\lvert p\rangle\\!\langle p\rvert_{j}.$ (64)

We have rewritten this expression to express the square explicitly as a sum of
products of bits, as

$T=\frac{2\pi^{2}}{\Omega^{2/3}}\sum_{j=1}^{\eta}\sum_{p\in
G}\sum_{w\in\\{x,y,z\\}}\sum_{r=0}^{n_{p}-2}\sum_{s=0}^{n_{p}-2}2^{r+s}p_{w,r}p_{w,s}\lvert
p\rangle\\!\langle p\rvert_{j},$ (65)

where $p_{w,r}$ is indicating bit $r$ of the $w$ component of $p$. To block-
encode this, we need to prepare the superposition state proportional to

$\lvert+\rangle\sum_{j=1}^{\eta}\lvert j\rangle\sum_{w=0}^{2}\lvert
w\rangle\sum_{r=0}^{n_{p}-2}2^{r/2}\lvert
r\rangle\sum_{s=0}^{n_{p}-2}2^{s/2}\lvert s\rangle.$ (66)

The preparation of the equal superposition over $j$ was addressed above. For
the equal superposition over $w$, we take $b_{r}=8$ and $n=2$, so the formula
$3n+2b_{r}-9$ for the complexity gives a cost of only $13$ Toffolis.

The superposition over $r$ and $s$ can be prepared with high chance of success
by using the following approach. First, we aim to prepare the following state
in unary

$2^{-n_{p}/2}\lvert n_{p}-1\rangle+\sum_{r=0}^{n_{p}-2}2^{-(r+1)/2}\lvert
r\rangle.$ (67)

Half the amplitude is on $\lvert 0\rangle$, then of the remaining states half
the amplitude is on $\lvert 1\rangle$, and so forth. This means that to
prepare the state in unary, we perform a Hadamard on the first qubit (labelled
0), then controlled on that we perform a Hadamard on the second qubit, and so
forth. There are in total $n_{p}-2$ controlled Hadamards, and each controlled
Hadamard can be performed with a single Toffoli and a catalytically used
$\lvert T\rangle$ state as in Figure 17 of Lee  _et al._ (2021). The qubit for
$n_{p}-1$ is being used to flag success of the state preparation (it flags
success with $\lvert 0\rangle$). Omitting the state $\lvert n_{p}-1\rangle$,
the state flagged for success is (up to normalization)

$\sum_{r=0}^{n_{p}-2}2^{-(r+1)/2}\lvert r\rangle.$ (68)

We then flip all these remaining bits so that the state becomes

$\sum_{r=0}^{n_{p}-2}2^{-(r+1)/2}\lvert
n_{p}-2-r\rangle=2^{-(n_{p}-1)/2}\sum_{r=0}^{n_{p}-2}2^{r/2}\lvert r\rangle.$
(69)

We will want the unary to be one-hot, but that conversion may be performed
with CNOTs so the complexity is omitted in the Toffoli count. Note that to
account for the case $r=0$, one can add one qubit in the state $\lvert
1\rangle$, so $r=0$ is encoded as $\lvert 10000\ldots\rangle$, $r=1$ is
encoded as $\lvert 11000\ldots\rangle$, and so forth. Then, the sequence of
CNOTs gives the desired one-hot unary. We therefore get a complexity of
$n_{p}-2$ for this state preparation. Together with the $13$ Toffolis used to
prepare the superposition over $w$, the overall complexity for the preparation
over the $w$, $r$, and $s$ registers is

$2n_{p}+9.$ (70)

The effect of the preparation having a small probability of failure is that
the effective value of $\lambda$ is increased to

$\lambda^{\prime}_{T}=\frac{6\eta\pi^{2}}{\Omega^{2/3}}2^{2(n_{p}-1)},$ (71)

where we have replaced $2^{n_{p}-1}-1$ with $2^{n_{p}-1}$.

In order to perform the $\mathrm{SEL}$ operation, we need to use the value of
$p$ in momentum register $j$. Because we need to perform arithmetic on the
momentum in register $j$ for $U$ and $V$ as well, the most efficient procedure
is to simply control a swap of the momentum in register $j$ into an ancilla
register. Similarly, the momentum in register $i$ will be swapped into an
ancilla register. The operations we perform on the momenta in these ancilla
registers will then be controlled by the first two qubits selecting between
$T,U,V$. That is more efficient than performing the operations in a controlled
way on every momentum register. The complexity of controlling the swaps of the
momenta out of and back into momentum registers $i$ and $j$ is $12\eta n_{p}$.
For the control, we perform unary iteration on $i$ and $j$. That requires
$\eta-2$ Toffolis for each of the four times (swapping the register out and
back for $i$, then again for $j$), for a total of $4\eta-8$. Together with
$12\eta n_{p}$ that gives a total of

$12\eta n_{p}+4\eta-8.$ (72)

We will consider that cost separately from the costs for implementing
$\mathrm{SEL}$ specifically for $T,U,V$. In the case of $T$, the other
operations needed to perform $\mathrm{SEL}$ are as follows.

  1. 1.

Use $w$ to control copying of component $w$ of $p$ into another ancilla. This
has cost $3(n_{p}-1)$.

  2. 2.

Copy bit $r$ of component $w$ of $p$ to an ancilla qubit, with cost $n_{p}-1$.

  3. 3.

Copy bit $s$ of component $w$ of $p$ to an ancilla qubit, with cost $n_{p}-1$.

  4. 4.

Controlled on those two bits, perform a phase flip on the $\lvert+\rangle$
state unless both bits are one. This gives an orthogonal $\lvert-\rangle$
state unless both bits are one. The orthogonal state gives no contribution to
the block-encoding. Only one Toffoli is needed here.

  5. 5.

Erase the two qubits containing $p_{w,r}$ and $p_{w,s}$, then the register
containing component $w$ of $p$. These erasures can be done with measurements
and Clifford gates.

The overall Toffoli complexity is only about as much as is needed to run
through the qubits representing $p$. It is $5(n_{p}-1)+1$. There is also some
complexity due to the fact that we also have a control qubit selecting $T$,
but control on one extra qubit needs only one more Toffoli. That gives a total
$\mathrm{SEL}$ cost for $T$

$5(n_{p}-1)+2,$ (73)

in addition to the $\mathrm{SEL}$ cost $12\eta n_{p}$ that is common between
$T,U,V$.

###  II.3 State preparation for $U$ and $V$

The $U$ and $V$ potential operators both have a $1/\|k_{\nu}\|^{2}$ weighting,
which means that in order to do the state preparation for the block encoding
we need to prepare a state proportional to

$\sum_{\nu\in G_{0}}\frac{1}{\left\lVert\nu\right\rVert}\lvert\nu\rangle,$
(74)

where $\nu\in
G_{0}=\big{[}-(N^{1/3}-1),N^{1/3}-1\big{]}^{3}\backslash\\{(0,0,0)\\}$ is a
three-component integer vector. In our implementation, we use eight registers
$\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle\lvert
m\rangle\lvert 0\rangle$ to hold this state. These subsystems are used as
follows.

  * •

$\lvert\mu\rangle$ flags which box was used in the series of nested boxes for
the state preparation.

  * •

$\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle$ are the
components of $\nu$ given as signed integers.

  * •

$\lvert m\rangle$ is an ancilla in an equal superposition used to give the
correct amplitude via an inequality test.

  * •

$\lvert 0\rangle$ flags that the state preparation was successful.

There will be four parts of the preparation where there is some amplitude for
failure. First, there is the preparation of $\mu$. Second, there is the
negative zero in the signed integer representation of $\nu$, which is not
allowed. Third, there is testing whether the position is inside the inner box
in the state preparation. Fourth, there is the success of the inequality test
used in the state preparation. Success for all four is indicated by $\lvert
0\rangle$ on the flag qubit above, which can then be used in amplitude
amplification.

The state that we will prepare has the form

$\frac{1}{\sqrt{\mathcal{M}2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\sum_{m=0}^{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil-1}\frac{1}{2^{\mu}}\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle\lvert
m\rangle\lvert 0\rangle+\lvert\Psi^{\bot}\rangle,$ (75)

where $B_{\mu}$ is a subset of $\nu$ to be defined later and
$\big{(}I\otimes\langle 111\rvert\big{)}\lvert\Psi^{\bot}\rangle=0$.
Conditioned on the qubit flagging success with $\lvert 0\rangle$, the
amplitude for each $\nu$ will then be

$\sqrt{\frac{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu}2^{n_{p}+2}}}\approx\frac{1}{8\sqrt{2^{n_{p}}}}\frac{1}{\left\lVert\nu\right\rVert}.$
(76)

The probability of the flag qubit in state $\lvert 0\rangle$ is around $1/4$,
and a single step of amplitude amplification boosts the success probability to
close to $1$ (Babbush  _et al._ , 2019, Fig. 1). Note that the probability is
slightly smaller here than in Babbush  _et al._ (2019), because we use a
preparation on $\mu$ with nonzero amplitude for failure.

The approach is to use a hierarchy of nested boxes in
$\big{[}{-}(N^{1/3}-1),N^{1/3}-1\big{]}$ indexed by $\mu$ and to prepare a set
of $\nu$ belonging to each box. We initially prepare a unary-encoded
superposition state

$\frac{1}{\sqrt{2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sqrt{2^{\mu}}\lvert\mu\rangle=\frac{1}{\sqrt{2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sqrt{2^{\mu}}\lvert
0\cdots 0\underbrace{1\cdots 1}_{\mu}\rangle.$ (77)

This state can be created using Hadamard and controlled Hadamard gates in a
similar way as we performed the state preparation for $r$ and $s$ in order to
implement the squaring of $p$. Here we have omitted the failure part of that
preparation, so this state is subnormalized, with the norm indicating the
probability of success. There is just one more bit than in the case where we
were preparing $r$ and $s$, so the cost is $n_{p}-1$ Toffolis.

In the next stage, we prepare

$\frac{1}{\sqrt{2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu_{x},\nu_{y},\nu_{z}=-(2^{\mu-1}-1),-0}^{2^{\mu-1}-1}\frac{1}{2^{\mu}}\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle,$
(78)

where the basis state is explicitly represented as

$\lvert 0\cdots 0\underbrace{1\cdots
1}_{\mu}\rangle\lvert\nu_{x,\text{sign}}0\cdots
0\nu_{x,\mu-2}\cdots\nu_{x,0}\rangle\lvert\nu_{y,\text{sign}}0\cdots
0\nu_{y,\mu-2}\cdots\nu_{y,0}\rangle\lvert\nu_{z,\text{sign}}0\cdots
0\nu_{z,\mu-2}\cdots\nu_{z,0}\rangle.$ (79)

Note that in this signed integer representation, the number zero appears
twice, with both a plus sign and a minus sign. The sum with “$,-0$” indicates
that the negative zero is included at this point. This state preparation can
be implemented using six Hadamards and $3(n_{p}-1)$ controlled Hadamards.
Specifically, we create the uniform superposition state in registers
$\lvert\nu_{x,\text{sign}},\nu_{x,0}\rangle\lvert\nu_{y,\text{sign}},\nu_{y,0}\rangle\lvert\nu_{z,\text{sign}},\nu_{z,0}\rangle$
using six Hadamards. We always have an equal superposition on $\nu_{x,0}$,
$\nu_{y,0}$ and $\nu_{z,0}$ because $\mu$ has a minimum value of $2$. Then we
generate uniform superpositions in $\lvert\nu_{x,j}\rangle$,
$\lvert\nu_{y,j}\rangle$, and $\lvert\nu_{z,j}\rangle$ controlled on the
($j+1$)th qubit in $\lvert\mu\rangle$ from the right for $j=2,\ldots,n_{p}$.
As before, the controlled Hadamards can be implemented with a Toffoli and a
catalytic state, so the Toffoli cost is $3(n_{p}-1)$.

Now, to remove the minus zero in this representation, we flag it as a failure.
This can be done by checking each $\lvert\nu_{x}\rangle$,
$\lvert\nu_{y}\rangle$, and $\lvert\nu_{z}\rangle$ whether the sign bit is one
and the remaining bits are zero, which requires controlled $\mathrm{Toffoli}$s
each having $n_{p}+1$ controls. We then use another two $\mathrm{Toffoli}$
gates to flag if any of these checkings returned true. Altogether, the Toffoli
cost is

$3n_{p}+2.$ (80)

In the next stage, we test whether all of $\nu_{x}$, $\nu_{y}$, and $\nu_{z}$
are smaller in absolute value than $2^{\mu-2}$. First, we convert the $\mu$
register to one-hot unary using $\mathrm{CNOT}$s. Then, conditioned on the
$\mu$th qubit being one, we flag $\lvert f_{\text{box}}\rangle$ depending on
whether the ($\mu-1$)th qubit $\lvert\nu_{x,\mu-2}\rangle$,
$\lvert\nu_{y,\mu-2}\rangle$, and $\lvert\nu_{z,\mu-2}\rangle$ are all zero,
using a controlled $\mathrm{Toffoli}$ gate with four controls, which has a
cost of three Toffolis. Because this must be done $n_{p}$ times, this gives a
cost of

$3n_{p}.$ (81)

The costs so far are then $4(n_{p}-1)$ controlled Hadamards and $6n_{p}+2$
Toffolis. Although the complexity of the controlled Hadamard is one Toffoli,
we will need to pay this cost again when inverting the preparation, but the
Toffolis can be erased with Clifford gates. The state excluding the failures
then becomes

$\frac{1}{\sqrt{2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\frac{1}{2^{\mu}}\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle,$
(82)

where $B_{\mu}$ is the set of $\nu$ such that $|\nu_{x}|$, $|\nu_{y}|$,
$|\nu_{z}|$ are less than $2^{\mu-1}$, excluding the case that they are all
less than $2^{\mu-2}$.

We now prepare an ancilla state in an equal superposition of $\lvert m\rangle$
for $m=0$ to $\mathcal{M}-1$. The parameter $\mathcal{M}$ is a power of two
and is chosen to be sufficiently large to obtain an accurate block-encoding.
This state can be prepared entirely using Hadamards, so no non-Clifford gates
are required. The resulting state excluding the failures is now

$\frac{1}{\sqrt{\mathcal{M}2^{n_{p}+2}}}\sum_{\mu=2}^{n_{p}+1}\sum_{\nu\in
B_{\mu}}\sum_{m=0}^{\mathcal{M}-1}\frac{1}{2^{\mu}}\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle\lvert
m\rangle.$ (83)

In the final stage, we test the inequality

$\big{(}2^{\mu-2}\big{)}^{2}\mathcal{M}>m\big{(}\nu_{x}^{2}+\nu_{y}^{2}+\nu_{z}^{2}\big{)}.$
(84)

We then need to use the four flag qubits for success to produce an overall
success flag. The final state will then be

$\frac{1}{\sqrt{\mathcal{M}2^{n_{p}+2}}}\sum_{\mu=2}^{n_{\mu}}\sum_{\nu\in
B_{\mu}}\sum_{m=0}^{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil-1}\frac{1}{2^{\mu}}\lvert\mu\rangle\lvert\nu_{x}\rangle\lvert\nu_{y}\rangle\lvert\nu_{z}\rangle\lvert
m\rangle\lvert 0\rangle+\lvert\Psi^{\bot}\rangle$ (85)

as desired. To achieve this, we compute
$\big{(}2^{\mu-2}\big{)}^{2}\mathcal{M}$ and
$m\big{(}\nu_{x}^{2}+\nu_{y}^{2}+\nu_{z}^{2}\big{)}$ separately in two ancilla
registers and compare the results. The unsigned number
$m\big{(}\nu_{x}^{2}+\nu_{y}^{2}+\nu_{z}^{2}\big{)}$ is bounded by
$3\mathcal{M}2^{2n_{\mu}-2}$, so we need at most $2n_{\mu}+n_{\mathcal{M}}$
qubits to represent it, where we use

$n_{\mathcal{M}}=\lceil\log(\mathcal{M})\rceil.$ (86)

The value $\big{(}2^{\mu-2}\big{)}^{2}\mathcal{M}$ is upper bounded by
$2^{2n_{\mu}-4}\mathcal{M}$, so we could represent it with
$2n_{\mu}-4+n_{\mathcal{M}}$ qubits.

Because $\mathcal{M}$ is a classically chosen value, and $\mu$ is given in
unary, $\big{(}2^{\mu-2}\big{)}^{2}\mathcal{M}$ can be computed in place with
_no_ gates, and simply a relabelling of qubits. We can place zeroed qubits in
between the qubits of the (one-hot) unary representation of $\mu$, to give
$2\mu$. For example, $\mu=1$ is changed from $000001$ to
$0\underline{0}0\underline{0}0\underline{0}0\underline{0}0\underline{0}1\underline{0}$,
$\mu=2$ is changed from $000010$ to
$0\underline{0}0\underline{0}0\underline{0}0\underline{0}1\underline{0}0\underline{0}$,
and so forth, where the underlined zeros are those added. Then this unary
representation of $2\mu$ is also a binary representation of $2^{2\mu+1}$.
Lastly, to multiply by $2^{-5}\mathcal{M}$, we simply pad with
$\log\mathcal{M}-5$ zeros. In practice these zeroed ancillae do not even need
to be included, and we simply relabel the qubits for $\mu$.

To compute $m\big{(}\nu_{x}^{2}+\nu_{y}^{2}+\nu_{z}^{2}\big{)}$, we first use
the result Lemma 8 given in Appendix G to cost the sum of three squares, which
is $3n^{2}-n-1$ for $n$-bit numbers. We can regard $\nu_{x}$, $\nu_{y}$, and
$\nu_{z}$ as unsigned numbers of $n_{p}$ bits, so squaring them costs
$3n_{p}^{2}-n_{p}-1$ $\mathrm{Toffoli}$s. Finally, we multiply two binary
numbers of length $\log(\mathcal{M})$ and $2n_{p}+2$ using a similar approach
(see Lemma 10) with

$2n_{\mathcal{M}}\big{(}2n_{p}+2\big{)}-n_{\mathcal{M}}$ (87)

$\mathrm{Toffoli}$ gates. To test the inequality
$\big{(}2^{\mu-2}\big{)}^{2}\mathcal{M}>m\big{(}\nu_{x}^{2}+\nu_{y}^{2}+\nu_{z}^{2}\big{)}$,
we use a ($2n_{p}+n_{\mathcal{M}}+2$)-bit comparator, which has cost Sanders
_et al._ (2019)

$2n_{p}+n_{\mathcal{M}}+2.$ (88)

The net cost of the squaring, multiplication and inequality test is then

$3n_{p}^{2}+n_{p}+1+4n_{\mathcal{M}}\big{(}n_{p}+1\big{)}$ (89)

$\mathrm{Toffoli}$ gates. There are a further 3 Toffolis needed to produce the
overall success flag qubit.

Now for the overall block encoding, the state preparation needs to be
inverted. Rather than paying this cost again, we can keep the ancillae used to
perform the arithmetic and inequality test. As a result, the inequality test
and arithmetic can be inverted with Clifford gates. The 3 Toffolis used for
the overall success flag qubit can be inverted with Clifford gates too.
Similarly, the $6n_{p}+2$ Toffolis used for the testing on $\nu$ can be
inverted with zero non-Clifford cost. The only operations in this state
preparation that have a non-Clifford cost in the inverse preparation are the
$4(n_{p}-1)$ controlled Hadamards used for the preparation of the $\mu$
register and superpositions on the $\nu$ register. As a result, the total
Toffoli cost for the preparation and inverse preparation is

$3n_{p}^{2}+15n_{p}-7+4n_{\mathcal{M}}\big{(}n_{p}+1\big{)}.$ (90)

In the case where amplitude amplification is performed, the cost will be three
times this, due to the need to invert the inequality test and calculation,
then perform it again.

The other key parts of the state preparation for $U$ and $V$ are to prepare
the superpositions over $i$, $j$ and $\ell$. The preparation over $i$ and $j$
are just simple preparations of equal superposition states which we have
addressed already. The final preparation we need to consider here is therefore
the preparation over $\ell$ with weightings $\sqrt{\zeta_{\ell}}$. In the
preparation over the second qubit selecting between $U$ and $V$, we have
already prepared a superposition, which will be an equal superposition over
$2\lambda_{\zeta}$ values for the case of $U$. That is the case where we need
to prepare this superposition state, so we can just prepare it in this case,
rather than preparing it for all cases as described above.

For the preparation, we can just use a QROM on this register, outputting the
value of $\ell$ in an output register. When there are $\zeta_{\ell}$ values in
the input register that give the same value of $\ell$ in the output register,
the result is obtaining an amplitude exactly proportional to
$\sqrt{\zeta_{\ell}}$. Moreover, we will need to compute $k_{\nu}\cdot
R_{\ell}$ in performing the $\mathrm{SEL}$ operation. The value of $R_{\ell}$
can be output in another ancilla register at the same time. The complexity of
the QROM, using the simple form of QROM, is just $\lambda_{\zeta}-1$. The QROM
should be controlled on the qubit which $U+V$, so $R_{\ell}$ are not output
there. The QROM automatically does not output $R_{\ell}$ for $V$, because this
input register is outside the range for $V$.

Although we have aimed to prepare a state with amplitudes
$\sqrt{\zeta_{\ell}}$, note that the only way that $\ell$ is used is in
outputting $R_{\ell}$. That means, if we are using QROM to output both $\ell$
and $R_{\ell}$ at once, then the $\ell$ register is not actually used. That
means we can simply omit the $\ell$ register, and output $R_{\ell}$. That
saves ancilla qubits, though does not change the Toffoli count.

In inverting this preparation, erasing the output of the QROM has a Toffoli
complexity of $2^{k}+\left\lceil 2^{-k}\lambda_{\zeta}\right\rceil$ where $k$
is taken to be a power of 2. The resulting complexity is approximately the
square root of the initial QROM. For convenience of describing the
complexities, we will define

${\rm Er}\left(x\right):=\min_{k\in\mathbb{N}}\left(2^{k}+\left\lceil
2^{-k}x\right\rceil\right),$ (91)

so this cost is ${\rm Er}\left(\lambda_{\zeta}\right)$. In terms of this
function, the cost of outputting the values of $R_{\ell}$ via QROM and erasing
them again can be given as

$\lambda_{\zeta}+{\rm Er}\left(\lambda_{\zeta}\right).$ (92)

###  II.4 The $\mathrm{SEL}$ operations for $U$ and $V$

We have described the block encoding for $T$, and the state preparation for
$U$ and $V$, and now give the method to perform the $\mathrm{SEL}$ operations
for $U$ and $V$. The general principle is as follows.

  1. 1.

Use the $i$ and $j$ registers to control a swap of the momentum registers $p$
and $q$ into ancillae.

  2. 2.

Controlled on the qubit that is $\lvert 0\rangle$ only for $V$, add $\nu$ (in
place) into the ancilla for $p$.

  3. 3.

Controlled on the qubit that is $\lvert 0\rangle$ for $U$ or $V$, subtract
$\nu$ into the ancilla for $q$.

  4. 4.

Check if either $p+\nu$ or $p-\nu$ is outside the box $G$, and if so perform a
phase gate on the $\lvert+\rangle$ state.

  5. 5.

Control swap the ancillae back into the momentum registers.

  6. 6.

Apply the phase $e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$ in the case of $U$
only.

The cost of the controlled swaps was discussed above, and is $12\eta n_{p}$.
The key cost to discuss now is the addition and subtraction of $\nu$ to
momenta. The difficulty is that the momentum and $\nu$ are encoded as signed
integers, whereas the usual methods for addition/subtraction are for integers
encoded in two’s complement.

In order to achieve this, just one of the numbers needs to be converted into
two’s complement. For example, since $p+\nu$ is computed in-place, we can
convert $p$ to two’s complement, then use the sign bit of $\nu$ to control
adding or subtracting $|\nu|$ into $p$. This control of addition versus
subtraction can be done with no more complexity than regular addition and
subtraction. This is done for each component, so the overall complexity of
converting between signed and two’s complement is approximately $12$ times
$n_{p}$ for the conversions. The addition and subtraction also need to be
controlled by the qubits selecting between $T,U,V$, so have complexity twice
uncontrolled additions and subtractions. That gives an overall complexity
about

$24n_{p}.$ (93)

More specifically, to convert an $n$ digit number from signed to two’s
complement, we take the sign qubit and use it as a control on CNOTs with the
remaining qubits. Then we add it to the remaining $n-1$ bits with complexity
$n-2$ Toffolis. We can just invert this process for converting in the opposite
direction. Because $p$ and $q$ have $n_{p}$ bits in each of their three
components, the complexity of converting them to two’s complement is
$6(n_{p}-2)$.

For the cost of the controlled addition (or subtraction), there is cost
$n_{p}+1$ to make each component controlled, which can be achieved by (for
example) controlling copying out of $\nu$ to an ancilla. (That can be erased
without Toffolis by measurements and phase fixups.) Then the addition has cost
$n_{p}+1$, giving a total cost of $2(n_{p}+1)$ for the controlled
addition/subtraction of each component. The total cost is $12(n_{p}+1)$ for
these controlled additions and subtractions.

Lastly we need to convert $p+\nu$ and $q-\nu$ back to signed form.
Adding/subtracting $\nu$ expands these numbers by 2 bits, so the cost is now
$6n_{p}$. The total cost is then

$6(n_{p}-2)+12(n_{p}+1)+6n_{p}=24n_{p}$ (94)

Toffolis. Thus we get exactly the same cost as approximately estimated above.

Another major part of the $\mathrm{SEL}$ operations is to check whether
$p+\nu$ or $q-\nu$ is outside the bounding box $G$. Because the complexity
only depends on $N$ through the number of bits used in the arithmetic $n_{p}$,
it is best to take $N$ as large as possible given $n_{p}$, which means that
the bounds of the box $\pm(N^{1/3}-1)/2$ are taken to be of the form
$2^{n_{p}-1}-1$. That means that, if $p+\nu$ or $q-\nu$ is outside the
bounding box, then at least one of the extra qubits introduced in the
arithmetic is set to $\lvert 1\rangle$. These qubits are additional, so are
not swapped from the ancilla back to the main momentum registers.

If we perform the checking as presented above, then we use these qubits to
control a $Z$ on the $\lvert+\rangle$ state, in order to remove parts of the
Hamiltonian where $p+\nu$ or $q-\nu$ are outside the bounding box, which are
excluded. However, note that in the block encoding we would select on $\lvert
0\rangle$ for these ancillae, which means that these parts of the Hamiltonian
are automatically eliminated. It is therefore unnecessary to perform this
controlled phase in the $\mathrm{SEL}$ operation.

The last thing we need to do for implementing $\mathrm{SEL}$ is to apply the
phase factor $-e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$. To achieve this we need
to compute $k_{\nu}\cdot R_{\ell}$ then add it into a phase gradient register
to apply the phase. We do not need to make the addition into the phase
gradient register controlled on the qubits that select $U$, because the
register with $R_{\ell}$ is nonzero only for $U$. The minus sign in this phase
does need to be performed in a controlled way, but that is just a controlled
phase on two qubits which is a Clifford gate.

In order to perform the phase shift, a subtlety is that it is given in terms
of $k_{\nu}$ rather than $\nu$. However, the conversion factor between
$k_{\nu}$ and $\nu$ can simply be absorbed into $R_{\ell}$, rather than
explicitly multiplied. This multiplication is most conveniently performed with
signed integers. The product of each pair of components will then have a sign
that is the product of the signs of the individual components. That product of
signs can be computed with a CNOT. However, if we were to then add the
components to compute $k_{\nu}\cdot R_{\ell}$, we would need to deal with the
problem of adding signed integers again. Instead, it is convenient to add the
products of each pair of components into the phase gradient state separately.
They can be added or subtracted into the phase gradient state controlled on
the sign with no more complexity than a regular addition.

The maximum value of a component of $k_{\nu}$ is $2\pi/\Omega^{1/3}$ times
$N^{1/3}-1$, which is the maximum value of a component of $\nu$, and the
maximum value of a component of $R_{\ell}$ is $\Omega^{1/3}/2$, since
$\Omega^{1/3}$ is the width of the spatial region. Therefore, a component of
$k_{\nu}$ times $R_{\ell}$ can be rewritten as a component of $\nu$ times
$2\pi R_{\ell}/\Omega^{1/3}$, which has a maximum value of $\pi$. There is one
initial CNOT used to find the product of the signs of the component of $\nu$
and $R_{\ell}$. That will be used to control addition or subtraction of
numbers into the phase gradient state. For the multiplication of the least
significant bit of $\nu$ times $2\pi R_{\ell}/\Omega^{1/3}$, we have cost
$n_{R}-1$ to compute the bit products, because we now exclude the sign bit.
The most significant bit of that product being added into the phase gradient
state corresponds to a phase shift of $\pi/2$, and the least significant bit
corresponds to a phase shift of $2\pi/2^{n_{R}}$. As a result, the cost of
adding these bits into the phase gradient state is $n_{R}-2$. For the next
most significant bit of $\nu$, the most significant bit of the product would
give phases of $\pi$, so we still need to compute the $n_{R}-1$ bit products,
but since all phases are multiplied by 2, the cost of adding into the phase
gradient state is reduced to $n_{R}-3$. After that, there is cost $n_{R}-2$
for the bit products, and $n_{R}-4$ for addition into the phase gradient
state. For the $k$’th bit of $\nu$, we have cost $n_{R}-k+1$ for the bit
products and $n_{R}-k-1$ for the addition for a total of $2(n_{R}-k)$. There
is a subtlety in that if $n_{R}\leq n_{p}$ then $n_{R}-k-1$ can be negative.
Ignoring that case for the moment, because we would typically expect
$n_{R}>n_{p}$, the Toffoli cost is

$-1+\sum_{k=1}^{n_{p}}2(n_{R}-k)=2n_{p}n_{R}-n_{p}(n_{p}+1)-1.$ (95)

In the case where $n_{R}=n_{p}$, then we get one more Toffoli, because we
simply have $n_{R}-k$ instead of $n_{R}-k-1$ for the maximum value of $k$.
Then the formula for the complexity can be given as
$2n_{p}n_{R}-n_{p}(n_{p}+1)=n_{R}(n_{R}-1)$. In the case where $n_{R}<n_{p}$,
then at $k=n_{R}+1$ we have both costs equal to zero, and so the sum can be
taken to $k=n_{R}$, to give

$\sum_{k=1}^{n_{R}}2(n_{R}-k)=n_{R}(n_{R}-1).$ (96)

Here we have removed the $-1$ at the beginning because the expression in the
sum overcounts the Toffolis by 1 at $k=n_{R}$. This cost is multiplied by 3,
because it is done once for each component. Therefore the Toffoli cost may be
given as

$\begin{cases}3[2n_{p}n_{R}-n_{p}(n_{p}+1)-1],&n_{R}>n_{p}\\\
3n_{R}(n_{R}-1),&n_{R}\leq n_{p}.\end{cases}$ (97)

For convenience of the final expression, we can give the Toffoli cost as just
$6n_{p}n_{R}$, but the cost can be computed more accurately with the
expression (97). The maximum number of ancilla qubits used is $n_{R}-1$ for
the bit products, plus $n_{R}-3$ for the addition into the phase gradient
register, for a total of $2n_{R}-4$.

The total Toffoli costs for the block encoding of the Hamiltonian are now as
in Table 2. Here $b_{r}$ is a number of bits used for ancilla qubit rotations
in preparations of equal superposition states, and can be taken to be $7$ for
a probability of success no less than $0.999$.

Procedure |  Toffoli cost
---|---
Preparing the superposition on the registers selecting between $T,U,V$, where $n_{\eta\zeta}$ is the number of bits needed to represent $\eta+\lambda_{\zeta}$; see (61). |  $2(n_{T}+4n_{\eta\zeta}+2b_{r}-12)$
Preparing equal superpositions over $\eta$ values of $i$ and $j$ in unary; see (62). |  $14n_{\eta}+8b_{r}-36$
The state preparation cost for the $w,r,s$ registers used for $T$, including the factor of 2 for inversion; see (70). |  $2(2n_{p}+9)$
Controlled swaps of the $p$ and $q$ registers into and out of ancillae (which is used for $T$, $U$, and $V$); see (72). |  $12\eta n_{p}+4\eta-8$
The $\mathrm{SEL}$ cost for $T$; see (73). |  $5(n_{p}-1)+2$
The preparation of the state with amplitude $1/{\left\lVert\nu\right\rVert}$ and inversion; see (90). |  $3n_{p}^{2}+15n_{p}-7+4n_{\mathcal{M}}\big{(}n_{p}+1\big{)}$
The QROM for $R_{\ell}$, combined with the state preparation with amplitudes $\sqrt{\zeta_{\ell}}$ (which may be implicit); see (92). |  $\lambda_{\zeta}+{\rm Er}\left(\lambda_{\zeta}\right)$
Additions and subtractions of $\nu$ into the momentum registers; see (93). |  $24n_{p}$
Phasing by $-e^{-\mathrm{i}k_{\nu}\cdot R_{\ell}}$; see (97). |  $6n_{p}n_{R}$.
Table 2: The costs involved in the block encoding of the Hamiltonian for
qubitization.

The complete step needed for the qubitization, to give the step that one would
perform phase estimation on, needs a reflection on the control ancilla. The
cost of this reflection corresponds to the numbers of qubits used in the state
preparation. The total qubit cost, and the qubits needed to be reflected upon,
are given in Appendix C. There are a total of

$n_{\eta\zeta}+2n_{\eta}+6n_{p}+n_{\mathcal{M}}+16$ (98)

qubits to reflect on. The corresponding Toffoli cost is two less. However,
this reflection needs to be controlled on the qubits used as control for the
phase measurement. Also, there is a single Toffoli cost for iterating those
control qubits. We bundle those extra two Toffolis together with this cost,
giving the cost as in (98).

###  II.5 Total cost of constructing the qubitization operator

This is a complete accounting for the cost of a single step of the block
encoding, except that we need to estimate the error due to the finite values
of $\cal{M}$ and $n_{R}$ in order to determine the appropriate values to use
for these quantities. The finite value of $n_{R}$ corresponds to replacing
$R_{\ell}$ with an $n_{R}$-digit approximation $\widetilde{R}_{\ell}$, so that
the potential term $U$ that the algorithm actually implements is close to

$U=-\frac{4\pi}{\Omega}\sum_{\ell=1}^{L}\sum_{j=1}^{\eta}\sum_{\begin{subarray}{c}p,q\in
G\\\ p\neq
q\end{subarray}}\bigg{(}\zeta_{\ell}\frac{e^{ik_{q-p}\cdot\tilde{R}_{\ell}}}{\left\lVert
k_{p-q}\right\rVert^{2}}\bigg{)}\lvert p\rangle\langle q\rvert_{j}.$ (99)

By choosing $n_{R}$ sufficiently large, this error can be at most an
appropriate fraction of the overall error budget $\epsilon$. To bound the
error due to this source and choose the value of $n_{R}$, we use the following
Lemma.

######  Lemma 2.

Given $\|\widetilde{R}_{\ell}-R_{\ell}\|\leq\delta_{R}$ for all $\ell$, if
$\lambda_{j}$ is the $j^{\rm th}$ eigenvalue of the true Hamiltonian and (99)
is used to approximate the nuclear potential, then there exists an eigenvalue
of the approximated Hamiltonian $\widetilde{\lambda}_{j}$ such that

$|\lambda_{j}-\widetilde{\lambda_{j}}|\leq\frac{2\delta_{R}\eta\lambda_{\zeta}}{\Omega^{2/3}}\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert}.$ (100)

###### Proof.

First, note from (Horn and Johnson, 2012, Corollary 6.3.4) that for any normal
matrices $M$ and $M^{\prime}$ that $\|M-M^{\prime}\|\leq\epsilon$ implies that
for every eigenvalue $\lambda_{j}$ of $M$ that there exists an eigenvalue of
$M^{\prime}$ within distance $\epsilon$. Thus it suffices to bound the
spectral norm of the difference between the true potential operator and the
approximated potential. Our claim then follows from the following bound on the
difference between the two Hamiltonians.

$\displaystyle\left\lVert U-\widetilde{U}\right\rVert$
$\displaystyle=\left\lVert\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{4\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\sum_{j=1}^{\eta}\big{(}e^{-\mathrm{i}k_{\nu}\cdot\tilde{R}_{\ell}}-e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}\big{)}\bigg{(}\sum_{p\in G}\lvert p-\nu\rangle\\!\langle
p\rvert_{j}\bigg{)}\right\rVert$ (101) $\displaystyle\leq\eta\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{4\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\big{|}e^{-\mathrm{i}k_{\nu}\cdot\tilde{R}_{\ell}}-e^{-\mathrm{i}k_{\nu}\cdot
R_{\ell}}\big{|}\leq\eta\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{4\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\big{|}k_{\nu}\cdot(\tilde{R}_{\ell}-R_{\ell})\big{|}$
(102) $\displaystyle\leq\eta\sum_{\nu\in
G_{0}}\sum_{\ell=1}^{L}\frac{4\pi\zeta_{\ell}}{\Omega\left\lVert
k_{\nu}\right\rVert^{2}}\left\lVert
k_{\nu}\right\rVert\left\lVert\tilde{R}_{\ell}-R_{\ell}\right\rVert\leq\frac{2\delta_{R}\eta\lambda_{\zeta}}{\Omega^{2/3}}\sum_{\nu\in
G_{0}}\frac{1}{\left\lVert\nu\right\rVert}.$ (103)

Note that $\sum_{\nu\in G_{0}}\frac{1}{\left\lVert\nu\right\rVert}$ can be
further upper bounded by a function of $N$. However, we use the current bound
as it is tighter and more versatile for numerical implementation. ∎

Next we consider the error due to the finite value of $\cal{M}$. As can be
seen from (85), the amplitude for each $\nu$ is proportional to

$\sqrt{\frac{\lceil{\mathcal{M}(2^{\mu-2}/\left\lVert\nu\right\rVert)^{2}}\rceil}{\mathcal{M}2^{2\mu}2^{n_{p}+2}}},$
(104)

where $\mu$ is chosen such that $\nu\in B_{\mu}$. From the definition of
$B_{\mu}$, this means that $\mu$ is the minimum integer such that $|\nu_{x}|$,
$|\nu_{y}|$ and $|\nu_{z}|$ are all less than $2^{\mu-1}$. That is,

$\mu=\lfloor\log\max(|\nu_{x}|,|\nu_{y}|,|\nu_{z}|)\rfloor+2.$ (105)

As a result, the success probability of this state preparation is given by