# Quantum simulation of battery materials using ionic pseudopotentials

Modjtaba Shokrian Zini  Xanadu, Toronto, ON, M5G 2C8, Canada  Alain Delgado
Xanadu, Toronto, ON, M5G 2C8, Canada  Roberto dos Reis  Xanadu, Toronto, ON,
M5G 2C8, Canada  Pablo A. M. Casares  Xanadu, Toronto, ON, M5G 2C8, Canada
Jonathan E. Mueller  Volkswagen AG, Berliner Ring 2, 38440 Wolfsburg, Germany
Arne-Christian Voigt  Volkswagen AG, Berliner Ring 2, 38440 Wolfsburg, Germany
Juan Miguel Arrazola  Xanadu, Toronto, ON, M5G 2C8, Canada

###### Abstract

Ionic pseudopotentials are widely used in classical simulations of materials
to model the effective potential due to the nucleus and the core electrons.
Modeling fewer electrons explicitly results in a reduction in the number of
plane waves needed to accurately represent the states of a system. In this
work, we introduce a quantum algorithm that uses pseudopotentials to reduce
the cost of simulating periodic materials on a quantum computer. We use a
qubitization-based quantum phase estimation algorithm that employs a first-
quantization representation of the Hamiltonian in a plane-wave basis. We
address the challenge of incorporating the complexity of pseudopotentials into
quantum simulations by developing highly-optimized compilation strategies for
the qubitization of the Hamiltonian. This includes a linear combination of
unitaries decomposition that leverages the form of separable pseudopotentials.
Our strategies make use of quantum read-only memory subroutines as a more
efficient alternative to quantum arithmetic. We estimate the computational
cost of applying our algorithm to simulating lithium-excess cathode materials
for batteries, where more accurate simulations are needed to inform strategies
for gaining reversible access to the excess capacity they offer. We estimate
the number of qubits and Toffoli gates required to perform sufficiently
accurate simulations with our algorithm for three materials: lithium manganese
oxide, lithium nickel-manganese oxide, and lithium manganese oxyfluoride. Our
optimized compilation strategies result in a pseudopotential-based quantum
algorithm with a total runtime four orders of magnitude lower than the
previous state of the art for a fixed target accuracy.

##  I Introduction

Quantum computing is being actively studied as a potential method to support
the development of next-generation lithium-ion batteries [1, 2, 3, 4, 5, 6].
The driving motivation is that quantum algorithms are uniquely positioned to
perform highly-accurate simulations without incurring prohibitive
computational costs [7]. Nevertheless, considerable progress still needs to
occur on both hardware and algorithms to make this promise a reality. The main
theoretical challenges are the high cost of implementing quantum algorithms
and the difficulty of identifying the applications.

There has been considerable progress in developing quantum algorithms for
simulating the properties of molecules and materials. A variety of different
strategies have been proposed, ranging from variational approaches designed
for noisy hardware with few qubits [8, 9, 10, 4, 11], algorithms tailored for
early fault-tolerant quantum computers [12, 13, 14, 15], and variants of
quantum phase estimation that require the full capabilities of large-scale
fault-tolerant quantum computers [16, 17, 18, 19, 20, 21]. Particular
attention has been devoted to improving the efficiency of quantum phase
estimation and Hamiltonian simulation algorithms, which has led to an overall
cost reduction of several orders of magnitude [17]. This progress has been
fueled by innovations such as qubitization [22, 23, 24], Hamiltonian
factorization techniques [18, 20, 25, 21], interaction-picture simulations
[26, 27, 28], improved Trotter bounds [29], and first quantization methods
[30, 19].

Simulating bulk materials presents additional challenges beyond those
associated with simulating finite molecules [31]. Arguably, among the most
pressing ones is the frequent need to use large unit cells (supercells). For
example, in the context of lithium-ion batteries, large supercells are needed
to predict the most stable phases of cathode materials. Their energies
calculated for different values of the lithium-ion concentration determine the
voltage profile of the battery cell [32, 33, 34]. The size of the supercell is
even more critical for the simulation of chemical reactions at the electrode-
electrolyte interface [35]. This results in systems with many hundreds of
electrons requiring a very large number of plane-wave basis functions to
achieve high-accuracy simulations [36].

In this work, we introduce a quantum algorithm that uses ionic
pseudopotentials (PPs) to reduce the cost of using quantum phase estimation to
simulate material, mirroring a strategy widely used for density functional
theory simulations [37, 38]. Replacing the bare Coulomb potential due to the
nucleus and the core electrons with an effective potential leads to a
substantial reduction in both the number of electrons and plane waves needed
to accurately represent the system. However, this is accomplished at the price
of a more complicated pseudopotential operator describing the interaction
between the valence electrons and the effective ionic cores. This greatly
complicates the implementation of the resulting quantum algorithm and can
negate the benefits of reducing the number of electrons and plane waves.

We tackle this challenge by deriving highly-optimized implementation
strategies for qubitization-based quantum phase estimation in first
quantization. We focus on the Hartwigsen-Goedecker-Hutter (HGH)
pseudopotentials [39] and develop a tailored linear combination of unitaries
decomposition for the pseudopotential term of the Hamiltonian. We also
carefully engineer the compilation of the qubitization encoding to reduce the
implementation cost. To avoid costly quantum arithmetic, our methods make
frequent use of quantum read-only memory (QROM) subroutines. To the best of
our knowledge, this is the first example of a quantum algorithm that can
incorporate pseudopotentials. Overall, simulations with our algorithm requires
orders-of-magnitude fewer plane waves to reach convergence for typical
materials than comparable all-electron calculations. This results in circuit-
depth reductions of many orders of magnitude for a given target accuracy.

However, while it is well-understood how one might use quantum computers to
simulate key properties of batteries, such as equilibrium voltages [5], it has
not yet been established which concrete class of battery simulation would
benefit the most from the capabilities of quantum computers. After all,
developing better lithium-ion batteries requires solving a multitude of
problems, and scientists are already equipped with sophisticated simulation
techniques that can be run on powerful supercomputers. It is therefore crucial
to identify problems where the limitations of classical methods are most
severe, and whose solution would be most impactful to battery development.

We propose an application of quantum computers for batteries that we contend
meet these criteria: simulating lithium-excess cathode materials [40]. These
materials offer an avenue to dramatically increase the energy density of
state-of-the-art cathodes. With theoretical capacities that are roughly twice
as high as commercial batteries, they would make the driving range and cost of
electric vehicles competitive with internal combustion engines [41].
Unfortunately, lithium-excess materials suffer from substantial capacity loss
even after a single charging cycle [42]. The rapid degradation of lithium-
excess materials responsible for the capacity loss has been attributed to
irreversible structural transformations of the material, but the relationship
between the proposed redox mechanisms and the observed transformations remains
a topic of active discussion as these processes are difficult to probe
experimentally [40]. To better understand these processes, researchers rely on
density functional theory simulations. However, the available density
functionals are not accurate enough to single out the dominant mechanisms
driving the materials structural changes [40]. This makes it difficult to
develop solutions for capacity loss in lithium-excess cathode materials using
classical computer simulations.

To assess the potential for quantum computers to provide a solution to this
problem, we perform a detailed estimation of the resources required to
implement our quantum algorithm as applied to three lithium-excess cathode
materials: lithium manganese oxide ($\text{Li}_{2}\text{MnO}_{3}$), lithium
nickel-manganese oxide
($\text{Li}[\text{Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}]\text{O}_{2}$),
and lithium manganese oxyfluoride ($\text{Li}_{0.75}\text{MnO}_{2}\text{F}$)
[43, 44, 34]. In each of these cases, the use of ionic pseudopotentials allows
a reduction of the number of electrons by a factor of two, and a reduction in
the number of plane waves by roughly three orders of magnitude compared to
all-electron simulations for a fixed target accuracy. This results in a total
runtime for the algorithm that is about four orders of magnitude lower than
the previous state-of-the-art [19].

The rest of this work is organized as follows. Section II provides background
information on first-quantization quantum algorithms for electronic structure,
the theory of ionic pseudopotentials, and the basic properties of quantum
read-only memories. Our quantum algorithm is described in Section III,
outlining the linear combination of unitaries and qubitization subroutines
that constitute the main technical contribution of this work. This is
complemented with a detailed error analysis in Section IV and a calculation of
the qubit and gate cost of the full algorithm in Section V. We then study the
application of the algorithm to the simulation of lithium-excess cathode
materials in Section VI.

##  II Background

This work is an interdisciplinary effort covering topics across computational
chemistry, quantum computing, and lithium-ion batteries. In an effort to make
this manuscript self-contained, this section provides background information
on key concepts that will be used throughout.

###  II.1 The plane-wave electronic Hamiltonian in first quantization

The ultimate goal of the quantum algorithm presented here is to solve the
electronic structure problem:

$H\Psi_{0}(\bm{r}_{1},\dots,\bm{r}_{\eta})=E_{0}\Psi_{0}(\bm{r}_{1},\dots,\bm{r}_{\eta}),$
(1)

where $H$ is the Hamiltonian of $\eta$ interacting electrons, and $\Psi_{0}$
and $E_{0}$ are the ground-state wave function and energy, respectively.

In the Born-Oppenheimer approximation [45], the electronic Hamiltonian, $H$,
is given by

$H=T+U+V,$ (2)

where $T$ is the total kinetic energy operator, $U$ is the Coulomb potential
due to the nuclei and $V$ is the electron-electron interaction term [46].

Plane-wave functions are a natural basis set to represent the electronic
states in periodic materials. They can be used to encode the translational
symmetry of crystal structures and allow us to derive closed-form expressions
for the Hamiltonian matrix elements. Plane-wave functions are defined as

$\varphi_{p}(\bm{r})=\frac{1}{\sqrt{\Omega}}e^{i\bm{G}_{p}\cdot\bm{r}},$ (3)

where $\Omega$ is the volume of the material’s unit cell and $\bm{G}_{p}$ is
the reciprocal lattice vector

$\bm{G}_{p}=\left[\sum_{i=1}^{3}p_{i}b_{i_{x}},\sum_{i=1}^{3}p_{i}b_{i_{y}},\sum_{i=1}^{3}p_{i}b_{i_{z}}\right],$
(4)

where $\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}$ are the primitive vectors of the
reciprocal lattice [47]. For a total number of plane waves $N$, the integer
vectors $\bm{p}$ contained in the set

$\mathcal{G}=\left[-\frac{N^{1/3}}{2}+1,\frac{N^{1/3}}{2}-1\right]^{3},$ (5)

define a uniform grid of points in the reciprocal lattice.

Previous works [30, 48, 19, 5] have argued in favour of using first
quantization techniques to accommodate the large number of plane waves that
are required for accurate simulations. We follow that strategy in this work.
In first quantization, the plane wave representation of the operators $T$, $U$
and $V$ are given by [19, 5]:

$\displaystyle T$
$\displaystyle=\sum_{i=1}^{\eta}\sum_{p\in\mathcal{G}}\frac{\|\bm{G}_{p}\|^{2}}{2}\ket{\bm{p}}\bra{\bm{p}}_{i},$
(6) $\displaystyle U$
$\displaystyle=-\frac{4\pi}{\Omega}\sum_{i=1}^{\eta}\sum_{q\in\mathcal{G}}\sum_{\begin{subarray}{c}\nu\in\mathcal{G}_{0}\\\
(\bm{q}-\bm{\nu})\in\mathcal{G}\end{subarray}}\frac{\sum_{I=1}^{L}Z_{I}e^{i\bm{G}_{\nu}\cdot\bm{R}_{I}}}{\|\bm{G}_{\nu}\|^{2}}\ket{\bm{q-\nu}}\bra{\bm{q}}_{i},$
(7) $\displaystyle V$ $\displaystyle=\frac{2\pi}{\Omega}\sum_{i\neq
j}^{\eta}\sum_{p,q\in\mathcal{G}}\sum_{\begin{subarray}{c}\nu\in\mathcal{G}_{0}\\\
(\bm{p}+\bm{\nu})\in\mathcal{G}\\\
(\bm{q}-\bm{\nu})\in\mathcal{G}\end{subarray}}\frac{1}{\|\bm{G}_{\nu}\|^{2}}\ket{\bm{p+\nu}}\bra{\bm{p}}_{i}\ket{\bm{q}-\bm{\nu}}\bra{\bm{q}}_{j},$
(8)

where $Z_{I}$ and $\bm{R}_{I}$ are respectively the atomic number and position
of the $I$th atomic species, and $L$ denotes the number of atoms in the unit
cell. In Eqs. 7 and 8, $\bm{G}_{\nu}=\bm{G}_{q}-\bm{G}_{p}$ and
$\bm{G}_{\nu}=\bm{G}_{p}-\bm{G}_{s}=\bm{G}_{r}-\bm{G}_{q}$, respectively, and
$\mathcal{G}_{0}=\mathcal{G}\setminus(0,0,0)$. The qubit representation of the
$N$ plane waves uses $n_{p}=\lceil\log(N^{1/3}+1)\rceil$ qubits to encode each
component of the plane wave vector. Thus a total of $3\eta n_{p}$ qubits are
required for the system register.

###  II.2 Ionic pseudopotentials

Performing accurate all-electron simulations of supercell structural models of
battery materials is hampered by the huge number of plane waves that are
needed to represent the core states and the valence states near the nucleus as
sketched in Fig. 1. Different strategies such as the augmented and
orthogonalized plane wave methods [49, 50, 51] have been proposed to overcome
this limitation. However, a key step to retain the advantages of plane waves
for materials simulations was taken by Phillips, Kleinman and Antoncik (PKA)
[52, 53]. Crucially, it follows from the PKA transformation that the nuclear
potential and the core electrons can be replaced by an effective potential
known as a pseudopotential that produces the same energies of the valence
states. Furthermore, the associated pseudo wave functions superimpose the true
valence wave functions outside the core region (see Fig. 1), and can be
accurately represented using a significantly smaller number of plane waves.

Figure 1: Representation of the oscillating character of the all-electron wave
function, $\phi$, versus the smooth behaviour of the pseudo wave function
,$\phi^{\text{PP}}$, in the core region $r<r_{c}$. The pseudo wave function is
a solution of the pseudopotential, $u^{\text{PP}}$, and superimposes the
valence wave function outside the core region.

Here, we focus on the use of _ionic pseudopotentials_ (PPs) which have proven
to be highly accurate for simulating materials [38]. They are typically
generated from all-electron atomic calculations performed using density
functional theory. Due to the spherical symmetry of the Coulomb potential, the
pseudo wave functions, $\phi_{lm}(\bm{r})=\phi_{l}(r)Y_{lm}(\phi,\theta)$, are
eigenstates of the angular momentum operator and characterized by the quantum
numbers $(l,m)$. Starting from a pseudo wave function ansatz, $\phi_{l}(r)$,
the effective model potential, $u_{l}(r)$, is found by inverting the radial
Schrödinger equation [54, 55, 56, 57]. Norm-conserving PPs [58] are obtained
from pseudo wave functions enclosing the same charge as the true wave function
in the core region $r<r_{c}$, where $r_{c}$ denotes a core radius around the
nucleus. Finally, the ionic PP is obtained by removing the screening effects
due to the valence electrons in the atom. This is referred to as “unscreening”
the pseudopotential, which makes the PP transferable to different atomic
environments [59, 60].

The general form of the pseudopotential operator is [46, 61]

$u^{\text{PP}}=u^{\text{loc}}+u^{\text{SL}},$ (9)

where $u^{\text{loc}}:=u^{\text{loc}}(r)$ is a local potential, i.e., obtained
by evaluating a simple function at point $r$. It represents a screened Coulomb
potential which joins smoothly the all-electron atomic potential at some
radius $r<r_{c}$. The second operator in Eq. 9 is defined as

$u^{\text{SL}}=\sum_{l=0}^{l_{\text{max}}}\sum_{m=-l}^{l}\Delta
u_{l}(r)\ket{lm}\bra{lm},$ (10)

where $l_{\text{max}}$ is the maximum angular momentum of the core electrons
and $\langle\bm{r}|lm\rangle=Y_{lm}(\phi,\theta)$ denotes the spherical
harmonics. The $l$-dependent potential, $\Delta
u_{l}(r)=u_{l}(r)-u^{\text{loc}}(r)$, is a short-ranged potential vanishing
beyond the core radius. Moreover, in the asymptotic limit,
$r\rightarrow\infty$, the full pseudopotential, $u^{\text{PP}}$, behaves as
$-Z_{\mathrm{ion}}/r$, where $Z_{\text{ion}}=Z-\eta_{\text{core}}$ is the
effective charge of the ionic core and $\eta_{\text{core}}$ is the number of
core electrons.

The operator $u^{\text{SL}}$ in Eq. 9 has a semi-local character since its
action on a given basis function,
$\langle\bm{r}|q\rangle=\varphi_{q}(r,\phi,\theta)$, is local in the radial
coordinates but involves an integral over the angular variables
$(\phi,\theta)$ [61]. Computing its plane wave matrix elements
$u_{pq}^{\text{SL}}$ requires evaluating the radial integral [46],

$I_{pq}=\int_{0}^{\infty}dr~{}r^{2}j_{l}(G_{p}r)\Delta u_{l}(r)j_{l}(G_{q}r),$
(11)

where $j_{l}$ denotes the spherical Bessel functions and
$G_{p}=\|\bm{G}_{p}\|$. The number of such integrals scales as
$\mathcal{O}(LN^{2})$, where $L$ is the number of atoms in the material’s unit
cell and $N$ is the total number of plane waves. Typically, $N$ can be very
large in actual simulations and computing these matrix elements becomes
computationally expensive [46].

To reduce this computational cost, Kleinman and Bylander (KB) [62] proposed
the separable pseudopotentials by expressing the radial operator $\Delta
u_{l}(r)$ in a form that is separable in the radial variables. The KB
construction was further modified by Blöch [63] to construct the non-local
(NL) pseudopotential

$\displaystyle u^{\mathrm{PP}}$ $\displaystyle=u^{\text{loc}}+u^{\text{NL}},$
(12)

with the operator $u^{\text{NL}}$ given by

$u^{\text{NL}}=\sum_{l=0}^{l_{\text{max}}}\sum_{m=-l}^{l}\frac{\ket{\xi_{lm}}\bra{\xi_{lm}}}{\bra{\xi_{lm}}\phi_{lm}\rangle}.$
(13)

In Eq. 13 $\langle\bm{r}|\phi_{lm}\rangle=\phi_{l}(r)Y_{lm}(\phi,\theta)$ is a
pseudo wave function solution of the model potential $u_{l}(r)$, and
$\langle\bm{r}|\xi_{lm}\rangle=\xi_{l}(r)Y_{lm}(\phi,\theta)$ are projectors
defined as

$\xi_{lm}(\bm{r})=\left\\{\varepsilon_{l}-\left[-\frac{1}{2}\nabla^{2}+u^{\mathrm{loc}}(r)\right]\right\\}\phi_{lm}(\bm{r}),$
(14)

where $\varepsilon_{l}$ is the all-electron reference energy associated with
the pseudo wave function $\phi_{lm}(\bm{r})$. Note from Eq. 13 that computing
the matrix elements of the separable potential $u_{pq}^{\text{NL}}$, as
opposed to the matrix elements of a semi-local operator (Eq. 11), requires
only to evaluate the product of the projection operations

$\bra{\xi_{lm}}\varphi_{q}\rangle=\int
d\bm{r}~{}\xi_{lm}^{*}(\bm{r})\varphi_{q}(\bm{r}),$ (15)

which scales linearly with the number of plane waves.

The non-local operator in Eq. 13 can be generalized to use two or more
projectors per angular momentum quantum numbers [64],

$u^{\mathrm{NL}}=\sum_{l=0}^{l_{\text{max}}}\sum_{m=-l}^{l}\left[\sum_{i,j}B_{ij}\ket{\beta_{i}}\bra{\beta_{j}}\right]_{lm},$
(16)

where the matrix $B_{ij}=\bra{\phi_{i}}\xi_{j}\rangle$ is used to define the
projectors $\ket{\beta_{i}}=\sum_{j}(B^{-1})_{ji}\ket{\xi_{j}}$.
Pseudopotential operators of this form are called _generalized separable
pseudopotentials_ , and are routinely used in classical electronic structure
calculations of materials [38]. A further generalization of the non-local
operator in Eq. 16 is obtained by relaxing the norm conservation condition.
This results in the so-called ultrasoft pseudopotentials [64], which we do not
consider here.

###  II.3 Qubitization-based quantum phase estimation

The quantum algorithm presented in this work is a qubitization-based quantum
phase estimation (QPE) algorithm for computing ground-state energies of
periodic materials. In contrast to other techniques, qubitization does not
introduce any further approximations in implementing the unitary that encodes
the eigenvalues of the Hamiltonian. This makes it particularly appealing when
the accuracy of the simulation is paramount. We briefly review the
qubitization method and refer to Refs. [65, 5, 19] for further details.

To qubitize a Hamiltonian, we first write it as a linear combinations of
unitaries (LCU),

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

where each $H_{\ell}$ is a unitary and the coefficients $\alpha_{\ell}>0$ are
referred to as (unnormalized) selection probabilities. The specific choice of
an LCU has a large impact on the cost of qubitization, especially through the
parameter $\lambda=\sum_{\ell}\alpha_{\ell}$ [66].

We then define the qubitization operator

$\displaystyle
Q=(2\ket{\bm{0}}\bra{\bm{0}}-\mathbbm{1})\text{PREP}_{H}^{\dagger}\text{SEL}_{H}\text{PREP}_{H},$
(18)

with the prepare and select unitaries given by

$\displaystyle\text{PREP}_{H}\ket{\bm{0}}\ket{\psi}$
$\displaystyle=\left(\sum_{\ell}\sqrt{\frac{\alpha_{\ell}}{\lambda}}\ket{\ell}\right)\ket{\psi},$
(19) $\displaystyle\text{SEL}_{H}$
$\displaystyle=\sum_{\ell}\ket{\ell}\bra{\ell}\otimes H_{\ell},$ (20)

Note that the reflection in $Q$ and PREP act on an auxiliary register
$\ket{\bm{0}}$. The latter prepares the so-called PREP state
$\sum_{\ell}\sqrt{\frac{\alpha_{\ell}}{\lambda}}\ket{\ell}$ with amplitudes
given by the selection probabilities. It does not alter the system register
$\ket{\psi}$, which is acted upon by SEL that applies the unitaries
$H_{\ell}$. These subroutines satisfy the block-encoding equation

$\displaystyle\bra{\bm{0}}\text{PREP}_{H}^{\dagger}\ \text{SEL}_{H}\
\text{PREP}_{H}\ket{\bm{0}}=\frac{H}{\lambda}.$ (21)

The operator $Q$ is block-diagonal, with each block $Q_{k}$ corresponding to a
two-dimensional subspace $W_{k}$ spanned by an orthonormal basis
$\text{span}\\{\ket{\bm{0},\Phi_{k}},\ket{\rho_{k}}\\}$, where
$\ket{\Phi_{k}}$ is an eigenstate of $H$ with eigenvalue $E_{k}$ and
$\ket{\rho_{k}}$ is orthogonal to $\ket{\bm{0},\Phi_{k}}$. The term
qubitization refers to these effective qubit subspaces. In its eigenbasis,
$Q_{k}$ can be written as

$\displaystyle Q_{k}$
$\displaystyle=e^{i\theta_{k}}\ket{\theta_{k}}\bra{\theta_{k}}+e^{-i\theta_{k}}\ket{-\theta_{k}}\bra{-\theta_{k}},$
(22)

where $\ket{\pm\theta_{k}}$ are the eigenstates, with
$\theta_{k}=\arccos(E_{k}/\lambda)$. Thus, by applying QPE on $Q$ with an
initial state $\ket{\bm{0}}\ket{\Psi_{0}}$, where
$\ket{\bm{0}}\ket{\Psi_{0}}=\alpha\ket{\theta_{0}}+\beta\ket{-\theta_{0}}$ for
some $\alpha,\beta$, we always recover the ground state and its energy since
$\cos(\pm\theta_{0})=E_{0}/\lambda$.

In the QPE algorithm, the unitary $Q$ is controlled on the state of auxiliary
qubits, which increases the Toffoli cost of the algorithm. To avoid this, as
shown in Ref. [67], one can use a reflection to have the inverse unitary
applied when the auxiliary qubit is in state $\ket{0}$. The only requirement
for this modification to work is for SEL to be self-inverse, which our
algorithm satisfies.

When acting on an initial state $\ket{\psi}$, the QPE algorithm outputs an
estimate of the ground-state energy $E_{0}$ with probability
$p_{0}=|\braket{\psi_{0}}{\psi}|^{2}$, where $\ket{\psi_{0}}$ is the ground
state. The QPE routine needs to be repeated $O(1/p_{0})$ times on average to
retrieve $E_{0}$ with high probability. It is thus necessary to prepare an
initial state with a sufficiently large overlap with the ground state. This is
the initial state preparation problem: a crucial and daunting challenge for
quantum algorithms. Strategies for preparing initial states have been studied
in Refs. [68, 69, 25, 70], and Ref. [5] described a method for preparing a
Hartree-Fock state for periodic materials in first quantization. While we
acknowledge the importance of developing better methods for initial state
preparation, in this work we focus on the problem of reducing the cost of QPE.

###  II.4 Quantum read-only memory (QROM)

Our main use of QROM in the quantum algorithm is to prepare arbitrary states
of few qubits. There are a variety of methods in the literature for this task;
we refer to [71] for an overview of such techniques. A QROM, or a data-lookup
oracle, is an operator $O$ that reads a register $\ket{x}$ and outputs a
corresponding bitstring $\ket{\theta_{x}}$ into an auxiliary register as:
$O\ket{x}\ket{\bm{0}}=\ket{x}\ket{\theta_{x}}$ [67]. The output bitstring
$\theta_{x}$ is precomputed and available in a data-lookup table. The non-
Clifford gate cost of QROM is exponential in the size of $x$ and polynomial in
the size of $\theta_{x}$. Parallelization can decrease the depth at the
expense of using more qubits [67, 72].

We follow [72] in our description of state preparation using QROM. Consider
the target state $\ket{\psi}=\sum_{x}a_{x}\ket{x}$, where we assume $a_{x}\geq
0$ for simplicity. For any bit-string $y$ of length $w\leq n$, where $n$ is
the total number of qubits, denote by $p_{y}$ the probability that the first
$w$ qubits of $\ket{\psi}$ are in state $\ket{y}$. Define also
$\cos(\theta_{y})=\sqrt{p_{y0}/p_{y}}$, and the QROM oracles
$O_{w}\ket{y}\ket{\bm{0}}=\ket{y}\ket{\theta_{y}}$, outputting the classically
precomputed and tabulated $\theta_{y}$ up to $b$ bits of precision for each
$w$. The QROM oracles can be used iteratively to prepare any state
$\ket{\psi}=\sum_{x}a_{x}\ket{x}$ as follows. By induction, for $1\leq w\leq
n$, do:

$\displaystyle\begin{split}&\ket{\psi_{w}}=\sum_{y\in\\{0,1\\}^{w}}\sqrt{p_{y}}\ket{y}\ket{\bm{0}}\ket{\bm{0}}\\\
&\underset{O_{w}}{\mapsto}\sum_{y\in\\{0,1\\}^{w}}\sqrt{p_{y}}\ket{y}\ket{\bm{0}}\ket{\theta_{y}}\\\
&\underset{R_{w+1}}{\mapsto}\sum_{y\in\\{0,1\\}^{w}}\sqrt{p_{y}}\ket{y}\left(\cos(\theta_{y})\ket{\bm{0}}+\sin(\theta_{y})\ket{1}\right)\ket{\bm{0}}\ket{\theta_{y}}\\\
&\underset{O^{\dagger}_{w}}{\mapsto}\sum_{y\in\\{0,1\\}^{w+1}}\sqrt{p_{y}}\ket{y}\ket{\bm{0}}\ket{\bm{0}}=\ket{\psi_{w+1}}\ket{\bm{0}}\ket{\bm{0}},\end{split}$
(23)

where $R_{w+1}$ is a one-qubit rotation on the $(w+1)$-th qubit controlled on
the state of the auxiliary system, and as a slight abuse of notation we use
$\ket{\bm{0}}$ to denote all-zero states of different number of qubits. While
we assumed non-negative amplitudes $a_{x}\geq 0$, the complex phases of
$\ket{\psi}$ can be implemented by storing
$\phi_{x}=\text{arg}[a_{x}/|a_{x}|]$ and applying it in the final iteration.
In our applications, $\phi_{x}=\pm 1$, meaning the phase is always real and
the last step is a simple $Z$ gate. For future reference throughout the text,
this entire state preparation procedure is called Algorithm 1.

The precision of the rotation angle is the main source of error, and we have
the following result for the total error, proved in Section E.1:

######  Lemma E.1.

The error in the state preparation Algorithm 1 is $2^{-b}\pi n$.

There are three different types of QROM oracles that can be used in Algorithm
1 as proposed in [72]. Two of these, called Select and SelSwapDirty, are of
interest to us (the latter is our terminology for the oracle described in [72,
Fig. 1d]). Select is mostly used when $n$ is small while $b$ is large. One
important property of SelSwapDirty is the space-depth trade-off that it
offers. With the help of dirty qubits, i.e., qubits that do not need to be
initialized to any specific state, the depth of the circuit can be lowered by
parallelizing the controlled-SWAP gates used in the Swap subroutine of
SelSwapDirty. Furthermore, the dirty qubits are returned to their initial
state, therefore any qubit not undergoing simultaneous computation in the
quantum circuit can be borrowed as a dirty qubit. We briefly overview the cost
of these routines in Table 4, and we refer the reader to Appendix E for more
details.

Operation |  Additional qubits |  Toffoli Depth |  Toffoli count $\leq\cdot+\mathcal{O}(\log\cdot)$
---|---|---|---
Select | $b+\lceil\log{N}\rceil$ | $N$ | $N$
SelSwapDirty | $b(\beta+1)+\lceil\log{N}\rceil$ | $2\lceil\frac{N}{\beta}\rceil+3\lceil\log{\beta\rceil}$ | $2\lceil\frac{N}{\beta}\rceil+3b\beta$

Table 1: ([72, Table II]) Gate and qubit cost of Select and SelSwapDirty
QROMs. The space-depth trade-off is determined by $\beta\in[1,N]$. Note that
$b\beta$ qubits of the Fig. 5b implementation are dirty, while
$b+\lceil\log(N)\rceil$ are clean qubits.

##  III Quantum algorithm

We now describe the pseudopotential-based quantum phase estimation algorithm.
We begin with the construction of the plane-wave representation of the
pseudopotential operators describing the ionic cores in the material, where we
derive closed-form expressions for the matrix elements of the local and non-
local potentials (see Appendix A). We then proceed to describe the LCU
decomposition of the Hamiltonian and the implementation of the qubitization
operator. These procedures exploit the structure of the pseudopotential to
optimize the qubitization of the modified operator $U$, as discussed in
Section III.3.

###  III.1 Plane wave matrix elements of the pseudopotential operator

We focus on the description of the separable Gaussian pseudopotentials
proposed by Hartwigsen, Goedecker and Hutter (HGH) [39]. The HGH
pseudopotentials are relatively easy to define and have proven to be
transferable and accurate. The total energies of different crystals computed
by using the HGH PPs exhibit typical deviations of the order of
$1~{}\text{meV}/\text{atom}$ with respect to all-electron calculations [38].

The HGH local potential is defined as

$u^{\mathrm{loc}}(r)=-\frac{Z_{\mathrm{ion}}}{r}\mathrm{erf(\alpha
r)}+\sum_{i=1}^{4}C_{i}(\sqrt{2}\alpha r)^{2i-2}~{}e^{-(\alpha r)^{2}},$ (24)

where $\mathrm{erf}(\alpha r)$ is the error function with
$\alpha=\frac{1}{\sqrt{2}r_{\mathrm{loc}}}$, $r_{\mathrm{loc}}$ is a local
radius parameter giving the charge distribution in the core, and the $C_{i}$
are tabulated coefficients [39]. On the other hand, the non-local part
$u^{\text{NL}}$ is given by a sum of separable terms [39]

$\displaystyle u^{\text{NL}}=$
$\displaystyle\sum_{lm}\sum_{i,j=1}^{3}|\beta_{i}^{(lm)}\rangle
B_{ij}^{(l)}\langle\beta_{j}^{(lm)}|,$ (25)

where
$\langle\bm{r}|\beta_{i}^{(lm)}\rangle=\beta_{i}^{(l)}(r)Y_{lm}(\phi,\theta)$.
The radial functions $\beta_{i}^{(l)}(r)$ are Gaussian-type projectors given
by

$\beta_{i}^{(l)}(r)=A_{l}^{i}r^{l+2i-2}\mathrm{exp}\left[-\frac{1}{2}\left(\frac{r}{r_{l}}\right)^{2}\right],$
(26)

where the $A_{l}^{i}$ is a constant defined in Section A.2 and the radii
$r_{l}$ give the range of the $l$-dependent projectors. The optimized values
of the coefficients $B_{ij}^{(l)}$ and the radii $r_{l}$ are reported in Ref.
[39].

The plane-wave matrix elements of the local and non-local components of the
HGH pseudopotentials are derived in Appendix A. For an ion located at the
coordinates $\bm{R}$, the matrix elements of the local potential are given by

$\displaystyle
u_{pq}^{\mathrm{loc}}(\bm{R})=\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}~{}e^{-(G_{\nu}r_{\text{loc}})^{2}/2}\bigg{\\{}-\frac{Z_{\mathrm{ion}}}{G_{\nu}^{2}}$
$\displaystyle+\frac{\sqrt{\pi}}{2}\Big{[}C_{1}r_{\text{loc}}^{3}+C_{2}(3r_{\text{loc}}^{3}-5r_{\text{loc}}^{5}G_{\nu}^{2})$
$\displaystyle+C_{3}(15r_{\text{loc}}^{3}-10r_{\text{loc}}^{5}G_{\nu}^{2}+r_{\text{loc}}^{7}G_{\nu}^{4})$
$\displaystyle+C_{4}(105r_{\text{loc}}^{3}-105r_{\text{loc}}^{5}G_{\nu}^{2}+21r_{\text{loc}}^{7}G_{\nu}^{4}-r_{\text{loc}}^{9}G_{\nu}^{6})\Big{]}\bigg{\\}},$
(27)

where $\bm{G}_{\nu}=\bm{G}_{q}-\bm{G}_{p}$.

The matrix elements of the non-local operator in Eq. 25 are derived in Section
A.2. Typically, electronic structure calculations are performed using one or
two projectors. For the sake of simplicity, we consider the case of one
projector per angular momentum $(i=j=1)$. In this case, the plane wave matrix
elements for the HGH non-local potential are given by

$\displaystyle u^{\mathrm{NL}}_{pq}(\bm{R})=$
$\displaystyle\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\Bigg{\\{}4B_{0}r_{0}^{3}~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{0}^{2}/2}$
$\displaystyle+\frac{16B_{1}r_{1}^{5}}{3}~{}(\bm{G}_{p}\cdot\bm{G}_{q})~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{1}^{2}/2}$
$\displaystyle+\left[\frac{32B_{2}r_{2}^{7}}{15}~{}(\bm{G}_{p}\cdot\bm{G}_{q})^{2}\right.$
$\displaystyle\left.+\frac{32B_{2}r_{2}^{7}}{45}~{}(G_{p}G_{q})^{2}\right]e^{-(G_{p}^{2}+G_{q}^{2})r_{2}^{2}/2}\Bigg{\\}},$
(28)

where the coefficient $B_{l}:=B_{11}^{(l)}$. To see roughly how the projector
expression in Eq. 25 could lead to the expression above, we show one example
of how the equation decomposes to a sum of projections. Consider the first
term in Eq. 28
$4B_{0}r_{0}^{3}~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{0}^{2}/2}\ket{\bm{p}}\bra{\bm{q}}$.
When summed over $\bm{p},\bm{q}$, this can be expressed as a scalar multiple
of the projection onto the Gaussian superposition state

$\displaystyle\sum_{\bm{p}\in\mathcal{G}}e^{-G_{p}^{2}r_{0}^{2}/2}\ket{\bm{p}}.$
(29)

A similar rewriting applies to other terms in Eq. 28, and involves projection
onto (derivatives) of Gaussian superpositions. See Appendix D for more
details.

###  III.2 The pseudopotential Hamiltonian

By including the ionic pseudopotentials, the all-electron problem defined by
the Hamiltonian in Eq. 2 transforms into a valence-only electron problem where
$\eta_{\text{val}}=\sum_{I=1}^{L}Z_{\text{ion}_{I}}$. This results in a
substantial reduction of the total number of electrons, plane waves $N$, and
the overall norm of the Hamiltonian describing the valence electrons.

In this approach, the expressions for the operators $T$ and $V$ (Eqs. 6 and 8)
remain formally identical. However, the operator $U$ accounting for the
electron-nuclei interactions in the all-electron case needs to be defined
using the effective pseudopotentials describing the ionic cores:

$\displaystyle
U=\sum_{i=1}^{\eta}\sum_{I=1}^{L}\frac{-Z_{I}}{|\bm{r}_{i}-\bm{R}_{I}|}\rightarrow
U=\sum_{i=1}^{\eta_{\text{val}}}\sum_{I=1}^{L}u^{\text{PP}}(\bm{r}_{i},\bm{R}_{I}).$
(30)

Using Eq. 12, the plane-wave representation of the operator $U$ becomes

$\displaystyle U=U_{\text{loc}}+U_{\text{NL}},$ (31) $\displaystyle
U_{\text{loc}}:=\sum_{i=1}^{\eta_{\text{val}}}\sum_{p,q=1}^{N}\sum_{I=1}^{L}u_{pq}^{\text{loc}}(\bm{R}_{I})|\bm{p}\rangle\langle\bm{q}|_{i},$
(32) $\displaystyle
U_{\text{NL}}:=\sum_{i=1}^{\eta_{\text{val}}}\sum_{p,q=1}^{N}\sum_{I=1}^{L}u_{pq}^{\text{NL}}(\bm{R}_{I})|\bm{p}\rangle\langle\bm{q}|_{i},$
(33)

where $u_{pq}^{\text{loc}}(\bm{R}_{I})$ and $u_{pq}^{\text{NL}}(\bm{R}_{I})$
are the matrix elements given by Eqs. 27 and 28, respectively. This yields the
pseudopotential Hamiltonian

$\displaystyle H=T+U_{loc}+U_{NL}+V.$ (34)

###  III.3 Linear Combination of Unitaries

We study the structure of the matrix entries of each of the four operators in
the pseudopotential Hamiltonian when deriving the LCU for $H$ as in Eq. 34.
There are two main differences from the setting in [19]: (i) We consider the
general case of non-cubic lattices, that is, the primitive vectors
$\bm{a_{i}}$ have different lengths and are not orthogonal. Hence, the
reciprocal lattice vector

$\displaystyle\bm{G}_{p}$
$\displaystyle=\sum_{\omega=1}^{3}p_{\omega}\bm{b}_{\omega}$
$\displaystyle=\frac{2\pi}{\Omega}\big{(}p_{1}(\bm{a}_{2}\times\bm{a}_{3})+p_{2}(\bm{a}_{3}\times\bm{a}_{1})+p_{3}(\bm{a}_{1}\times\bm{a}_{2})\big{)},$
(35)

can no longer be substituted by $\bm{G}_{p}=2\pi\Omega^{-1/3}\bm{p}$, and (ii)
While the local term $U_{loc}$ resembles the operator $U$ in the all-electron
setting, the non-local term $U_{NL}$ has a radically different structure.
Dealing with the complexity of this non-local term is one of the biggest
challenges we face in deriving an efficient decomposition.

We briefly show how the LCUs are derived for $T,V,$ and $U_{loc}$. The case of
$U_{NL}$ is explained at a high level, with details appearing in Appendix D.
Hereafter, we make use of the following conventions. We use
$\eta:=\eta_{\text{val}}$ to denote the number of valence electrons in the
pseudopotential Hamiltonian. To make our future discussions more precise, we
define the compound index $I=(t,t_{j})$ where $t$ indicates the atomic species
of the $I$-th ionic core, and $t_{j}$ enumerate the cores of that type
(reading for example as the ‘third oxygen atom’). To make the reference to
nuclei $I$ more explicit, we may sometimes use $t_{I}$ and $t_{I,j}$.

LCU for $T$. Recall that $T$ is a diagonal operator with entries
$\frac{G_{p}^{2}}{2}$. We rewrite

$G_{p}^{2}=\bm{G}_{p}\cdot\bm{G}_{p}=\sum_{\omega,\omega^{\prime}}\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle\bm{p}_{\omega}\bm{p}_{\omega^{\prime}},$
(36)

and take the binary expansion

$\displaystyle\bm{p}_{\omega}$
$\displaystyle=(-1)^{\bm{p}_{\omega,n_{p}-1}}\sum_{r}2^{r}\bm{p}_{\omega,r},$
$\displaystyle\bm{p}_{\omega^{\prime}}$
$\displaystyle=(-1)^{\bm{p}_{\omega^{\prime},n_{p}-1}}\sum_{s}2^{s}\bm{p}_{\omega^{\prime},s},$
(37)

where $\bm{p}_{\omega^{\prime},s}$ is the $s$-th bit of $\bm{p}_{\omega}$.
Note the signed integer representation, where the $(n_{p}-1)$-th bit
determines the sign of the coordinate. One can apply the same trick as in the
orthonormal lattice case [19], rewriting
$\bm{p}_{\omega,r}\bm{p}_{\omega^{\prime},s}=(1+(-1)^{\bm{p}_{\omega,r}\bm{p}_{\omega^{\prime},s}+1})/2$.
By doing so, we reach the following decomposition:

$\displaystyle\begin{split}T&=\sum_{j=1}^{\eta}\sum_{\bm{p}\in\mathcal{G}}\frac{G_{p}^{2}}{2}\ket{\bm{p}}_{j}\bra{\bm{p}}_{j}=\\\
&\frac{1}{4}\sum_{j=1}^{\eta}\sum_{1\leq\omega,\omega^{\prime}\leq
3}|\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle|\sum_{1\leq r,s\leq
n_{p}-2}2^{r+s}\sum_{b\in\\{0,1\\}}\\\
&\sum_{\bm{p}\in\mathcal{G}}(-1)^{\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle\bm{p}_{\omega}\bm{p}_{\omega^{\prime}})}(-1)^{b(\bm{p}_{\omega,r}\bm{p}_{\omega^{\prime},s}+1)}\ket{\bm{p}}_{j}\bra{\bm{p}}_{j}.\end{split}$
(38)

Notice the appearance of the inner product in
$(-1)^{\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle)}$,
which yields one for an orthogonal lattice. The unitaries are in the last line
of the equation. The LCU in Eq. 38 is also used for the kinetic term of the
all-electron calculations as the lattices in our case studies are no longer
orthonormal.

######  Remark III.1.

When the lattice is orthogonal, the sum corresponding to $b=0$ is a multiple
of the identity, which can be omitted in the qubitization of $H$. This
shifting is used in our case studies with orthogonal lattices, giving an
improvement over the LCU proposed in Ref. [19], as it decreases the value of
$\lambda$ and thus the simulation cost.

LCU for $V$. The unitaries used here are signed translations of the lattice by
the momentum vector $\bm{\nu}$. This strategy applies, given that the matrix
entry $V_{(\bm{p}^{\prime}_{i},\bm{q}^{\prime}_{j}),(\bm{p}_{i},\bm{q}_{j})}$
of $V$ is nonzero only when
$\bm{p}^{\prime}_{i}-\bm{p}_{i}=\bm{\nu}=\bm{q}_{j}-\bm{q}^{\prime}_{j}\in\mathcal{G}_{0}$,
and its value depends only on $\bm{\nu}$ (see [5, App. E.2]). This term
identically appears in the all-electron setting. Furthermore, in that case, as
the term $\frac{1}{G_{\nu}^{2}}$ is shared by $U$ and $V$ selection
probabilities, the LCU for $U$ is similarly derived. See Appendix F for more
details. The LCU is given by

$\displaystyle\begin{split}V&=\sum_{\bm{\nu}\in\mathcal{G}_{0}}\frac{\pi}{\Omega
G_{\nu}^{2}}\sum_{i\neq j=1}^{\eta}\sum_{b\in\\{0,1\\}}\\\
&\sum_{\bm{p},\bm{q}\in\mathcal{G}}\left((-1)^{b([\bm{p}+\bm{\nu}\notin\mathcal{G}]\vee[\bm{q}-\bm{\nu}\notin\mathcal{G}])}\ket{\bm{p}+\bm{\nu}}_{i}\bra{\bm{p}}_{i}\ket{\bm{q}-\bm{\nu}}_{j}\bra{\bm{q}}_{j}\right),\end{split}$
(39)

where as before the unitaries are denoted in the last line of the equation.

LCU for $U_{loc}$. The LCU for the local term is derived similarly to $V$,
given its symmetry with respect to any translation of the lattice by
$\bm{\nu}\in\mathcal{G}_{0}$. However, the selection probabilities are
different, as $U_{loc}$ entries include an additional exponentially decaying
term in the numerator, denoted by $\gamma_{I}(G_{\nu})$ and defined below. The
resulting LCU, with unitaries corresponding to the last line of the equation,
is given by

$\displaystyle\begin{split}&U_{loc}=\sum_{j=1}^{\eta}\sum_{\bm{\nu}\in\mathcal{G}_{0}}\sum_{I=1}^{L}\frac{2\pi|\gamma_{I}(G_{\nu})|}{\Omega
G_{\nu}^{2}}\sum_{b\in\\{0,1\\}}\\\
&\sum_{\bm{q}\in\mathcal{G}}\left((-1)^{\text{sgn}(\gamma_{I}(G_{\nu}))}e^{i\bm{G}_{\nu}\cdot\bm{R}_{I}}(-1)^{b[(\bm{q}-\bm{\nu})\notin\mathcal{G}]+1}\ket{\bm{q}-\bm{\nu}}\bra{\bm{q}}_{j}\right),\end{split}$
(40)

where

$\displaystyle\begin{split}\gamma_{I}(G_{\nu})&:=\Big{(}e^{-(G_{\nu}^{2}r_{loc}^{2})/2}\big{(}-Z_{\mathrm{ion}}+\\\
&(C_{1}+3C_{2})\frac{\sqrt{\pi}r_{loc}^{3}G_{\nu}^{2}}{2}-C_{2}\frac{\sqrt{\pi}r_{loc}^{5}G_{\nu}^{4}}{2}\big{)}\Big{)}.\end{split}$

The parameters $r_{loc},C_{1},C_{2}$ were previously introduced in Eq. 27.
They depend on the atomic type $t_{I}$, so the function $\gamma_{I}$ is
determined by $t_{I}$. We use this fact to change the notation to $\gamma_{t}$
when the context is clear.

LCU for $U_{NL}$. The main insight in deriving an LCU for the non-local term
is to exploit the projector representation of the operator as in Eq. 25. We
then break down each projection $P$ into a sum of the identity and a
reflection operator $\frac{1}{2}\mathbbm{1}-\frac{1}{2}(\mathbbm{1}-2P)$. This
leads to a smaller $\lambda$ and more efficient strategies for the
qubitization of $U_{NL}$. This decomposition also gives rise to the identity
terms that lead to a shifted Hamiltonian. For convenience, we still refer to
$U_{NL}$ and $H$ by the same name after the shifting. A full derivation can be
found in Appendix D, and the resulting LCU is given by

$\displaystyle\begin{split}&U_{NL}=\sum_{j=1}^{\eta}\sum_{I,\sigma}|c_{I,\sigma}|\\\
&(-1)^{\text{sgn}(c_{I,\sigma})}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,\sigma}}\bra{\Psi_{I,\sigma}})R(\bm{R}_{I}),\end{split}$
(42)

where

$\displaystyle
R(\bm{R}_{I})\ket{\bm{p}}=e^{i\bm{G}_{p}\cdot\bm{R}_{I}}\ket{\bm{p}}.$ (43)

As before, the second line in the equation for the LCU are the unitaries in
the decomposition. The coefficients $c_{I,\sigma}$ are complicated expressions
representing sums over Gaussian terms, which we define fully in Appendix D.
For example, the simplest one corresponds to

$\displaystyle c_{I,0}:=\frac{-8\pi
r_{0}^{3}B_{0}(\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}})}{\Omega}.$ (44)

These coefficients depend on the atomic specie $t_{I}$ and the label $\sigma$
that denotes the type of Gaussian superposition defining the states
$\ket{\Psi_{I,\sigma}}$. There are eleven different choices for
$\sigma\in\\{0,\ldots,10\\}$. All corresponding Gaussian states are derived in
Appendix D, with $\ket{\Psi_{I,\sigma}}$ for $\sigma=0$ shown in Eq. 29.
Finally, $\text{sgn}(x)=0$ if $x\geq 0$ and is equal to $1$ otherwise. Notice
we need to use this sign function as the selection probabilities must be
positive.

We make a few comments on the generalizability of the LCUs. Our preference in
Section III.3 to only include $C_{1},C_{2}$, and not the higher order terms in
Eq. 27 is motivated by the fact that $C_{3}=C_{4}=0$ for all of our case
studies. Nevertheless, the algorithm subroutines that are relevant to
$\gamma_{I}(G_{\nu})$ apply without changes to any material for which
$C_{3},C_{4}\neq 0$. In fact, the LCUs above are applicable to any
pseudopotential as defined in Sections A.1 and 90. Indeed, for the local term,
the LCU is identical to LABEL:eq:lcu_loc. For the non-local term, if one
includes more than one projector, then the LCU will involve more Gaussian
superpositions $\ket{\Psi_{I,\sigma}}$, some of which do not appear in our
explicit derivation in Appendix D. The generalization of our subroutines is
straightforward as well, especially given the fact that we employ QROM for the
preparation of superpositions, where the functional defining the amplitudes
could change according to the specific parameters of the pseudopotential.

###  III.4 Breakdown of the qubitization subroutines

The qubitization of $H=T+U_{loc}+U_{NL}+V$ is roughly broken down into a
qubitization of each term. This means defining PREP and SEL subroutines for
each of the four operators, for instance SELT for the kinetic energy term. We
use this notation from now on. Henceforth, ‘AE’ refers to the all-electron
setting, with Hamiltonian $H=T+U+V$ given by Eqs. 6, 7 and 8 and no assumption
on the lattice. In our discussion of PREP and SEL, we mention the needed
adjustments to the AE algorithm when the lattice is not orthonormal. The all-
electron setting with an orthonormal lattice, hereafter referred to by ‘OAE’,
has been qubitized in [19, 5].

####  III.4.1 Prepare operator (PREP)

Below, the target PREP state for the pseudopotential Hamiltonian is shown, as
implied by their LCUs in Eqs. 38, 39, LABEL:eq:lcu_loc and LABEL:eq:lcu_NL.
The error analysis is done in Section IV. We refer to [19, Eq. (48)] for the
similar equation in the OAE setting.

$\displaystyle\text{PREP}_{H}\ket{\bm{0}}=$ (45)
$\displaystyle\big{(}\sum_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}\in\\{0,1\\}^{2}}\sqrt{\frac{\lambda_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}}{\lambda}}\ket{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\big{)}_{\mathcal{X}}$
(46)
$\displaystyle\otimes\big{(}\sum_{\omega,\omega^{\prime}\in\\{1,2,3\\}}\frac{\sqrt{|\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle|}}{\|\sum\limits_{\omega,\omega^{\prime}\in\\{1,2,3\\}}|\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle|\|}\ket{\omega,\omega^{\prime},\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle)}\big{)}_{f}$
(47)
$\displaystyle\otimes\left(\frac{1}{2^{n_{p}-1}-1}\sum_{r,s=0}^{n_{p}-2}2^{(r+s)/2}\ket{r}_{g}\ket{s}_{h}\right)$
(48)
$\displaystyle\otimes\ket{+}_{b}\otimes\frac{1}{\sqrt{\eta}}\left(\sqrt{\eta-1}\ket{0}_{c}\sum_{i\neq
j=1}^{\eta}\ket{i}_{d}\ket{j}_{e}+\ket{1}_{c}\sum_{j=1}^{\eta}\ket{j}_{d}\ket{j}_{e}\right)$
(49)
$\displaystyle\otimes\left(\sqrt{\frac{P_{\nu,V}}{\lambda_{\nu,V}}}\ket{0}_{j_{V}}\sum_{\bm{\nu}\in\mathcal{G}_{0}}\frac{1}{G_{\nu}}\ket{\bm{\nu}}_{k_{V}}+\sqrt{1-P_{\nu,V}}\ket{1}_{j_{V}}\ket{\bm{\nu}^{\perp}}_{k_{V}}\right)$
(50)
$\displaystyle\otimes\Big{(}\frac{1}{\sqrt{\lambda_{\nu,loc}}}\sum_{I,\bm{\nu}\in\mathcal{G}_{0}}\frac{|\gamma_{I}(G_{\nu})|^{1/2}}{G_{\nu}}\ket{\bm{\nu}}_{k_{loc}}\ket{I}_{k^{\prime}_{loc}}\ket{\text{sgn}(\gamma_{I}(G_{\nu}))}_{s_{loc}}\Big{)}$
(51)
$\displaystyle\otimes\frac{1}{\sqrt{\sum|c_{I,\sigma}|^{2}}}\left(\sum_{I,\sigma}\sqrt{|c_{I,\sigma}|}\ket{I}_{k^{\prime}_{NL}}\ket{\sigma}_{k_{NL}}\ket{\text{sgn}(c_{I,\sigma})}_{s_{NL}}\right).$
(52)

Each register is denoted by a subscript, such as $\mathcal{X}$ which labels
the first register. While we borrow techniques from [19] to prepare some
registers, there are also adjustments and new registers:

  * •

We have to qubitize four operators instead of three.

  * •

As the lattice is no longer orthonormal, the selection probabilities are more
general, and in cases like $T$, they include the inner product of the
reciprocal lattice vectors.

  * •

When using pseudopotentials, the selection probabilities of $U_{loc}$ and $V$
are not proportional to each other, hence their corresponding momentum state
superposition cannot be shared. This is in contrast to the AE case, where the
momentum state superposition of $U$ and $V$ is shared.

  * •

The preparation of the state in (51) for the local term is via QROM, instead
of the usual inequality test involving quantum arithmetic [19, Sec. II.C.].
The same holds for the register storing the amplitudes $\sqrt{|c_{I,\sigma}|}$
for $U_{NL}$.

We give a sketch of the preparation of the state of each register, with more
details provided in Appendix H.

  1. 1.

The state of register $\mathcal{X}$ is made of two qubits, and is a
superposition prepared by QROM with amplitudes
$\sqrt{\lambda_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}/\lambda}$,
where
$\lambda_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$
is the sum of the selection probabilities of the LCU for $\textstyle\chi$. We
enumerate $T,V$ with
$\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}=00,01$
and $U_{loc},U_{NL}$ with
$\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}=10,11$.
This state enables us to end up with the desired
$\frac{H}{\lambda}=(T+U_{loc}+U_{NL}+V)/\lambda$ instead of
$\frac{T}{\lambda_{T}}+\frac{U_{loc}}{\lambda_{loc}}+\frac{U_{NL}}{\lambda_{NL}}+\frac{V}{\lambda_{V}}$.

  2. 2.

The state $f$ is made of five qubits, two for each of
$\omega,\omega^{\prime}$, and one storing
$\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle)$. This is
prepared using QROM and is a part of the PREP state of $T$. There is a failure
probability in the preparation of this state, where ineligible states are
flagged by an additional qubit. This is not shown above to avoid cluttering.

  3. 3.

The superposition in registers $g,h$ is the same one from OAE. It is prepared
by implementing a unary state using controlled-Hadamards and bit flips [5,
Fig. II.]. These two registers along with $f,b,e$ form the PREP state of $T$.

  4. 4.

The state of register $b$ is given by the $\ket{+}$ eigenstate of Pauli $X$.
It is used in the PREP state of $T,U_{loc},V$, where the unitaries in the LCUs
are signed permutations. Note that $b$ is not part of the PREP state of $T$
when the lattice is orthogonal (Remark III.1).

  5. 5.

The states of registers $c,d,e$ also appear in the OAE setting, and form the
superposition over the electrons, which is needed in first quantization as
part of the PREP state of all four operators. The technique to prepare the
uniform superpositions is available in [21, App. A.2], and is followed by
checking whether $i=j$, with the result stored in the additional qubit of
register $c$.

  6. 6.

The superposition over the momentum $\bm{\nu}$ corresponding to $V$ is given
in registers $j_{V},k_{V}$. This is prepared using an inequality test
(Appendix F) followed by amplitude amplification to increase the probability
of success, flagged by $\ket{0}_{j_{V}}$. The procedure is very similar to the
OAE setting [19, Sec. II.C.], except that we use QROM instead of quantum
arithmetic to compute one side of the inequality test.

  7. 7.

The superposition with registers $k_{loc},k^{\prime}_{loc},s_{loc}$ correspond
to the momentum state of $U_{loc}$. The register $k^{\prime}_{loc}$ enumerates
all nuclei $I=(t,t_{j})$. While for $V$ the similar superposition is created
using an inequality test, according to our simulations, the exponentially
small term $\gamma_{I}(G_{\nu})$ makes the rejection of the inequality test
happen too often, leading to excessive rounds of amplitude amplifications that
increase the Toffoli cost. Therefore QROM is used to prepare almost the entire
superposition. Note that $\gamma_{I}$ only depends on the atomic species, so
QROM prepares the superposition over register $k_{loc},s_{loc}$ and the atomic
species index $t$ in the register $k^{\prime}_{loc}$. What remains to be done
is a uniform superposition over all nuclei $t_{j}$ of each specie $t$, which
is created using the technique in [21, App. A.2].

  8. 8.

The registers $k_{NL},k^{\prime}_{NL},s_{NL}$ give a superposition for the
PREP state of $U_{NL}$. The coordinate $I$ in register $k^{\prime}_{NL}$
enumerates all nuclei, with the same splitting mentioned above for
$k^{\prime}_{loc}$. To prepare Eq. 52, recall that $c_{I,\sigma}$ only depends
on the atomic type $t_{I}$ and the Gaussian superposition type $\sigma$.
Therefore, a QROM produces the superposition over
$\ket{t}_{k^{\prime}_{NL}}\ket{\sigma}_{k_{NL}}\ket{\text{sgn}(c_{t,\sigma})}$.
Then a uniform superposition over the $t_{j}$ nuclei of type $t$ is
implemented, giving the desired Eq. 52.

####  III.4.2 Select operator (SEL)

We mostly borrow the corresponding implementation in the OAE case [19] for
every SEL operator except $\text{SEL}_{NL}$, while adjusting for general
lattices and the different registers holding the momentum state for $U_{loc}$
and $V$. We devote more explanation to $\text{SEL}_{NL}$ as it is the operator
with no similar precedent in the literature. Nevertheless, this section is not
a detailed compilation, especially for $\text{SEL}_{NL}$ which is the most
involved; we refer to Appendix H for more details.

There are some commonalities among all SEL operators, which we briefly
discuss. The action of each operator is controlled on a register that flags
the success of the corresponding state preparation. For example for $T$, we
need to check three conditions:

  * •

The state
$\ket{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}_{\mathcal{X}}$
is equal to $\ket{00}_{\mathcal{X}}$,

  * •

The register $c$ flags the success of $i=j$ in registers $d,e$ (i.e.,
$\ket{1}_{c}$),

  * •

The ancilla attached to register $f$ flags the meaningful basis states
($\omega\neq 4$ and $\omega^{\prime}\neq 4$) in the superposition prepared by
QROM.

If all the above conditions hold, then we get $\ket{0}_{T}$, and get
$\ket{1}_{T}$ otherwise, in which case SELT acts as the identity. Checking
correctness of the PREP states can be performed using a few logical gates
(Toffoli, CNOT, X). Note that these checks are performed as part of the PREP
procedure, but we find it more informative to introduce them here.

The second design shared by all select operators is a common CSWAP circuit and
its inverse. This circuit first copies $\ket{\bm{p}}_{i}$ (controlled on
$\ket{i}_{d}$) into an auxiliary register, swaps it back into its place after
the relevant SEL operations are carried out, then does the same for
$\ket{\bm{p}}_{j}$ controlled on $\ket{j}_{e}$. This ensures that the SEL
operations are all done on a single auxiliary register, obviating the need for
the far more costly controlled operations directly on the system register [5,
Eq. (E27)]. In the AE case, this is actually the costliest part of
implementing the select operator. We now discuss how to implement each select
operator, the sum of which is the desired $\text{SEL}_{H}$.

1\. $\text{SEL}_{T}$: The transformation is given by

$\displaystyle\begin{split}\text{SEL}_{T}:&\ket{0}_{T}\ket{+}_{b}\ket{j}_{e}\ket{\omega,\omega^{\prime},\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle)}_{f}\ket{r}_{g}\ket{s}_{h}\ket{\bm{p}}_{j}\to\\\
&(-1)^{b(\bm{p}_{\omega,r}\bm{p}_{\omega^{\prime},s}+1)+\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle\bm{p}_{\omega}\bm{p}_{\omega^{\prime}})}\ket{0}_{T}\ket{+}_{b}\ket{j}_{e}\\\
&\otimes\ket{\omega,\omega^{\prime},\text{sgn}(\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle)}_{f}\ket{r}_{g}\ket{s}_{h}\ket{\bm{p}}_{j}.\end{split}$
(53)

The action above is essentially a phase, which is a combination of multi-
controlled Z gates. The sign of the inner products appearing in the phase is
the term making the distinction from the OAE case ([19, Eq. (49)]).

2\. $\text{SEL}_{V}$: The implementation follows the OAE setting [19, Eq.
(51)]:

$\displaystyle\begin{split}\text{SEL}_{V}:&\ket{0}_{V}\ket{b}_{b}\ket{i}_{d}\ket{j}_{e}\ket{0}_{c}\ket{0}_{j_{V}}\ket{\bm{\nu}}_{k_{V}}\ket{\bm{p}}_{i}\ket{\bm{q}}_{j}\to\\\
&(-1)^{b([(\bm{p}+\bm{\nu})\not\in\mathcal{G}]\lor[(\bm{q}-\bm{\nu})\not\in\mathcal{G}])}\ket{0}_{V}\ket{b}_{b}\ket{i}_{d}\ket{j}_{e}\ket{0}_{c}\\\
&\otimes\ket{0}_{j_{V}}\ket{\bm{\nu}}_{k_{V}}\ket{\bm{p}+\bm{\nu}}_{i}\ket{\bm{q}-\bm{\nu}}_{j}.\end{split}$
(54)

The addition and subtraction along with the phase implementation are all
controlled on $\ket{0}_{V}$. The non-diagonal superposition
$\ket{i}_{d}\ket{j}_{e}$ over the electrons, flagged by $\ket{0}_{c}$, is
acted upon by $\text{SEL}_{V}$, while the rest of the SEL operators only act
on the $e$ register.

3\. $\text{SEL}_{loc}$: The operator acts in two main steps:

$\displaystyle\begin{split}&\ket{0}_{loc}\ket{b}_{b}\ket{j}_{e}\ket{\bm{\nu}}_{k_{loc}}\ket{I}_{k^{\prime}_{loc}}\ket{\text{sgn}(\gamma_{I}(G_{\nu}))}_{s_{loc}}\ket{\bm{0}}_{\bm{R}}\ket{\bm{p}}_{j}\to\\\
&(-1)^{b[(\bm{p}-\bm{\nu})\not\in\mathcal{G}]}\ket{0}_{loc}\ket{b}_{b}\ket{j}_{e}\ket{\bm{\nu}}_{k_{loc}}\ket{I}_{k^{\prime}_{loc}}\ket{\text{sgn}(\gamma_{I}(G_{\nu}))}_{s_{loc}}\\\
&\otimes\ket{\bm{R}_{I}}_{\bm{R}}\ket{\bm{p}-\bm{\nu}}_{j}.\end{split}$ (55)

This first step illustrates a controlled sign, along with a controlled
subtraction $\ket{\bm{p}}_{j}\rightarrow\ket{\bm{p}-\bm{\nu}}_{j}$. It also
shows a QROM that reads $\ket{x}_{loc}\ket{I}_{k^{\prime}_{loc}}$ and outputs
$\bm{R}_{I}$ into the register $\ket{\cdot}_{\bm{R}}$ only if $x=0$. This is
then used to implement the phase $e^{i\bm{G}_{\nu}\cdot\bm{R}_{I}}$, following
the techniques in [19, Eq. (50)]. The second step maps the previous state to

$\displaystyle\begin{split}&\to
e^{i\bm{G}_{\nu}\cdot\bm{R}_{I}}(-1)^{b[(\bm{p}-\bm{\nu})\not\in\mathcal{G}]+\text{sgn}(\gamma_{I}(G_{\nu}))}\ket{0}_{loc}\ket{b}_{b}\\\
&\otimes\ket{j}_{e}\ket{\bm{\nu}}_{k_{loc}}\ket{I}_{k^{\prime}_{loc}}\ket{\text{sgn}(\gamma_{I}(G_{\nu}))}_{s_{loc}}\ket{0}_{\bm{R}}\ket{\bm{p}-\bm{\nu}}_{j}.\end{split}$
(56)

Here we erase the register $\bm{R}$ by taking the inverse of the QROM, and
finish by applying $\text{Z}_{s_{loc}}$ controlled on $\ket{0}_{loc}$. As a
remark, the AE circuit for $\text{SEL}_{U}$ is exactly the same as its OAE
implementation in [19].

4\. $\text{SEL}_{NL}$: The unitary from the LCU in LABEL:eq:lcu_NL is

$\displaystyle(-1)^{\text{sgn}(c_{I,\sigma})}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,\sigma}}\bra{\Psi_{I,\sigma}})R(\bm{R}_{I}),$
(57)

which can be broken down as a series of transformations, which we describe
below. To avoid cluttering, we present only the registers involved in each
stage. The first transformation applies a phase as follows:

$\displaystyle\begin{split}&\ket{0}_{NL}\ket{\bm{q}}_{j}\ket{I}_{k^{\prime}_{NL}}\ket{\bm{0}}_{\bm{R}}\to\\\
&\ket{0}_{NL}\ket{\bm{q}}_{j}\ket{I}_{k^{\prime}_{NL}}\ket{\bm{R}_{I}}_{\bm{R}}\to\\\
&e^{i\bm{G}_{q}\cdot\bm{R}_{I}}\ket{0}_{NL}\ket{\bm{q}}_{j}\ket{I}_{k^{\prime}_{NL}}\ket{\bm{R}_{I}}_{\bm{R}}.\end{split}$
(58)

To implement this transformation, first a QROM reads
$\ket{x}_{NL}\ket{I}_{k^{\prime}_{NL}}$ and outputs $\bm{R}_{I}$ into
$\ket{\cdot}_{\bm{R}}$ only if $x=0$. Then the phase
$e^{i\bm{G}_{q}\cdot\bm{R}_{I}}$ is applied, similar to how
$e^{i\bm{G}_{\nu}\cdot\bm{R}_{I}}$ was applied for $U_{loc}$. The next stage
is to apply the reflection onto a Gaussian state $\ket{\Psi_{I,\sigma}}$. To
do so, we need to prepare the state and apply a reflection:

$\displaystyle
U_{I,\sigma}(\mathbbm{1}-2\ket{0}_{NL}\ket{\bm{0}}\bra{\bm{0}}\bra{0}_{NL})U_{I,\sigma}^{\dagger},$
(59)

where $U_{I,\sigma}\ket{\bm{0}}=\ket{\Psi_{I,\sigma}}$ acts on the same
register as $\ket{\bm{q}}_{j}$. Note that the reflection also includes the
flag qubit $\ket{0}_{NL}$, ensuring that $\text{SEL}_{NL}$ acts by identity if
the basis state has not been successfully prepared for the non-local term.

This reflection is the most expensive part of $\text{SEL}_{H}$, so it is
worthwhile to discuss strategies to reduce its cost. Many materials have
either orthogonal or partially orthogonal lattices, i.e., when a lattice
vector is orthogonal to the other two. This crystallographic feature is
prevalent among many materials of interest, including the ones utilized as
cathode materials. For instance, about half of the crystal structures
available in the Materials Project database [73] have orthogonal or partially
orthogonal lattices. Assuming this, the Gaussian state can always be
decomposed into the tensor product of three one-dimensional (1D) Gaussian
states, or a 1D+2D Gaussian state, respectively. Since QROM cost rises
exponentially with the number of read qubits, it is important to exploit this
decomposition. As a result, for orthogonal (and partially orthogonal
lattices), one has three (two) QROMs acting in parallel and all reading
$n_{p}$ ($n_{p}$ and $2n_{p}$) qubits, instead of one QROM reading $3n_{p}$
qubits. As an example of the decomposition, we have the following for Eq. 62
when the lattice is orthogonal:

$\displaystyle\ket{\Psi_{I,10}}=$
$\displaystyle\sum_{\bm{p}_{1}}e^{-G_{p,1}^{2}r_{2}^{2}/2}\ket{\bm{p}_{1}}\otimes\sum_{\bm{p}_{2}}G_{p,2}e^{-G_{p,2}^{2}r_{2}^{2}/2}\ket{\bm{p}_{2}}$
$\displaystyle\otimes\sum_{\bm{p}_{3}}G_{p,3}e^{-G_{p,3}^{2}r_{2}^{2}/2}\ket{\bm{p}_{3}}.$
(60)

where $G_{p,\omega}$ is the $\omega$ coordinate of $\bm{G}_{p}$. Once the
reflection in Eq. 59 is implemented on $\ket{\bm{q}}_{j}\ket{0}_{j}$, it leads
to a superposition of the form
$\sum_{\bm{p}\in\mathcal{G}}\phi_{I,\sigma}(\bm{p})\ket{\bm{p}}_{j}$ for some
amplitudes $\phi_{I,\sigma}(\bm{p})$.

The final steps are (i) the application of the phase
$e^{-i\bm{G}_{p}\cdot\bm{R}_{I}}$, which is done similar to its inverse at the
beginning, (ii) the erasure of $\bm{R}_{I}$ from the register $\bm{R}$ by the
inverse of the QROM that created it, and (iii) the Z gate $\text{Z}_{s_{NL}}$
on $s_{NL}$ controlled on $\ket{0}_{NL}$. Focusing on the notable registers,
the final result is:

$\displaystyle(-1)^{\text{sgn}(c_{I,\sigma})}\ket{0}_{NL}\ket{\sigma}_{k_{NL}}\ket{I}_{k^{\prime}_{NL}}\ket{\text{sgn}(c_{I,\sigma})}_{s_{NL}}\otimes$
$\displaystyle\Big{(}\sum_{\bm{p}\in\mathcal{G}}e^{i(\bm{G}_{q}-\bm{G}_{p})\cdot\bm{R}_{I}}\phi_{I,\sigma}(\bm{p})\ket{\bm{p}}_{j}\Big{)}.$
(61)

Finally, we explain our choice of QROM for $U_{I,\sigma}$ in further detail.
First, while the preparation of a discrete Gaussian state such as
$\ket{\Psi_{I,0}}=\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}/2}\ket{\bm{p}}$ has been
specifically treated in earlier [74], we have to also prepare higher-order
derivatives of such states, for example:

$\displaystyle\ket{\Psi_{I,10}}=\sum_{\bm{p}}(G_{p,2}G_{p,3})e^{-G_{p}^{2}r_{2}^{2}/2}\ket{\bm{p}},$
(62)

Even for a diagonal covariance matrix corresponding to an orthogonal lattice,
the preparation method in [74] is inefficient compared to QROM as it assumes
an arithmetic oracle.

Furthermore, the specific case of a discrete Gaussian state with a non-
diagonal covariance matrix has not been properly investigated and optimized.
While the work in [74] provides some ideas like a simple basis change, the
details regarding non-orthogonal lattices are far more complicated and our
estimates show that the algorithm does not yield the actual state in a cost-
efficient way. In our range of applications, more recent methods like
inequality test coupled with quantum arithmetics [19] fail at providing a good
balance of the product of number of qubits and number of gates. A more
sophisticated preparation method called state preparation without coherent
arithmetic [71] is promising. However, it is also more complicated for cost
and error analysis and crucially provides less parallelization and qubit/gate
trade-off opportunities, which we frequently exploit to reduce overall cost of
the algorithm.

##  IV Error analysis and the effective value of $\lambda$

The quantum phase estimation algorithm targets a maximum total error that we
denote by $\varepsilon$. To achieve this, we need to identify all individual
sources of error in the algorithm. We use $\varepsilon_{X}$ to denote each
source of error, where $X$ will be replaced by a label describing the type of
error. These errors are ultimately related to finite precision operations
respectively using $n_{X}$ bits. The choice of $n_{X}$ further determines the
number of qubits and non-Clifford gates used in the algorithm, and is involved
in identifying the effective value of the normalization factor $\lambda$ after
qubitization. Below, we review the different sources of error, and discuss how
to compute $\lambda$. We relegate the detailed derivations to Appendices I and
J.

###  IV.1 Overview of errors and finite size approximations

Many of our errors are related to the precision of the rotation angles
$\theta$ in Algorithm 1 when building a superposition using QROM. Below, we
show the complete list of all qubitization errors:

  1. 1.

$\varepsilon_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$
is the error due to using
$n_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$
bits for the precision of the rotation angles necessary to build the
superposition of register $\mathcal{X}$ (Eq. 46) using QROM.

  2. 2.

$\varepsilon_{B}$ is a similar error, due to using $n_{B}$ bits for preparing
the superposition of register $f$ in Eq. 47.

  3. 3.

$\varepsilon_{NL}$ is the error due to using $n_{NL}$ bits in the QROM for
building the PREP state for $U_{NL}$ in Eq. 52.

  4. 4.

$\varepsilon_{M_{V}}$ is the error due to using $n_{M_{V}}$ bits in the QROM
computing the inequality test for preparing the PREPV state in Eq. 50.

  5. 5.

$\varepsilon_{M_{loc}}$ is the error due to using $n_{M_{loc}}$ bits for
preparing the local PREP state Eq. 51 with QROM.

  6. 6.

$\varepsilon_{\Psi}$ is the error due to using $n_{\Psi}$ bits for building
the superpositions $\ket{\Psi_{I,\sigma}}$ using QROM.

  7. 7.

$\varepsilon_{R}\leq\varepsilon_{R,loc}+\varepsilon_{R,NL}$ is the error due
to the finite size register $n_{R}$ used to represent $\bm{R}_{I}$ in register
$\bm{R}$ for the implementation of $\text{SEL}_{loc}$ and $\text{SEL}_{NL}$.

Finally, let $\varepsilon_{\text{QPE}}$ be the error from quantum phase
estimation. To achieve an approximation $\varepsilon$ of the ground state
energy, it is necessary and sufficient that

$\displaystyle\varepsilon^{2}\geq\varepsilon_{\text{QPE}}^{2}+(\varepsilon_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}+\varepsilon_{B}+\varepsilon_{NL}+\varepsilon_{R}+\varepsilon_{M_{V}}+\varepsilon_{M_{loc}}+\varepsilon_{\Psi})^{2}.$
(63)

The proof follows a similar argument to the OAE setting [19, Thm. 4].

Each $\varepsilon_{X}$ is upper bounded by an expression depending on $n_{X}$,
which thus determines the value of $n_{X}$ needed to obtain an error
$\varepsilon_{X}$. We estimate all errors in Appendix J, with results
summarized in Table 6. As an example, following the Lemma in Section II.4, it
can be shown that

$\displaystyle\varepsilon_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\leq\frac{4\pi}{2^{n_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}}}\lambda,$
(64)

which implies that

$\displaystyle
n_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}=\lceil\log(\frac{4\pi\lambda}{\varepsilon_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}})\rceil,$
(65)

gives an error less than or equal to
$\varepsilon_{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$.

###  IV.2 The effective value of $\lambda$

The value of $\lambda$ determines the scaling in the block-encoding of $H$ and
is a significant factor in the cost that is highly dependent on the chosen LCU
and compilation strategy. Therefore, its accurate computation is important.
The algorithm targets a PREP state that is different from the theoretical one
implied by the LCUs, and the resulting effective $\lambda$ is closely related
but technically different from the one implied by the LCUs. Note that this is
also a feature of previous work ([19, Thm. 4]).

Computing the effective $\lambda$ requires finding the effective one for each
of the four operators $T,U_{loc},U_{NL},V$. To do so, we identify any
failure/success probability embedded in the PREP implementation. Here, failure
refers to a basis state itself being inadmissible. This appears when preparing
a uniform superposition over a basis that is not a power of two, or when an
inequality test is rejected. Even upon success, the obtained amplitudes are
usually different from the desired ones. For example, we use an inequality
test to prepare amplitudes $G_{\nu}^{-1}$ for Eq. 50, and after success, we
get some complicated expression approximating $G_{\nu}^{-1}$ (see Eq. 166).

Overall, the $\lambda$ for each operator will roughly look like a sum of
effective selection probabilities $\sum\tilde{\alpha_{\ell}}$ divided by a
product of success probabilities $\prod_{i}P_{i}$. Then the effective total
value $\lambda=\lambda_{T}+\lambda_{loc}+\lambda_{NL}+\lambda_{V}$ can be
derived. We do this in Appendix I and summarize the values in Table 5. As an
example, we show $\lambda_{T}$ below. This is computed after the compilation
of the LCU of $T$ in Eq. 38. In the numerator of the expression below, we have
the explicit value for $\sum\tilde{\alpha_{\ell}}$. In the denominator, we
have a success probability related to the preparation of the uniform
superposition over pairs of electrons in Eq. 49:

$\displaystyle\lambda_{T}=\frac{\eta
2^{2n_{p}-3}\sum_{\omega,\omega^{\prime}\in\\{1,2,3\\}}|\langle\bm{b}_{\omega},\bm{b}_{\omega^{\prime}}\rangle|}{\text{Ps}(\eta,b_{r})^{2}}.$
(66)

##  V Gate and qubit cost

Having described the main steps of the algorithm and identified the sources of
error, we now quantify the number of gates and qubits needed to run the full
procedure. We follow standard practice in the literature and focus on Toffoli
gates since they constitute the bulk of non-Clifford gates used in the
algorithm.

The compilation of the algorithm involves numerous strategies, each requiring
its own separate cost estimate. This leads to an extensive analysis that
cannot be summarized in simple expressions. We have thus gathered all Toffoli
costing calculations in the Appendix Table 7 and demonstrate the derivation of
each in Appendix K. These formulas are implemented in the code that is used to
perform resource estimation calculations, publicly available at
https://github.com/XanaduAI/pseudopotentials. Here we focus on highlighting
the most expensive steps of the algorithm.

First, we recall the rough estimate for the cost of a qubitization-based QPE
algorithm:

$\displaystyle\left\lceil\frac{\pi\lambda}{\varepsilon_{\text{QPE}}}\right\rceil(2\text{PREP}_{cost}+\text{SEL}_{cost}),$
(67)

where we explicitly include the qubitization costs of PREP and SEL. The factor
of two appears because in the qubitization operator we apply PREP and its
complex conjugate. The multiplicative factor $\lambda/\varepsilon$ is the
largest contributor to the total cost. For example, for a system with
$N=10^{5}$ plane waves and a hundred electrons (Table 9), this fraction gives
a factor of about $10^{9}$, while the qubitization cost contributes roughly a
factor of $10^{5}$.

The two most costly subroutines in the qubitization part of the algorithm are:

  * •

Preparing
$\sum_{I,\bm{\nu}}\frac{|\gamma_{I}(G_{\nu})|^{1/2}}{G_{\nu}}\ket{\bm{\nu},I,\text{sgn}(\gamma_{I}(G_{\nu}))}$,
which is part of the PREP of the local term.

  * •

Performing the reflection on the Gaussian superpositions
$\ket{\Psi_{I,\sigma}}$, which is part of SEL of the non-local term.

The expressions for the Toffoli cost of these steps are roughly given by

$\displaystyle
2\big{(}2\left\lceil\frac{2^{3n_{p}+\tau+1}}{\beta_{loc}}\right\rceil+3\beta_{loc}n_{M_{loc}}(3n_{p})\big{)},$
(68)
$\displaystyle\begin{split}&6\big{(}8\left\lceil\frac{2^{n_{\text{QROM}}+1}-2^{n_{\text{QROM}}-n_{p}}}{\beta_{\Psi}}\right\rceil+3\beta_{\Psi}n_{\Psi}n_{p}\big{)}\\\
&+30\big{(}2\left\lceil\frac{2^{n_{\text{QROM}}+1}-2^{n_{\text{QROM}}-n_{p}}}{\beta_{\Psi}^{\prime}}\right\rceil+3\beta_{\Psi}^{\prime}n_{\Psi}n_{p}\big{)},\end{split}$
(69)

respectively, where $n_{\text{QROM}}=n_{p}+\lceil\log(\mathcal{T})\rceil+2$,
where $\mathcal{T}$ is the number of atomic species in the cell. Given how the
parameters $\beta_{loc},\beta_{\Psi},\beta_{\Psi}^{\prime}$ divide the leading
terms, their values are quite important in lowering the cost. However, these
scalars have to satisfy constraints based, for example, on the number of
available dirty qubits. The space-time trade-offs of QROM during the resource
estimation is explained further in Appendix M.

In our calculations, each of these constitutes the largest share of
$\text{PREP}_{cost}$ and $\text{SEL}_{cost}$, respectively. While
$\text{SEL}_{cost}$ is slightly cost-dominant for a smaller number of plane
waves, the trend is reversed as $N$ grows beyond $10^{5}$. Consequently, one
significant obstacle in further optimizing the qubitization algorithm is the
preparation of Eq. 51, where $\gamma_{I}(G_{\nu})$ is the arithmetically
involved term in Section III.3.

The derivation of the qubit cost of the algorithm is of similar complexity to
the gate cost. We provide a full analysis in Appendix L, where we reuse any
uncomputed and clean qubits whenever possible. We also make the distinction
between clean and dirty qubits as the latter is implicated in QROMs. Although
our counting is different from [19, 5], our results in Section VII follow the
same behaviour: the overwhelming contribution to the qubit cost is the size of
the system register, which requires $3\eta n_{p}$ qubits. As an example, the
total clean qubit cost for simulating a system of 408 electrons with
$N=10^{5}$ plane waves is about 9,892, while the system register uses $3\eta
n_{p}=3\cdot 408\cdot 6=7,344$ clean qubits (Table 12).

##  VI Application: lithium-excess cathode materials

In this section, we focus on the quantum simulation of lithium-excess (Li-
excess) cathode materials, which have been recently proposed for higher-
capacity cathodes [40, 75, 76]. By replacing a fraction of the transition
metals with Li atoms, Li-excess materials can potentially offer up to twice
the specific capacity, i.e., the total amount of charge stored per unit mass,
of conventional cathodes [77]. For example,
$\text{Li}_{2}\text{Mn}\text{O}_{3}$ has a theoretical capacity of 460
$\text{mAh}~{}\text{g}^{-1}$ responsible for the voltage of $\sim 4$ V [78].
This yields a specific energy of $(460~{}\text{mAh}~{}\text{g}^{-1}\times
4~{}\text{V})=1840~{}\text{Wh}~{}\text{Kg}^{-1}$ which is well above the
specific energy ($800~{}\text{Wh}~{}\text{Kg}^{-1}$ at the cathode material
level) required to enable full driving performance and significantly reduce
the cost of electric vehicles [41]. However experiments reveal average
voltages of $\sim 3.8$ V and a much smaller capacity of about
$180~{}\text{mAh}~{}\text{g}^{-1}$ following the first charging cycle [79].
This voltage decay and the abrupt capacity loss, which is common to other Li-
excess materials, is an important obstacle in designing higher-capacity
batteries.

The voltage profile of a cathode material, i.e., the voltage $V(x)$ measured
as a function of the concentration of the Li ions $x$, provides relevant
information about the Li insertion process [80]. The typical voltage profile
of Li-excess materials is sketched in Fig. 2. A distinguishing feature of this
profile is a long plateau region following the initial sloping curve during
charge. This is indicative of the cathode material transitioning from a solid
solution phase to a region where two phases of the material coexist [80].
Moreover, the hysteresis loop reveals that the removal of Li ions is
accompanied by irreversible structural transformations causing the capacity
loss observed experimentally [40].

Figure 2: a) Sketch of a typical voltage profile for Li-excess materials. The
curves depict the charging (orange) and discharging (dark-blue) processes. In
the transition metal (TM) oxidation region, the extraction of Li ions is
charge-compensated by electrons provided by the transition metals. Further Li
extraction is thought to be possible via the oxidation of oxygen atoms (excess
oxidation). Phase transformations occur within the plateau region. b) A
typical formation energy plot with respect to the Li-ion concentration $x$ in
a Li-excess cathode. The brighter (dimmed) blue circles correspond to the more
(less) stable phases in this representation. The dotted line indicates the
convex hull consisting of the most stable phases. The orange circles
correspond to the fully lithiated and delithiated phases where ab initio
calculations are strongly supported by experimental data. Differences between
stable, unstable, and inaccurate calculations (orange dimmed circle) are in
the energy scale of millielectronvolts (meV), as depicted in the inset.

The capacity gain of Li-excess materials is most frequently attributed to the
oxidation of oxygen anions ($\text{O}^{2-}\rightarrow\text{O}^{n-},n<2$)
(anion redox) [33]. This redox mechanism and alternative mechanisms have been
proposed and discussed extensively in Ref. [40] along with their relation to
structural transforms resulting in materials degradation and performance loss.
While important questions remain open, there seems to be a consensus that the
increase of oxygen hole states and lithium vacancies destabilizes the metal-
oxide chemical bonds and leads to the formation of oxygen dimers [81, 34].
Furthermore, oxygen dimerization cooperatively favors the migration of
transition metals to Li-vacancy sites in the structure, which is the main
process driving the structural transformations in Li-excess materials [82,
83].

It follows that identifying the most stable phases of the cathode material for
compositions with excess Li is crucial for determining the dominant redox
mechanisms and, more importantly, for deriving potential strategies (e.g.,
doping, modifying the crystal structure) to retain more reversible capacity.
The relevant quantity for ranking the stability of possible phases of the
delithiated material is the formation energy computed for a given Li ion
concentration [31]. For example, in the case of
$\text{Li}_{2}\text{Mn}\text{O}_{3}$, this is given by [43]

$\displaystyle E_{f}(x)=E(\text{Li}_{x}\text{Mn}\text{O}_{3})$
$\displaystyle-\frac{x}{2}E(\text{Li}_{2}\text{Mn}\text{O}_{3})$
$\displaystyle-\left(1-\frac{x}{2}\right)E(\text{Li}_{0}\text{Mn}\text{O}_{3}),$
(70)

where $E(\text{Li}_{x}\text{Mn}\text{O}_{3})$ is the ground-state energy of
the material in a given phase, and $E(\text{Li}_{2}\text{Mn}\text{O}_{3})$ and
$E(\text{Li}_{0}\text{Mn}\text{O}_{3})$ are respectively the total energies of
the stable reference materials at the end points of the charging curve ($x=2$
and $x=0$). The difference between the formation energies of different phases
is typically on the scale of 1 meV $\approx 0.04~{}\text{mHa}$ [31, 84, 43,
34], which sets the required accuracy for computing the energies needed to
evaluate Eq. 70.

The release of molecular oxygen, $\text{O}_{2}$, upon Li extraction has been
observed in experiments and also theoretically predicted [40]. Quite recently,
McColl et al. [34] investigated the stability of the delithiated phases of the
disordered rocksalt structures of $\text{Li}_{2}\text{MnO}_{2}\text{F}$ and
found that the formation of $\text{O}_{2}$, following transition metal
migration, is the dominant redox mechanism. The thermodynamical driving force
for this reaction is the oxygen-vacancy formation energy [84, 43], defined as

$\displaystyle
E_{\text{O}_{2}}(x)=E(\text{Li}_{x}\text{MnO}_{2-\delta}\text{F})$
$\displaystyle-E(\text{Li}_{x}\text{MnO}_{2-\delta}\text{F})$
$\displaystyle+\frac{\delta}{2}E(\text{O}_{2}),$ (71)

where $\delta$ denotes the number of oxygen atoms removed per formula unit and
$E(\text{O}_{2})$ is the ground-state energy of the oxygen molecule.

Different strategies have been proposed to mitigate the capacity loss of Li-
excess materials [85, 86, 87, 88, 89, 90, 91, 44]. Most of them rely on
modifying the composition and/or the atomic structure of the material to
suppress transition metal migration in the delithiated cathode. To that end,
computing accurate site energies, i.e., the ground-state energy of the
material as the transition metal occupies different lattice sites, has proven
to be useful [44]. Furthermore, from the ground-state energy of the transition
state (TS) of the material along the TM migration path [92], it is possible to
compute the activation energy $E_{\text{TS}}-E_{0}$ to describe the kinetic
pathways and, crucially, determine the reversibility of the structural
transformations of the material.

The accuracy of the electronic structure method used to simulate Li-excess
materials is key as it impacts the entire computational methodology used to
investigate these materials [40, 33, 44, 34]. Due to the practical
computational limitations of more accurate methods, density functional theory
simulations are typically employed; however, these are not able to reliably
predict stable phases of the delithiated material [40, 84]. In summary,
relevant electronic structure calculations for Li-excess materials include
formation energies of delithiated phases, oxygen-vacancy formation energies,
site energies, and activation energies for kinetic pathways. Any of these
simulations can be reduced to a series of ground-state energy calculations,
each of which can be performed using our quantum algorithm. We now study the
resources required to implement our quantum algorithm for the ground-state
energy calculations for Li-excess materials.

Figure 3: Representation of the atomic structures of the Li-excess materials selected to perform the resource estimation of the quantum algorithm. These structural models corresponds to the supercells described in Table 2. The polyhedra depict the octahedrally coordinated transition metals. These figures were generated using the VESTA package [93]. Material |  $\text{Li}_{0.5}\text{MnO}_{3}$ |  $\text{Li}_{0.75}\text{[Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}\text{]O}_{2}$ |  $\text{Li}_{0.75}\text{MnO}_{2}\text{F}$
---|---|---|---
Crystal system / space group | Monoclinic / [C 2/m] | Hexagonal / [$\text{P}6_{3}\text{mc}$] | Cubic / [Fm3m]
Supercell size | $2\times 2\times 1$ | $2\times 3\times 2$ | $3\times 2\times 2$
Lattice vectors (Å) |  $\vec{a}_{1}=(10.02,0,0)$ $\vec{a}_{2}=(0,17.32,0)$ $\vec{a}_{3}=(-1.6949,0,4.7995)$ |  $\vec{a}_{1}=(5.7081,0,0)$ $\vec{a}_{2}=(-4.2811,7.4151,0)$ $\vec{a}_{3}=(0,0,19.6317)$ |  $\vec{a}_{1}=(12.48,0,0)$ $\vec{a}_{2}=(0,8.32,0)$ $\vec{a}_{3}=(0,0,8.32)$
Supercell volume ($\text{\AA}^{3}$) | 832.9405 | 830.9604 | 863.8955
$N_{\text{atoms}}$ | 72 (8 Li, 16 Mn, 48 O) | 90 (22 Li, 14 Mn, 6 Ni, 48 O) | 76 (12 Li, 16 Mn, 16 F, 32 O)
$\eta^{\text{PP}}~{}(\eta^{\text{AE}})$ | 408 (808) | 468 (968) | 428 (836)
$N^{\text{PP}}~{}(N^{\text{AE}})$ | 55,473 ($5.8\times 10^{8}$) | 67,767 ($8.7\times 10^{8}$) | 57,655 ($6.4\times 10^{8}$)
Table 2: Crystal lattice parameters for the materials selected to perform
resource estimation. Recall that $\eta$ denotes the number of electrons in the
supercell structural models and $N$ is the number of plane waves required to
converge the ground-state energy of the material at the level of density
functional theory. The superscripts PP and AE are used to differentiate
between the pseudopotential and all-electron cases, respectively.

##  VII Resource estimation

We now estimate the number of qubits and Toffoli gates needed to implement our
pseudopotential-based algorithm for three different Li-excess materials:
lithium manganese oxide ($\text{Li}_{2}\text{MnO}_{3}$), lithium nickel-
manganese oxide
($\text{Li}[\text{Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}]\text{O}_{2}$),
denoted as LLNMO, and lithium manganese oxyfluoride
($\text{Li}_{2}\text{MnO}_{2}\text{F}$). To that end we have built supercell
structural models of these materials as described in detail in Refs. [43, 44,
34]. The atomic models are shown in Fig. 3. The lattice parameters of the
supercells are summarized in Table 2, and the procedure to delithiate the
pristine materials is described in Appendix B.

The resource estimation results targeting chemical accuracy are shown in Table
3. All calculations are performed using our resource estimation software,
which is available at https://github.com/XanaduAI/pseudopotentials. This code
allows us to compute the resources of both the pseudopotential and the all-
electron algorithms for any material. In all cases, we see that the
pseudopotential-based algorithm has gate counts that are roughly four orders
of magnitude lower than in the all-electron case. Qubit numbers are also
considerably lower, roughly needing less than half as many qubits.

Material | Qubit cost | Toffoli cost
---|---|---
PP | AE | PP | AE
$\text{Li}_{0.5}\text{MnO}_{3}$ | 9808 | 24974 | $\mathbf{5.00\times 10^{15}}$ | 2.16$\times 10^{19}$
$\text{Li}_{0.75}\text{[Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}\text{]O}_{2}$ | 11130 | 29784 | $\mathbf{4.84\times 10^{15}}$ | 3.67$\times 10^{19}$
$\text{Li}_{0.75}\text{MnO}_{2}\text{F}$ | 10260 | 26629 | $\mathbf{3.87\times 10^{15}}$ | 1.18$\times 10^{19}$
$\text{Li}_{2}\text{FeSiO}_{4}$ | 2650 | 4859 | 6.38$\mathbf{\times 10^{13}}$ | 1.59$\times 10^{17}$
Table 3: Resource estimation using $N^{\text{PP}}$ and $N^{\text{AE}}$ plane
waves as specified in Table 2 for the pseudopotential (PP) and the all-
electron (AE) algorithm, respectively. For $\text{Li}_{2}\text{FeSiO}_{4}$, we
have $N^{\text{PP}}=19,549,N^{\text{AE}}=5.46\times 10^{7}$ and
$\eta^{\text{PP}}=100,\eta^{\text{AE}}=150$. The rest of parameters needed are
supplied in [5]. Better numbers are indicated in bold. Qubit cost corresponds
to logical qubits.

We note that the cost of the pseudopotential and all-electron algorithms are
comparable even when using the same number of planewaves. For
$\text{Li}_{0.5}\text{MnO}_{3}$ and $N=10^{5}$, the Toffoli cost for PP is
$2.68\times 10^{14}$, and for the all-electron is $1.44\times 10^{14}$. We
refer to the tables in Appendix M for comparisons for each material. The cost
competitiveness for simulations using the same number of plane waves is not
just the result of only having to use half the number of electrons with the
pseudopotential, but also a result of our tailored compilation strategy, which
involves a suitable LCU to lower the value for $\lambda$, and the use of QROM
for the state preparation method. When targeting the same energy accuracy, the
cost of the all-electron algorithm is significantly higher as the number of
electrons increases, and a very large number of plane waves are needed to
properly describe the quantum system (see Table 2).

The number of plane waves used to perform resource estimation are reported in
Table 2. These numbers were estimated by performing a convergence analysis of
the total energy of the materials using density functional theory calculations
with the Perdew-Burke-Ernzerhof (PBE) exchange–correlation functional [94]
(see Appendix C). Using pseudopotentials, we achieve convergence with roughly
$6\times 10^{4}$ plane waves. For the all-electron calculations, we need
approximately $7\times 10^{8}$ plane waves.

Calculations using the HGH pseudopotentials were performed using the Quantum
ESPRESSO package [95], and the total energy was computed for the Gamma
($\Gamma$)-point only, which is sufficient to represent the large systems
studied in this work [96]. A cutoff energy of $70~{}\text{Ry}$ $\approx$ 950
eV was set for both $\text{Li}_{0.5}\text{MnO}_{3}$ and
$\text{Li}_{0.75}\text{MnO}_{2}\text{F}$, while a cutoff of 80 Ry $\approx$
1000 eV was set for
$\text{Li}[\text{Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}]\text{O}_{2}$ (see
Appendix Fig. 4). For the all-electron case, we used the localized augmented
plane-wave plus local orbitals (LAPW+lo) method as implemented in the WIEN2k
code [97]. In this case, a larger number of plane waves are needed to achieve
convergence of the total energy. The most important parameters that have to be
considered are the muffin-tin radii ($R_{\text{MT}}$) and the plane-wave
cutoff ($K_{\text{max}}$) for the expansion of the wave function in the
interstitial. The basis set cutoff parameter is defined by the product
$R_{\text{MT}}K_{\text{max}}$, which was set to 7 for all the structures in
this work. In all cases, the convergence energy was calculated without
applying structure relaxation as this would not significantly affect the total
number of plane waves needed to achieve convergence.

Using the data in Table 2, we can measure the Toffoli and logical qubit cost
of our algorithm. These estimates are benchmarked against the all-electron
setting for the same materials. Results are also shown for the dilithium iron
silicate cathode $\text{Li}_{2}\text{FeSiO}_{4}$, studied earlier in [5]. We
only report the cost of performing one round of quantum phase estimation. A
more comprehensive analysis of total cost should include the overhead due to a
limited overlap of the initial state and the cost of initial state
preparation. We can also consider the T-gate cost of the algorithm, but these
are at least an order of magnitude lower than that of the qubitization-based
QPE [5].

We also report the depth of the circuits along with the clean portion of the
qubit cost in Appendix M illustrating further the advantage of using
pseudopotentials and QROMs in our algorithm. We compute this for the following
reasons. First, studying depth is the first step towards answering the more
complex and interesting question of algorithmic runtime. Second, compared to
other state preparation methods, QROM offers far more space-depth trade-off
flexibility, something we would like to leverage. Also, given the high costs
of the algorithm, it is reasonable to assume that capable hardware would also
offer the possibility of simultaneous Toffoli applications.

To estimate the Toffoli depth, we have to specify certain parameters such as
the maximum number $n_{\text{dirty}}$ of dirty qubits available for QROM, the
maximum allowed number $n_{\text{tof}}$ of simultaneous Toffoli applications.
More details are provided in Appendix M. Notice that due to the relatively
significant use of QROMs, the PP-based algorithm stands to benefit more than
the AE algorithm on circuit depth reduction. As an example, for
$\text{Li}_{0.75}\text{[Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}\text{]O}_{2}$,
the clean qubit cost and Toffoli depth for PP are respectively 11,130 and
9.59$\times 10^{14}$, while for AE they are 29,784 and 3.59$\times 10^{19}$.

##  VIII Conclusions

This work introduced the first example of a quantum algorithm that makes use
of ionic pseudopotentials to simulate periodic materials. Our main technical
contributions are a collection of highly-optimized compilation strategies for
qubitization-based quantum phase estimation in first quantization. They are
designed to reduce the cost of implementing the algorithm despite the
mathematical complexity of pseudopotential operators. A key ingredient is the
use of quantum read-only memories to avoid performing complicated arithmetic
operations on a quantum computer. This also helps navigate tradeoffs between
the number of qubits and number of gates in the algorithm, which can be
exploited to reduce costs. Overall, we reduce the cost of the quantum
algorithm by about four orders of magnitude compared to the all-electron
approach described in Ref. [19] when applied to simulating lithium-excess
materials for a fixed target accuracy of the ground-state energy.

Using these tailored quantum algorithms, we estimated the number of qubits and
Toffoli gates needed to simulate Li-excess cathode materials proposed in the
literature. In each of these cases, a large supercell is necessary to ensure
the quality of the simulation, which considerably increases the cost of the
all-electron approach, making the use of pseudopotentials even more important.
On the other hand, even though our algorithm is applicable to any lattice, it
benefits from the orthogonality of the lattice vectors, as this reduces the
cost of quantum read-only memory techniques for preparing relevant Gaussian
superposition states. It is desirable to study strategies where the cost of
implementing the quantum algorithm is less dependent on the specific structure
of the lattice.

The accuracy of the quantum algorithm depends on the quality of the
pseudopotentials. Therefore, in selecting the Hartwigsen-Goedecker-Hutter
pseudopotential we rely on extensive benchmarking by the density functional
theory community that have identified it as accurate and transferable [38].
Additionally, we have constrained our analysis to separable pseudopotentials
with one projector per angular momentum channel in the non-local component of
the pseudopotential operator. Possible extensions of this work might consider
more than one projector and the use of ultrasoft pseudopotentials [64] to
further reduce the number of plane-wave basis functions as observed in
classical algorithms.

Quantum computing offers the potential to perform high-accuracy simulations of
strongly-correlated systems of unprecedented size. This is an outstanding
challenge for classical methods, that either cannot offer the same accuracy
guarantees, or incur prohibitive costs. Still, to realize the potential
advantages of quantum computing, further efforts are invaluable to reduce the
cost of quantum algorithms. For example, the quality of the initial state
preparation method should be addressed, since a single round of quantum phase
estimation may need to be repeated too many times for states with poor
overlap.

Overall, we have demonstrated that the method of pseudopotentials can be
effectively applied to reduce the cost of quantum algorithms to simulate
battery materials that require large supercell structural models. This can
potentially unlock the application of quantum computing to address more
complicated processes involving doped materials as well as chemical reactions
at the electrode-electrolyte interface.

##  IX Acknowledgments

We thank Matthew Kiser, Stepan Fomichev, Soran Jahangiri, and Nathan Killoran
for their fruitful comments.

## References

  * Ho et al. [2018] A. Ho, J. R. McClean, and S. P. Ong, Joule 2, 810 (2018).
  * Rice et al. [2021] J. E. Rice, T. P. Gujarati, M. Motta, T. Y. Takeshita, E. Lee, J. A. Latone, and J. M. Garcia, The Journal of Chemical Physics 154, 134115 (2021).
  * Kim et al. [2022] I. H. Kim, Y.-H. Liu, S. Pallister, W. Pol, S. Roberts, and E. Lee, Physical Review Research 4, 023019 (2022).
  * Clinton et al. [2022] L. Clinton, T. Cubitt, B. Flynn, F. M. Gambetta, J. Klassen, A. Montanaro, S. Piddock, R. A. Santos, and E. Sheridan, arXiv:2205.15256 (2022).
  * Delgado et al. [2022] A. Delgado, P. A. M. Casares, R. dos Reis, M. S. Zini, R. Campos, N. Cruz-Hernández, A.-C. Voigt, A. Lowe, S. Jahangiri, M. A. Martin-Delgado, et al., Physical Review A 106, 032428 (2022).
  * Rubin et al. [2023] N. C. Rubin, D. W. Berry, F. D. Malone, A. F. White, T. Khattar, A. E. DePrince III, S. Sicolo, M. Kühn, M. Kaicher, J. Lee, et al., arXiv:2302.05531 (2023).
  * McArdle et al. [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Kandala et al. [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Yoshioka et al. [2020] N. Yoshioka, Y. O. Nakagawa, Y. ya Ohnishi, and W. Mizukami, arXiv (2020), eprint 2008.09492.
  * Anselmetti et al. [2021] G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, and R. M. Parrish, New Journal of Physics 23, 113010 (2021).
  * Arrazola et al. [2022] J. M. Arrazola, O. Di Matteo, N. Quesada, S. Jahangiri, A. Delgado, and N. Killoran, Quantum 6, 742 (2022).
  * Lin and Tong [2022] L. Lin and Y. Tong, PRX Quantum 3, 010318 (2022).
  * Wan et al. [2022] K. Wan, M. Berta, and E. T. Campbell, Physical Review Letters 129, 030503 (2022).
  * Wang et al. [2022] G. Wang, D. Stilck-França, R. Zhang, S. Zhu, and P. D. Johnson, arXiv:2209.06811 (2022).
  * Ding and Lin [2022] Z. Ding and L. Lin, arXiv:2211.11973 (2022).
  * Babbush et al. [2016] R. Babbush, D. W. Berry, I. D. Kivlichan, A. Y. Wei, P. J. Love, and A. Aspuru-Guzik, New Journal of Physics 18, 033032 (2016).
  * Reiher et al. [2017] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proceedings of the national academy of sciences 114, 7555 (2017).
  * Kivlichan et al. [2018] I. D. Kivlichan, J. R. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan, and R. Babbush, Physical Review Letters 120, 110501 (2018).
  * Su et al. [2021] Y. Su, D. W. Berry, N. Wiebe, N. Rubin, and R. Babbush, PRX Quantum 2, 040332 (2021).
  * Motta et al. [2021] M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, and G. K. Chan, npj Quantum Information 7, 1 (2021).
  * Lee et al. [2021] J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, PRX Quantum 2, 030305 (2021).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Physical Review Letters 118, 010501 (2017).
  * Low and Chuang [2019a] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019a).
  * Berry et al. [2019] D. W. Berry, C. Gidney, M. Motta, J. R. McClean, and R. Babbush, Quantum 3, 208 (2019).
  * von Burg et al. [2021] V. von Burg, G. H. Low, T. Häner, D. S. Steiger, M. Reiher, M. Roetteler, and M. Troyer, Physical Review Research 3, 033055 (2021).
  * Low and Wiebe [2018] G. H. Low and N. Wiebe, arXiv:1805.00675 (2018).
  * Kieferová et al. [2019] M. Kieferová, A. Scherer, and D. W. Berry, Physical Review A 99, 042314 (2019).
  * Rajput et al. [2022] A. Rajput, A. Roggero, and N. Wiebe, Quantum 6, 780 (2022).
  * Childs et al. [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
  * Babbush et al. [2018a] R. Babbush, N. Wiebe, J. R. McClean, J. McClain, H. Neven, and G. K.-L. Chan, Physical Review X 8, 011044 (2018a).
  * Urban et al. [2016] A. Urban, D.-H. Seo, and G. Ceder, npj Computational Materials 2, 1 (2016).
  * Van der Ven et al. [1998] A. Van der Ven, M. Aydinol, G. Ceder, G. Kresse, and J. Hafner, Physical Review B 58, 2975 (1998).
  * Seo et al. [2016] D.-H. Seo, J. Lee, A. Urban, R. Malik, S. Kang, and G. Ceder, Nature Chemistry 8, 692 (2016).
  * McColl et al. [2022] K. McColl, R. A. House, G. J. Rees, A. G. Squires, S. W. Coles, P. G. Bruce, B. J. Morgan, and M. S. Islam, Nature communications 13, 1 (2022).
  * Wang et al. [2018] A. Wang, S. Kadam, H. Li, S. Shi, and Y. Qi, npj Computational Materials 4, 1 (2018).
  * Hine et al. [2009] N. Hine, K. Frensch, W. Foulkes, and M. Finnis, Physical Review B 79, 024112 (2009).
  * Schwerdtfeger [2011] P. Schwerdtfeger, ChemPhysChem 12, 3143 (2011).
  * Lejaeghere et al. [2016] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, et al., Science 351, aad3000 (2016).
  * Hartwigsen et al. [1998] C. Hartwigsen, S. Gœdecker, and J. Hutter, Physical Review B 58, 3641 (1998).
  * Zhang et al. [2022] M. Zhang, D. A. Kitchaev, Z. Lebens-Higgins, J. Vinckeviciute, M. Zuba, P. J. Reeves, C. P. Grey, M. S. Whittingham, L. F. Piper, A. Van der Ven, et al., Nature Reviews Materials pp. 1–19 (2022).
  * Li et al. [2020] W. Li, E. M. Erickson, and A. Manthiram, Nature Energy 5, 26 (2020).
  * Gent et al. [2017] W. E. Gent, K. Lim, Y. Liang, Q. Li, T. Barnes, S.-J. Ahn, K. H. Stone, M. McIntire, J. Hong, J. H. Song, et al., Nature communications 8, 1 (2017).
  * Lim et al. [2015] J.-M. Lim, D. Kim, Y.-G. Lim, M.-S. Park, Y.-J. Kim, M. Cho, and K. Cho, Journal of Materials Chemistry A 3, 7066 (2015).
  * Eum et al. [2020] D. Eum, B. Kim, S. J. Kim, H. Park, J. Wu, S.-P. Cho, G. Yoon, M. H. Lee, S.-K. Jung, W. Yang, et al., Nature Materials 19, 419 (2020).
  * Born and Oppenheimer [1927] M. Born and J. R. Oppenheimer, Annalen der Physik 84, 457 (1927).
  * Kohanoff [2006] J. Kohanoff, _Electronic structure calculations for solids and molecules: theory and computational methods_ (Cambridge university press, 2006).
  * Ashcroft and Mermin [1976] N. Ashcroft and N. Mermin, _Solid state physics_ (1976).
  * Babbush et al. [2019] R. Babbush, D. W. Berry, J. R. McClean, and H. Neven, npj Quantum Information 5, 1 (2019).
  * Slater [1937] J. C. Slater, Physical Review 51, 846 (1937).
  * Wills et al. [2010] J. M. Wills, M. Alouani, P. Andersson, A. Delin, O. Eriksson, and O. Grechnyev, _Full-Potential Electronic Structure Method: energy and force calculations with density functional and dynamical mean field theory_ , vol. 167 (Springer Science & Business Media, 2010).
  * Herring [1940] C. Herring, Physical Review 57, 1169 (1940).
  * Phillips and Kleinman [1959] J. C. Phillips and L. Kleinman, Physical Review 116, 287 (1959).
  * Antončík [1959] E. Antončík, Journal of Physics and Chemistry of Solids 10, 314 (1959).
  * Troullier and Martins [1990] N. Troullier and J. Martins, Solid State Communications 74, 613 (1990).
  * Hamann et al. [1979] D. Hamann, M. Schlüter, and C. Chiang, Physical Review Letters 43, 1494 (1979).
  * Kerker [1980] G. Kerker, Journal of Physics C: Solid State Physics 13, L189 (1980).
  * Bachelet et al. [1982] G. B. Bachelet, D. R. Hamann, and M. Schlüter, Physical Review B 26, 4199 (1982).
  * Troullier and Martins [1991] N. Troullier and J. L. Martins, Physical Review B 43, 1993 (1991).
  * Louie et al. [1982] S. G. Louie, S. Froyen, and M. L. Cohen, Physical Review B 26, 1738 (1982).
  * Engel et al. [2001] E. Engel, A. Höck, R. Schmid, R. Dreizler, and N. Chetty, Physical Review B 64, 125111 (2001).
  * Martin [2020] R. M. Martin, _Electronic structure: basic theory and practical methods_ (Cambridge university press, 2020).
  * Kleinman and Bylander [1982] L. Kleinman and D. Bylander, Physical Review Letters 48, 1425 (1982).
  * Blöchl [1990] P. E. Blöchl, Physical Review B 41, 5414 (1990).
  * Vanderbilt [1990] D. Vanderbilt, Physical Review B 41, 7892 (1990).
  * Low and Chuang [2019b] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019b).
  * Loaiza et al. [2022] I. Loaiza, A. M. Khah, N. Wiebe, and A. F. Izmaylov, arXiv:2208.08272 (2022).
  * Babbush et al. [2018b] R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. R. McClean, A. Paler, A. Fowler, and H. Neven, Physical Review X 8, 041015 (2018b).
  * Sugisaki et al. [2016] K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, and T. Takui, The Journal of Physical Chemistry A 120, 6459 (2016).
  * Tubman et al. [2018] N. M. Tubman, C. Mejuto-Zaera, J. M. Epstein, D. Hait, D. S. Levine, W. Huggins, Z. Jiang, J. R. McClean, R. Babbush, M. Head-Gordon, et al., arXiv:1809.05523 (2018).
  * Lee et al. [2022] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, et al., arXiv:2208.02199 (2022).
  * McArdle et al. [2022] S. McArdle, A. Gilyén, and M. Berta, arXiv:2210.14892 (2022).
  * Low et al. [2018] G. H. Low, V. Kliuchnikov, and L. Schaeffer, arXiv:1812.00954 (2018).
  * Jain et al. [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, et al., APL materials 1, 011002 (2013).
  * Kitaev and Webb [2008] A. Kitaev and W. A. Webb, arXiv:0801.0342 (2008).
  * Campéon and Yabuuchi [2021] B. D. L. Campéon and N. Yabuuchi, Chemical Physics Reviews 2, 041306 (2021).
  * Jiao et al. [2022] S. Jiao, Q. Li, X. Xiong, X. Yu, H. Li, L. Chen, and X. Huang, Applied Physics Letters 121, 070501 (2022).
  * Yabuuchi [2019] N. Yabuuchi, The Chemical Record 19, 690 (2019).
  * Yabuuchi et al. [2011] N. Yabuuchi, K. Yoshii, S.-T. Myung, I. Nakai, and S. Komaba, Journal of the American Chemical Society 133, 4404 (2011).
  * Wang et al. [2013] R. Wang, X. He, L. He, F. Wang, R. Xiao, L. Gu, H. Li, and L. Chen, Advanced Energy Materials 3, 1358 (2013).
  * Van der Ven et al. [2020] A. Van der Ven, Z. Deng, S. Banerjee, and S. P. Ong, Chemical Reviews 120, 6977 (2020).
  * Chen and Islam [2016] H. Chen and M. S. Islam, Chemistry of Materials 28, 6656 (2016).
  * Mohanty et al. [2014] D. Mohanty, J. Li, D. P. Abraham, A. Huq, E. A. Payzant, D. L. Wood III, and C. Daniel, Chemistry of Materials 26, 6272 (2014).
  * Kleiner et al. [2018] K. Kleiner, B. Strehle, A. R. Baker, S. J. Day, C. C. Tang, I. Buchberger, F.-F. Chesneau, H. A. Gasteiger, and M. Piana, Chemistry of Materials 30, 3656 (2018).
  * Lee and Persson [2014] E. Lee and K. A. Persson, Advanced Energy Materials 4, 1400498 (2014).
  * Li et al. [2014] Q. Li, G. Li, C. Fu, D. Luo, J. Fan, and L. Li, ACS Applied Materials & Interfaces 6, 10330 (2014).
  * Qiu et al. [2016] B. Qiu, M. Zhang, L. Wu, J. Wang, Y. Xia, D. Qian, H. Liu, S. Hy, Y. Chen, K. An, et al., Nature Communications 7, 1 (2016).
  * Shi et al. [2016] J.-L. Shi, J.-N. Zhang, M. He, X.-D. Zhang, Y.-X. Yin, H. Li, Y.-G. Guo, L. Gu, and L.-J. Wan, ACS Applied Materials & Interfaces 8, 20138 (2016).
  * Shin et al. [2018] Y. Shin, W. H. Kan, M. Aykol, J. K. Papp, B. D. McCloskey, G. Chen, and K. A. Persson, Nature Communications 9, 1 (2018).
  * Ku et al. [2018] K. Ku, J. Hong, H. Kim, H. Park, W. M. Seong, S.-K. Jung, G. Yoon, K.-Y. Park, H. Kim, and K. Kang, Advanced Energy Materials 8, 1800606 (2018).
  * Pei et al. [2020] Y. Pei, Q. Chen, M. Wang, B. Li, P. Wang, G. Henkelman, L. Zhen, G. Cao, and C.-Y. Xu, Nano Energy 71, 104644 (2020).
  * Chen et al. [2020] Q. Chen, Y. Pei, H. Chen, Y. Song, L. Zhen, C.-Y. Xu, P. Xiao, and G. Henkelman, Nature Communications 11, 1 (2020).
  * Henkelman et al. [2000] G. Henkelman, B. P. Uberuaga, and H. Jónsson, Journal of Chemical Physics 113, 9901 (2000).
  * Momma and Izumi [2011] K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Perdew et al. [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Giannozzi et al. [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Payne et al. [1992] M. C. Payne, M. P. Teter, D. C. Allan, T. Arias, and a. J. Joannopoulos, Reviews of modern physics 64, 1045 (1992).
  * Blaha et al. [2001] P. Blaha, K. Schwarz, G. K. Madsen, D. Kvasnicka, J. Luitz, et al., An augmented plane wave plus local orbitals program for calculating crystal properties 60 (2001).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , vol. 55 (US Government printing office, 1964).
  * Sharpe et al. [2020] R. Sharpe, R. A. House, M. J. Clarke, D. Förstermann, J.-J. Marie, G. Cibin, K.-J. Zhou, H. Y. Playford, P. G. Bruce, and M. S. Islam, Journal of the American Chemical Society 142, 21799 (2020).
  * Perdew et al. [1992] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Physical Review B 46, 6671 (1992).

Appendix

##  Appendix A Plane wave matrix elements of the pseudopotential operator

In this section we derive the expressions for the matrix elements of the local
and non-local components of the pseudopotential operator described in Section
III.1 in a plane-wave basis.

###  A.1 Matrix elements of the local potential

For an ion located at the coordinates $\bm{R}$, the plane-wave matrix element
of the local potential is defined by the integral

$u_{pq}^{\mathrm{loc}}(\bm{R})=\frac{1}{\Omega}\int
d\bm{r}^{3}u^{\text{loc}}(||\bm{r}-\bm{R}||)e^{i\bm{G}_{\nu}\cdot\bm{r}},$
(72)

where $\bm{G}_{\nu}=\bm{G}_{q}-\bm{G}_{p}$. By changing the variable
$\bm{r}\to\bm{r}+\bm{R}$ we obtain

$u_{pq}^{\mathrm{loc}}(\bm{R})=\frac{1}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\int
d\bm{r}^{3}u^{\text{loc}}(r)e^{i\bm{G}_{\nu}\cdot\bm{r}},$ (73)

with $r=\|\bm{r}\|$. In spherical coordinates, the integral above transforms
as

$\displaystyle u_{pq}^{\mathrm{loc}}(\bm{R})$
$\displaystyle=\frac{2\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\int_{0}^{\infty}dr~{}u^{\text{loc}}(r)\int_{0}^{\pi}d\theta\cos(G_{\nu}r\cos(\theta))~{}r^{2}\sin(\theta).$
(74)

To integrate the angular variable we use

$\displaystyle\mathrm{cos}[G_{\nu}r\mathrm{cos}(\theta)]r^{2}\mathrm{sin}(\theta)=r\frac{d}{d\theta}\left[-\frac{\mathrm{sin}(G_{\nu}r\mathrm{cos}(\theta))}{G_{\nu}}\right],$
(75)

which simplifies Eq. 74 to the integral over the radial variable

$\displaystyle u_{pq}^{\mathrm{loc}}(\bm{R})=\frac{4\pi}{\Omega}$
$\displaystyle
e^{i\bm{G}_{\nu}\cdot\bm{R}}\left[\frac{1}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}u^{\text{loc}}(r)~{}\text{sin}(G_{\nu}r)\right].$
(76)

Note that Eq. 76 is general and it can be used to compute the matrix element
of any given local potential $u^{\text{loc}}(r)$. For the case of the
Hartwigsen-Goedecker-Hutter (HGH) pseudopotential we insert Eq. 24 into Eq. 76
to obtain the following expression:

$\displaystyle u_{pq}^{\mathrm{loc}}(\bm{R})$
$\displaystyle=\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\left[\frac{-Z_{\mathrm{ion}}}{G_{\nu}}\int_{0}^{\infty}dr~{}\text{erf}(\alpha
r)~{}\text{sin}(G_{\nu}r)+\sum_{i=1}^{4}\frac{C_{i}}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}(\sqrt{2}\alpha
r)^{2i-2}~{}e^{-(\alpha r)^{2}}~{}\text{sin}(G_{\nu}r)\right],$ (77)

where $\alpha=\frac{1}{\sqrt{2}r_{\text{loc}}}$. The first integral is
computed as follows:

$\kern-6.0pt-\frac{Z_{\mathrm{ion}}}{G_{\nu}}\int_{0}^{\infty}dr~{}\text{erf}(\alpha
r)~{}\text{sin}(G_{\nu}r)=-\frac{Z_{\mathrm{ion}}e^{-G_{\nu}^{2}\alpha^{-2}/4}}{G_{\nu}^{2}}.$
(78)

To evaluate the second term we compute the integrals:

$\displaystyle\frac{C_{1}}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}e^{-(\alpha
r)^{2}}~{}\text{sin}(G_{\nu}r)=C_{1}\frac{\sqrt{\pi}}{2}r_{\text{loc}}^{3}e^{-(G_{\nu}r_{\text{loc}})^{2}/2},$
(79) $\displaystyle\frac{2C_{2}}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}(\alpha
r)^{2}~{}e^{-(\alpha
r)^{2}}~{}\text{sin}(G_{\nu}r)=C_{2}\frac{\sqrt{\pi}}{2}\left[3r_{\text{loc}}^{3}-r_{\text{loc}}^{5}G_{\nu}^{2}\right]e^{-(G_{\nu}r_{\text{loc}})^{2}/2},$
(80) $\displaystyle\frac{4C_{3}}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}(\alpha
r)^{4}~{}e^{-(\alpha
r)^{2}}~{}\text{sin}(G_{\nu}r)=C_{3}\frac{\sqrt{\pi}}{2}\left[15r_{\text{loc}}^{3}-10r_{\text{loc}}^{5}G_{\nu}^{2}+r_{\text{loc}}^{7}G_{\nu}^{4}\right]e^{-(G_{\nu}r_{\text{loc}})^{2}/2},$
(81) $\displaystyle\frac{8C_{4}}{G_{\nu}}\int_{0}^{\infty}dr~{}r~{}(\alpha
r)^{6}~{}e^{-(\alpha
r)^{2}}~{}\text{sin}(G_{\nu}r)=C_{4}\frac{\sqrt{\pi}}{2}\left[105r_{\text{loc}}^{3}-105r_{\text{loc}}^{5}G_{\nu}^{2}+21r_{\text{loc}}^{7}G_{\nu}^{4}-r_{\text{loc}}^{9}G_{\nu}^{6}\right]e^{-(G_{\nu}r_{\text{loc}})^{2}/2}.$
(82)

Using Eq. 78 and Eqs. (79) to (82) we obtain a final expression for the matrix
elements of the local potential

$\displaystyle
u_{pq}^{\mathrm{loc}}(\bm{R})=\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}~{}e^{-(G_{\nu}r_{\text{loc}})^{2}/2}\bigg{\\{}-\frac{Z_{\mathrm{ion}}}{G_{\nu}^{2}}+\frac{\sqrt{\pi}}{2}\Big{[}$
$\displaystyle
C_{1}r_{\text{loc}}^{3}+C_{2}(3r_{\text{loc}}^{3}-5r_{\text{loc}}^{5}G_{\nu}^{2})+C_{3}(15r_{\text{loc}}^{3}-10r_{\text{loc}}^{5}G_{\nu}^{2}+r_{\text{loc}}^{7}G_{\nu}^{4})$
$\displaystyle+C_{4}(105r_{\text{loc}}^{3}-105r_{\text{loc}}^{5}G_{\nu}^{2}+21r_{\text{loc}}^{7}G_{\nu}^{4}-r_{\text{loc}}^{9}G_{\nu}^{6})\Big{]}\bigg{\\}}.$
(83)

###  A.2 Matrix elements of the non-local potential

The matrix elements of the non-local potential $u^{\text{NL}}$ defined in Eq.
25 are given by

$\displaystyle
u^{\mathrm{NL}}_{pq}=\sum_{lm}\sum_{i,j}~{}\langle\varphi_{p}|\beta_{i}^{(lm)}\rangle~{}B_{ij}^{(l)}~{}\langle\beta_{j}^{(lm)}|\varphi_{q}\rangle,$
(84)

where
$\langle\bm{r}|\beta_{i}^{(lm)}\rangle=\beta_{i}^{(lm)}(\bm{r}-\bm{R})=\beta_{i}^{(l)}(||\bm{r}-\bm{R}||)Y_{lm}(\phi,\theta)$,
and $\bm{R}$ denotes the coordinates of an ion in the material. Here
$Y_{lm}(\phi,\theta)$ denote spherical harmonics. By changing the variable
$\bm{r}\to\bm{r}+\bm{R}$ the overlap integral
$\langle\varphi_{p}|\beta_{i}^{(lm)}\rangle$ writes as

$\langle\varphi_{p}|\beta_{i}^{(lm)}\rangle=e^{-i\bm{G}_{p}\cdot\bm{R}}\int
d\bm{r}\varphi_{p}(\bm{r})\beta_{i}^{(lm)}(\bm{r}).$ (85)

The plane wave function $\varphi_{p}(\bm{r})$ can be expanded in terms of
spherical harmonics $Y_{lm}(\phi,\theta)$ and Bessel functions $j_{l}(G_{p}r)$
[98]

$\varphi_{p}(\bm{r})=\frac{4\pi}{\sqrt{\Omega}}\sum_{lm}i^{l}j_{l}(G_{p}r)Y_{lm}^{*}(\hat{\bm{G}}_{p})Y_{lm}(\hat{\bm{r}}),$
(86)

where $\hat{\bm{r}}$ and $\hat{\bm{G}}_{p}$ are unitary vectors. Using this
expansion we compute the overlap in Eq. 85 as follows:

$\displaystyle\langle\varphi_{p}|\beta_{i}^{(lm)}\rangle$
$\displaystyle=\frac{4\pi}{\sqrt{\Omega}}e^{-i\bm{G}_{p}\cdot\bm{R}}\sum_{l^{\prime}m^{\prime}}(-i)^{l^{\prime}}Y_{l^{\prime}m^{\prime}}(\hat{\bm{G}}_{p})\int_{0}^{\infty}dr~{}r^{2}j_{l^{\prime}}(G_{p}r)\beta_{i}^{(l)}(r)\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta~{}\mathrm{sin}(\theta)Y_{l^{\prime}m^{\prime}}^{*}(\phi,\theta)Y_{lm}(\phi,\theta)$
$\displaystyle=\frac{4\pi}{\sqrt{\Omega}}e^{-i\bm{G}_{p}\cdot\bm{R}}(-i)^{l}Y_{lm}(\hat{\bm{G}}_{p})\int
dr~{}r^{2}j_{l}(G_{p}r)\beta_{i}^{(l)}(r).$ (88)

Similarly, the projection $\langle\beta_{j}^{(lm)}|\varphi_{q}\rangle$ is
given by

$\langle\beta_{j}^{(lm)}|\varphi_{q}\rangle=\frac{4\pi}{\sqrt{\Omega}}e^{i\bm{G}_{q}\cdot\bm{R}}i^{l}Y_{lm}^{*}(\hat{\bm{G}}_{q})\int
dr~{}r^{2}j_{l}(G_{q}r)\beta_{j}^{(l)}(r).$ (89)

Plugging Eqs. (88)-(89) into Eq. 84 we obtain an expression for computing the
matrix elements in Eq. 84 for a given set of projectors $\beta_{i}^{l}(r)$

$\displaystyle u^{\mathrm{NL}}_{pq}(\bm{R})=$
$\displaystyle\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\sum_{l=0}^{l_{\mathrm{max}}}(2l+1)P_{l}(\hat{\bm{G}}_{p}\cdot\hat{\bm{G}}_{q})\sum_{i,j}~{}\langle
p|i\rangle_{l}~{}B_{ij}^{(l)}~{}\langle j|q\rangle_{l},$ (90)

where $\bm{G}_{\nu}=\bm{G}_{q}-\bm{G}_{p}$, $\langle p|i\rangle_{l}$ denotes
the overlap integral over the radial coordinates, and $P_{l}$ is a Legenre
polynomial.

$\langle p|i\rangle_{l}=\int dr~{}r^{2}j_{l}(G_{p}r)\beta_{i}^{(l)}(r).$ (91)

To compute the matrix elements of the non-local component of the Hartwigsen-
Goedecker-Hutter (HGH) pseudopotenial, we use the HGH projectors defined as:

$\beta_{i}^{(l)}(r)=A_{l}^{i}r^{l+2i-2}\mathrm{exp}\left[-\frac{1}{2}\left(\frac{r}{r_{l}}\right)^{2}\right],$
(92)

where the radii $r_{l}$ give the range of the $l$-dependent projectors, and

$A_{l}^{i}=\frac{\sqrt{2}}{r_{l}^{l+(4i-1)/2}\Gamma(l+\frac{4i-1}{2})},$ (93)

with $\Gamma$ denoting the gamma function. For the case of one projector
($i=1$) and $l_{\text{max}}\leq 2$, the overlap integrals defined by Eq. 91
are calculated as:

$\displaystyle\langle p|1\rangle_{0}$
$\displaystyle=\int_{0}^{\infty}dr~{}r^{2}\frac{\mathrm{sin}(G_{p}r)}{(G_{p}r)}~{}\beta_{1}^{(0)}(r)=2r_{0}^{3/2}~{}e^{-(G_{p}r_{0})^{2}/2},$
(94) $\displaystyle\langle p|1\rangle_{1}$
$\displaystyle=\int_{0}^{\infty}dr~{}r^{2}\left[\frac{\mathrm{sin}(G_{p}r)}{(G_{p}r)^{2}}-\frac{\mathrm{cos}(G_{p}r)}{(G_{p}r)}\right]~{}\beta_{1}^{(1)}(r)=\frac{4}{3}r_{1}^{5/2}~{}G_{p}~{}e^{-(G_{p}r_{1})^{2}/2},$
(95) $\displaystyle\langle p|1\rangle_{2}$
$\displaystyle=\int_{0}^{\infty}dr~{}r^{2}\left[\frac{3\mathrm{sin}(G_{p}r)}{(G_{p}r)^{3}}-\frac{\mathrm{sin}(G_{p}r)}{(G_{p}r)}-\frac{3\mathrm{cos}(G_{p}r)}{(G_{p}r)^{2}}\right]~{}\beta_{1}^{(2)}(r)=\frac{8}{15}r_{2}^{7/2}~{}G_{p}^{2}~{}e^{-(G_{p}r_{2})^{2}/2}.$
(96)

Using (94)-(96) we compute the product of the projections
$\bra{p}i\rangle_{l}\bra{i}q\rangle_{l}$ entering Eq. 90 to obtain the final
expression for the matrix elements

$\displaystyle
u^{\mathrm{NL}}_{pq}(\bm{R})=\frac{4\pi}{\Omega}e^{i\bm{G}_{\nu}\cdot\bm{R}}\Bigg{\\{}$
$\displaystyle
4B_{0}r_{0}^{3}~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{0}^{2}/2}+\frac{16B_{1}r_{1}^{5}}{3}~{}(\bm{G}_{p}\cdot\bm{G}_{q})~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{1}^{2}/2}+$
$\displaystyle\left[\frac{32B_{2}r_{2}^{7}}{15}~{}(\bm{G}_{p}\cdot\bm{G}_{q})^{2}+\frac{32B_{2}r_{2}^{7}}{45}~{}(G_{p}G_{q})^{2}\right]e^{-(G_{p}^{2}+G_{q}^{2})r_{2}^{2}/2}\Bigg{\\}},$
(97)

where the coefficient $B_{l}:=B_{11}^{(l)}$ in Eq. 90.

##  Appendix B Structural models of lithium-excess materials

Here we provide more details on how structural models were built for the
lithium-excess materials studied in this work. These models are used for
determining the number of plane waves needed for convergence of the total
energy, computed using density functional theory, and for resource estimation
of the quantum algorithm.

For $\text{Li}_{0.5}\text{MnO}_{3}$, we first modify the position of the Mn
atoms in the pristine structure to generate the spinel-like phase simulated in
Ref [43]. Thus, we delithiated the material following the results reported in
Ref. [84]. The structural model for the LLNMO material corresponds to the most
stable structure of the
$\text{Li}_{0.75}[\text{Li}_{0.17}\text{Ni}_{0.25}\text{Mn}_{0.58}]\text{O}_{2}$
material predicted in Ref. [44]. Similarly, we removed one lithium layer from
the supercell to estimate the resources to compute the site energies reported
in Ref. [44]. Finally, the structural model for the delithiated
$\text{Li}_{0.75}\text{MnO}_{2}\text{F}$ material was obtained by randomly
taking Li ions off the pristine material until we match the Li-ion
concentration of $x=0.75$. The relaxed structure of the fully lithiated
material was taken from Ref. [99].

##  Appendix C Number of plane waves used to perform resources estimation of
the quantum algorithm

The purpose of estimating the number of plane waves is to have an appropriate
number that we can use to perform resource estimations of the quantum
algorithm. To this aim, we estimate the number of plane waves needed to
perform a classical calculation using the structures studied in this work. We
conducted an energy convergence test by changing the cutoff energy and
calculating the total energy at the $\Gamma$-point (Fig. 4). Each point in
Fig. 4a plot corresponds to the total energy value, calculated by using the
Quantum ESPRESSO package, for a given kinetic energy cutoff. Fig. 4b) depicts
a similar analysis conducted using the WIEN2k program, an all-electron
approach, where the cutoff in the reciprocal space is given by the
($R_{\text{MT}}$ and $K_{\text{max}}$) parameter. The number of plane waves
was chosen to ensure that the total energy variation was less than 1 kcal/mol
= 0.043 eV (chemical accuracy).

Figure 4: Convergence of the ground-state energy for the selected materials
computed using DFT at the level of the PBE-GGA exchange-correlation functional
[100]. a) Energies computed using plane waves and the HGH pseudopotentials as
implemented in the Quantum ESPRESSO package. b) Energies obtained for the all-
electron calculations using WIEN2k as a function of the cutoff parameters
$R_{\text{MT}}$ and $K_{\text{max}}$. The number of plane waves used to
perform the resource estimations for each material is indicated with arrows.

##  Appendix D Linear combination of unitaries for the non-local potential

This section provides details on the construction of the linear combination of
unitaries (LCU) decomposition for the non-local term of the pseudopotential
operator. The LCU for the non-local term exploits the projector representation
of the operator, and breaks down each projection $P$ to the sum of unitaries
$\frac{1}{2}\mathbbm{1}-\frac{1}{2}(\mathbbm{1}-2P)$. Below we demonstrate how
this is done in details. Recall

$\displaystyle U_{NL}=$
$\displaystyle\sum_{j=1}^{\eta}\sum_{I=1}^{L}f_{I}(\bm{p},\bm{q})\ket{\bm{p}}_{j}\bra{\bm{q}}_{j},$
(98) $\displaystyle
f_{I}(\bm{p},\bm{q})=\frac{4\pi}{\Omega}e^{i\bm{G}_{\bm{\nu}}\cdot\bm{R}_{I}}\Bigg{\\{}$
$\displaystyle
4r_{0}^{3}B_{0}~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{0}^{2}/2}+\frac{16r_{1}^{5}B_{1}}{3}~{}(\bm{G}_{p}\cdot\bm{G}_{q})~{}e^{-(G_{p}^{2}+G_{q}^{2})r_{1}^{2}/2}$
$\displaystyle+\left[\frac{32r_{2}^{7}B_{2}}{15}~{}(\bm{G}_{p}\cdot\bm{G}_{q})^{2}+\frac{32r_{2}^{7}B_{2}}{45}~{}(G_{p}G_{q})^{2}\right]e^{-(G_{p}^{2}+G_{q}^{2})r_{2}^{2}/2}\Bigg{\\}}.$
(99)

We decompose $U_{NL}$ by focusing on each individual term in the function
$f_{I}$. Define:

$\displaystyle\ket{\Psi_{I,0}}$
$\displaystyle:=\frac{1}{\sqrt{\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}}}}\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}/2}\ket{\bm{p}}.$
(100)

Now for the term in Appendix D including $r_{0}$, we define

$\displaystyle U_{NL,0}$
$\displaystyle:=\frac{4\pi}{\Omega}\sum_{j=1}^{\eta}\sum_{I=1}^{L}\left(\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}}\right)4r_{0}^{3}B_{0}R(\bm{R}_{I})^{\dagger}\ket{\Psi_{I,0}}\bra{\Psi_{I,0}}R(\bm{R}_{I})$
$\displaystyle=\sum_{I=1}^{L}\frac{8\pi
r_{0}^{3}B_{0}\eta(\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}})}{\Omega}\mathbbm{1}+\sum_{j=1}^{\eta}\sum_{I=1}^{L}\frac{-8\pi
r_{0}^{3}B_{0}(\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}})}{\Omega}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,0}}\bra{\Psi_{I,0}})R(\bm{R}_{I}),$
(101)

where $R(\bm{R}_{I})\ket{\bm{p}}=e^{i\bm{G}_{p}\cdot\bm{R}_{I}}\ket{\bm{p}}$
is a phase action and we have suppressed the electron index $j$ for the phase
and for the register in state $\ket{\Psi_{I,0}}$. Thus, in the formula above,
using the same notations of LABEL:eq:lcu_NL in the main text, for $\sigma=0$
we have

$\displaystyle c_{I,0}:=\frac{-8\pi
r_{0}^{3}B_{0}(\sum_{\bm{p}}e^{-G_{p}^{2}r_{0}^{2}})}{\Omega}.$ (102)

The sign of $c_{I,0}$ only depends on $B_{0}$, i.e., on the atomic type of
$I$. For the next term in $U_{NL}$ involving $r_{1}$, we adopt the notation
$\sigma=(1,\omega)$. For $\omega\in\\{x,y,z\\}$, we define

$\displaystyle\ket{\Psi_{I,1,\omega}}$
$\displaystyle:=\frac{1}{\sqrt{\sum_{\bm{p}}G_{p,\omega}^{2}e^{-G_{p}^{2}r_{1}^{2}}}}\sum_{\bm{p}}G_{p,\omega}e^{-G_{p}^{2}r_{1}^{2}/2}\ket{\bm{p}},$
(103)

where $\bm{G}_{p}=(G_{p,1},G_{p,2},G_{p,3})$, and we have the following
decomposition depending on the coordinate $\omega$:

$\displaystyle U_{NL,1,\omega}$
$\displaystyle:=\frac{4\pi}{\Omega}\sum_{j=1}^{\eta}\sum_{I=1}^{L}\left(\sum_{\bm{p}}G_{p,\omega}^{2}e^{-G_{p}^{2}r_{1}^{2}}\right)\frac{16r_{1}^{5}B_{1}}{3}R(\bm{R}_{I})^{\dagger}\ket{\Psi_{I,1,\omega}}\bra{\Psi_{I,1,\omega}}R(\bm{R}_{I})$
$\displaystyle=\sum_{I=1}^{L}\frac{32\pi
r_{1}^{5}B_{1}\eta(\sum_{\bm{p}}G_{p,\omega}^{2}e^{-G_{p}^{2}r_{1}^{2}})}{3\Omega}\mathbbm{1}+\sum_{j=1}^{\eta}\sum_{I=1}^{L}\frac{-32\pi
r_{1}^{5}B_{1}(\sum_{\bm{p}}G_{p,\omega}^{2}e^{-G_{p}^{2}r_{1}^{2}})}{3\Omega}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,1,\omega}}\bra{\Psi_{I,1,\omega}})R(\bm{R}_{I})$
(104)

implying

$\displaystyle c_{I,1,\omega}=\frac{-32\pi
r_{1}^{5}B_{1}(\sum_{\bm{p}}G_{p,\omega}^{2}e^{-G_{p}^{2}r_{1}^{2}})}{3\Omega}.$
(105)

For the term including the scalar $32r_{2}^{7}B_{2}/45$, corresponding to the
index $\sigma=(2,0)$, we define

$\displaystyle\ket{\Psi_{I,2,0}}$
$\displaystyle:=\frac{1}{\sqrt{\sum_{\bm{p}}G_{p}^{4}e^{-G_{p}^{2}r_{2}^{2}}}}\sum_{\bm{p}}G_{p}^{2}e^{-G_{p}^{2}r_{2}^{2}/2}\ket{\bm{p}},$
(106)

and note the following decomposition:

$\displaystyle U_{NL,2,0}$
$\displaystyle:=\frac{4\pi}{\Omega}\sum_{j=1}^{\eta}\sum_{I=1}^{L}(\sum_{\bm{p}}G_{p}^{4}e^{-G_{p}^{2}r_{2}^{2}})\frac{32r_{2}^{7}B_{2}}{45}R(\bm{R}_{I})^{\dagger}\ket{\Psi_{I,2,0}}\bra{\Psi_{I,2,0}}R(\bm{R}_{I})$
$\displaystyle=\sum_{I=1}^{L}\frac{64\pi
r_{2}^{7}B_{2}\eta(\sum_{\bm{p}}G_{p}^{4}e^{-G_{p}^{2}r_{2}^{2}})}{45\Omega}\mathbbm{1}+\sum_{j=1}^{\eta}\sum_{I=1}^{L}\frac{-64\pi
r_{2}^{7}B_{2}(\sum_{\bm{p}}G_{p}^{4}e^{-G_{p}^{2}r_{2}^{2}})}{45\Omega}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,2,0}}\bra{\Psi_{I,2,0}})R(\bm{R}_{I}),$
(107)

implying

$\displaystyle c_{I,2,0}=\frac{-64\pi
r_{2}^{7}B_{2}(\sum_{\bm{p}}G_{p}^{4}e^{-G_{p}^{2}r_{2}^{2}})}{45\Omega}.$
(108)

Finally, for the term including the scalar $\frac{32r_{2}^{7}B_{2}}{15}$ in
Appendix D, corresponding to the index $\sigma=(2,(\omega,\omega^{\prime}))$
where $\omega\neq\omega^{\prime}\in\\{x,y,z\\}$, we first define

$\displaystyle\ket{\Psi_{I,2,(\omega,\omega^{\prime})}}$
$\displaystyle:=\frac{1}{\sqrt{\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})^{2}e^{-G_{p}^{2}r_{2}^{2}}}}\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})e^{-G_{p}^{2}r_{2}^{2}/2}\ket{\bm{p}}$
(109)

and derive the following decomposition

$\displaystyle U_{NL,2,(\omega,\omega^{\prime})}$
$\displaystyle:=\frac{4\pi}{\Omega}\sum_{j=1}^{\eta}\sum_{I=1}^{L}\left(\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})^{2}e^{-G_{p}^{2}r_{2}^{2}}\right)\frac{32r_{2}^{7}B_{2}}{15}\ket{\Psi_{I,2,(\omega,\omega^{\prime})}}\bra{\Psi_{I,2,(\omega,\omega^{\prime})}}$
$\displaystyle=\sum_{I=1}^{L}\frac{64\pi
r_{2}^{7}B_{2}\eta\left(\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})^{2}e^{-G_{p}^{2}r_{2}^{2}}\right)}{15\Omega}\mathbbm{1}+$
$\displaystyle\sum_{j=1}^{\eta}\sum_{I=1}^{L}\frac{-64\pi
r_{2}^{7}B_{2}\left(\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})^{2}e^{-G_{p}^{2}r_{2}^{2}}\right)}{15\Omega}R(\bm{R}_{I})^{\dagger}(\mathbbm{1}-2\ket{\Psi_{I,2,(\omega,\omega^{\prime})}}\bra{\Psi_{I,2,(\omega,\omega^{\prime})}})R(\bm{R}_{I}),$
(110)

implying

$\displaystyle c_{I,2,(\omega,\omega^{\prime})}=\frac{-64\pi
r_{2}^{7}B_{2}\left(\sum_{\bm{p}}(G_{p,\omega}G_{p,\omega^{\prime}})^{2}e^{-G_{p}^{2}r_{2}^{2}}\right)}{15\Omega}.$
(111)

This finishes the LCU for $U_{NL}$:

$\displaystyle
U_{NL}=U_{NL,0}+\sum_{\omega}U_{NL,1,\omega}+U_{NL,2,0}+\sum_{\omega\neq\omega^{\prime}}U_{NL,2,(\omega,\omega^{\prime})}.$
(112)

In each of the above decompositions, there is a multiple of the identity,
which are not considered in the qubitization as we can shift the Hamiltonian
by the appropriate scalar. Throughout the text, we use the alternative
indexing $0\leq\sigma\leq 10$ to index the operators above in the order they
were derived.

| a) |
---|---
| b) |
---|---
Figure 5: ([72, Fig. 1]) (a) The Select $\sum_{x=0}^{N-1}\ket{x}\bra{x}\otimes
X^{\theta_{x}}$ with $N=4$. The symbol $\oslash$ indicates control by a number