# Fault-tolerant quantum simulation of materials using Bloch orbitals

Nicholas C. Rubin<EMAIL_ADDRESS>Google Research, Venice, CA 90291,
United States  Dominic W. Berry<EMAIL_ADDRESS>School of
Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109,
Australia  Fionn D. Malone  Google Research, Venice, CA 90291, United States
Alec F. White  Quantum Simulation Technologies, Inc., Boston, 02135, United
States  Tanuj Khattar  Google Research, Venice, CA 90291, United States  A.
Eugene DePrince, III  Google Research, Venice, CA 90291, United States
Department of Chemistry and Biochemistry, Florida State University,
Tallahassee, FL, USA  Sabrina Sicolo  BASF SE, Carl-Bosch-Strasse 38, 67063
Ludwigshafen, Germany  Michael Kühn  BASF SE, Carl-Bosch-Strasse 38, 67063
Ludwigshafen, Germany  Michael Kaicher  BASF SE, Carl-Bosch-Strasse 38, 67063
Ludwigshafen, Germany  Joonho Lee  Google Research, Venice, CA 90291, United
States  Department of Chemistry, Columbia University, New York, NY, USA  Ryan
Babbush<EMAIL_ADDRESS>Google Research, Venice, CA 90291, United States

###### Abstract

The simulation of chemistry is among the most promising applications of
quantum computing. However, most prior work exploring algorithms for block-
encoding, time-evolving, and sampling in the eigenbasis of electronic
structure Hamiltonians has either focused on modeling finite-sized systems, or
has required a large number of plane wave basis functions. In this work, we
extend methods for quantum simulation with Bloch orbitals constructed from
symmetry-adapted atom-centered orbitals so that one can model periodic ab
initio Hamiltonians using only a modest number of basis functions. We focus on
adapting existing algorithms based on combining qubitization with tensor
factorizations of the Coulomb operator. Significant modifications of those
algorithms are required to obtain an asymptotic speedup leveraging
translational (or, more broadly, Abelian) symmetries. We implement block
encodings using known tensor factorizations and a new Bloch orbital form of
tensor hypercontraction. Finally, we estimate the resources required to deploy
our algorithms to classically challenging model materials relevant to the
chemistry of Lithium Nickel Oxide battery cathodes within the surface code.

###### Contents

  1. I Introduction
  2. II Electronic structure Hamiltonian of materials in Bloch orbitals
    1. II.1 Basis functions and matrix elements
    2. II.2 The second-quantized Hamiltonian
  3. III Qubitization of materials Hamiltonians
    1. III.1 The sparse Hamiltonian representation
    2. III.2 The single-factorization Hamiltonian representation
    3. III.3 The double-factorization Hamiltonian representation
    4. III.4 The tensor hypercontraction Hamiltonian representation
  4. IV Scaling comparison and runtimes for diamond
  5. V Classical and quantum simulations of LNO
    1. V.1 LNO background
    2. V.2 Correlated $k$-point calculations
    3. V.3 Single shot density matrix embedding theory
    4. V.4 Quantum resource estimates for LNO
  6. VI Conclusion
  7. A Sparse representation derivations
    1. A.1 The Pauli operator representation of the one-body term
    2. A.2 One-body correction for sparse case
    3. A.3 Complexity for sparse implementation
  8. B Single-factorization derivations
    1. B.1 One-body correction for single factorization
    2. B.2 Complexity for single-factorized representation
  9. C Double-factorization derivations
    1. C.1 One-body correction
    2. C.2 Complexity of the double-factorized representation
  10. D Tensor hypercontraction derivations
    1. D.1 THC symmetries
    2. D.2 Complexity of the tensor hypercontraction representation
  11. E Correlation diagnostics for LNO
  12. F Classical timing benchmarks
  13. G Generating the THC factors

##  I Introduction

Recently, first quantization quantum algorithms and constant factor resource
estimation analysis for molecular systems [1, 2] have been adapted to
materials [3]. While the first quantization approach using a plane wave
representation is attractive due to the smooth convergence to the continuum
limit [4, 5] a local basis representation such as atom-centered basis sets has
other advantages. Similar to the molecular simulation setting, local basis
functions can be advantageous when describing spatially localized phenomena
such as heterogeneous catalysis or efficiently describing cusps [6].

The desire for systematically improvable electronic structure methods to treat
the many examples of strongly correlated phenomena [7, 8, 9] in the condensed
phase has recently driven the application of ab initio wavefunction theories
to the periodic setting [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].
Standard treatments of symmetry in wavefunction theories [23, 24] can be used
to exploit the translational symmetry of periodic systems, thus enabling the
application of post-Hartree-Fock methods to material systems. Despite these
advantages, classical ab initio treatment of such problems is limited due to
the large simulation cells needed to converge to the thermodynamic limit. This
drawback has further driven the use of embedding theories [25, 26, 27, 28, 29,
30] and downfolding [31]. Naturally, one may ask if fault-tolerant quantum
computers can alleviate the computational burden associated with ab initio
simulation of solids within the local basis framework.

In this paper, we describe how to extend molecular quantum simulation
algorithms of second quantization Hamiltonians represented in local basis sets
to periodic systems using the qubitization framework [32, 33]. Though the
general structure of the algorithms is largely unchanged, introducing
symmetry–i.e. symmetry-adapting the block encodings–requires non-trivial
modifications to realize an improvement in the asymptotic complexity. The
first steps in this direction were taken in Ref. [34] using the “sparse”
Hamiltonian representation. We provide an alternative derivation for block
encodings using this representation and introduce symmetry-adapted block
encodings for three other more performant tensor factorizations of the
Hamiltonian: single factorization (SF), double factorization (DF) and tensor
hypercontraction (THC). The result is orders of magnitude improvement in the
quantum resources required to simulate materials.

For each of the four Hamiltonian representations we describe the origin of the
asymptotic speedup (or lack thereof in one case), provide compiled algorithms
for constant factor resource estimates, and compare the performance to non-
symmetry-adapted block encodings. We note that the derived symmetry-adapted
block encodings apply to any Abelian point group symmetry with minor
modifications. For SF, sparse, and DF the symmetry-adapted block encodings
provide an asymptotic speedup for walk operator construction proportional to
the square root of the number of $k$-points used to sample the Brillouin zone.
For THC, there is no asymptotic improvement due to the linear cost of unary
iteration in the block encoding. Going beyond asymptotic analysis and
compiling to total Toffolis, we find that for DF and THC using symmetry-
adapted block encodings provides no asymptotic speedup over their non-
symmetry-adapted counterparts due to the increased number of applications of
the walk operator for fixed precision phase estimation. DF and THC are
sensitive to the numerical compression of the Hamiltonian, and thus we expect
the number of walk operator applications can be decreased. Furthermore, there
are classical advantages to using the symmetry-adapted block encodings coming
from the reduced classical complexity of representing the Hamiltonian as the
system size is increased towards the thermodynamic limit.

In parallel with recent studies estimating quantum resources required to
simulate high-value molecular targets [35, 32, 36, 37], we estimate the
quantum resources required to simulate an open materials science problem
related to the cathode structure of Lithium Nickel Oxide (LNO) batteries. The
LNO systems are universally observed in the high symmetry $\mathrm{R\bar{3}m}$
structure which is at odds with the predicted Jahn-Teller activity of low-spin
trivalent Ni [38]; more background can be found in Section V.1. This
discrepancy combined with the difficulty of synthesizing pure LNO, the size of
the unit cells [39], and potential strong correlation at the high symmetry
structure [40] makes the LNO problem an interesting application target for
quantum simulation advantage. This realistic problem frames the algorithmic
improvements articulated in this paper and the prospects of the quantum
advantage given modern electronic structure methods. We find that the required
resource estimates for simulating a set of benchmark systems and the LNO
problem before reaching the thermodynamic limit are already substantial. In
fact, the large simulation cells required to converge these calculations to
the thermodynamic limit is ultimately a significant hurdle for ab initio
simulations.

The layout of the rest of the paper is as follows: Section II describes the
atom-centered basis sets and the Hamiltonian that we use, Section III
describes the qubitization algorithm and the origin of the asymptotic speedup
in constructing walk operators using each of the four Hamiltonian
representations. Each subsection is dedicated to a particular Hamiltonian
factorization and describes the qubitization algorithm and how to calculate
associated parameters. Section IV compares all methods and extrapolates
quantum resources required to simulate a diamond crystal converged towards the
thermodynamic limit, and Section V reports the accuracy and correlation
analysis of various electronic structure methods for LNO while providing
estimates of quantum computing resources and runtimes. We close with prospects
for this class of methods.

##  II Electronic structure Hamiltonian of materials in Bloch orbitals

Though plane-wave basis sets are used in most periodic Density Functional
Theory (DFT) calculations, there is a long history of local-basis methods as
well. The use of a localized basis set has a number of advantages over plane
waves: 1) 0D (molecular), 1D, 2D, and 3D systems can be treated on an equal
computational footing, 2) Calculations on low-density systems with large unit
cells can be more efficient[41, 42, 43] 3) Hartree-Fock exchange can be more
efficiently computed in the smaller, local-orbital basis[44, 45, 43, 42] and
4) The local-orbital representations can lower the computational cost of
correlation corrections with a more compact representation of the virtual
space. (1) - (3) have spurred the development of local-orbital DFT and
Hartree-Fock methods with Gaussian orbitals [46, 43, 42, 44] and numerical
atomic orbitals [47], while (4) has been behind recent work to apply
correlated electronic structure theory to periodic solids [10, 11, 13, 14, 15,
16, 17]. In the following subsection we describe the symmetry-adapted periodic
sum of Gaussian-type orbitals used in this work.

###  II.1 Basis functions and matrix elements

A local basis function, $\tilde{\chi}_{p}$, can be adapted to the
translational symmetry of a lattice to form a periodized function

$\displaystyle\chi_{p,\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{T}}e^{i\mathbf{k}\cdot\mathbf{T}}\tilde{\chi}_{p}(\mathbf{r}-\mathbf{T}),$
(1)

where $\mathbf{T}$ represents a lattice translation vector and $\mathbf{k}$ is
a crystal momentum vector lying in the first Brillouin zone. The lattice
momentum $\mathbf{k}$ labels an irreducible representation of the group of
translations defined by the translational symmetry of the material. Functions
of this form are easily verified to be Bloch functions in that

$\chi_{p,\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}u_{p,\mathbf{k}}(\mathbf{r})$
(2)

where $u_{p,\mathbf{k}}(\mathbf{r})$ has the same periodicity as the lattice.

Orbitals are constructed from a linear combination of the underlying Bloch
orbitals,

$\displaystyle\phi_{i\mathbf{k}}(\mathbf{r})=N_{k}^{-1/2}\sum_{p}c_{p,i}(\mathbf{k})\chi_{p,\mathbf{k}}(\mathbf{r}),$
(3)

where $N_{k}$ is the total number of $k$ points. The expansion coefficients,
$c_{p,i}(\mathbf{k})$ are determined from the appropriate periodic self-
consistent field procedure, usually Hartree-Fock or Kohn-Sham DFT. The
resulting orbitals are normally constrained to be orthogonal by convention and
can serve as a basis for representing the second-quantized Hamiltonian.

The matrix elements of a one-electron operator,

$\displaystyle T_{p\mathbf{k}_{p},q\mathbf{k}_{q}}=\int
dr\,\phi^{*}_{p\mathbf{k}_{p}}(r)\mathcal{O}_{1}\phi_{q\mathbf{k}_{q}}(\mathbf{r})$
(4)

are non-zero only when $\mathbf{k}_{p}=\mathbf{k}_{q}$ as long as
$\mathcal{O}_{1}$ has the translational symmetry of the lattice. We can use a
similar strategy to derive the structure of the two-electron integrals which
are given by

$\displaystyle
V_{p\mathbf{k}_{p},q\mathbf{k}_{q},r\mathbf{k}_{r},s\mathbf{k}_{s}}$
$\displaystyle=\int\int
dr_{1}\,dr_{2}\,\phi_{p\mathbf{k}_{p}}^{*}(\mathbf{r}_{1})\phi_{q\mathbf{k}_{q}}(\mathbf{r}_{1})\mathcal{O}_{2}\phi_{r\mathbf{k}_{r}}^{*}(\mathbf{r}_{2})\phi_{s\mathbf{k}_{s}}(\mathbf{r}_{2}).$
(5)

The translational symmetry of the Bloch orbitals implies the 2-electron
operator $\mathcal{O}_{2}$ matrix elements can only be nonzero when
$\left(\mathbf{k}_{p}+\mathbf{k}_{r}-\mathbf{k}_{q}-\mathbf{k}_{s}\right)=\mathbf{G}$
where $\mathbf{G}$ is a reciprocal lattice vector. We note that this
expression for nonzero matrix elements by symmetry is a specific instance of
the more general expression. More generally, given a group $\mathfrak{g}$ with
its irreducible representations labeled by $\\{\Gamma_{i}\\}$, the two-
electron integral is nonzero by symmetry whenever
$\Gamma_{p}\otimes\Gamma_{q}\otimes\Gamma_{r}\otimes\Gamma_{s}$ contains the
complete symmetric representation [23]. For periodic systems $\mathfrak{g}$ is
the set of translational symmetries. Despite this sparsity, the evaluation of
the nonzero matrix elements for all basis functions is often a major
computational bottleneck whenever local basis sets are used.

Local orbitals provide a more compact representation than plane waves, so
fewer basis functions are needed. Unfortunately, there are $O(N_{k}^{3}N^{4})$
generally nonzero two-electron matrix elements for $N_{k}$ $k$-points and $N$
basis functions in the primitive cell. For very large calculations locality
can be exploited to yield asymptotically linear-scaling DFT methods [48, 49].
Linear scaling Hartree-Fock is also possible for insulators [50]. This linear
regime is almost never reached in practice, and it is usually advantageous to
instead reduce the cost by tensor factorization.

Though our discussion has been thus far general with regard to the choice of
local basis functions, Gaussian basis functions are by far the most popular
choice in molecular calculations, and crystalline Gaussian orbitals are also a
popular choice for periodic calculations. This popularity is due to the
existence of analytic formulas which allow for fast, numerically exact
evaluation of the matrix elements of most common operators. Despite the
existence of efficient numerical techniques, the large number of two-electron
integrals that must be evaluated in periodic calculations requires a more
efficient procedure. Traditionally, this is accomplished with the Gaussian
plane wave (GPW) method [51, 41] which only requires storage of
$O(N_{k}^{2}N^{2}n_{\mathrm{pw}})$ integrals where $n_{\mathrm{pw}}$ is the
number of plane waves used to evaluate the integrals. In molecular
calculations, the most common decomposition is called the resolution of the
identity (RI) or sometimes density fitting (DF) [52, 53, 54, 55]. This
procedure requires the storage of $O(N_{k}^{2}N^{2}n_{\mathrm{aux}})$
integrals where $n_{\mathrm{aux}}$ is the size of the auxiliary basis set.
Both the GPW and the RI method can be considered as density fitting approaches
where the former uses a plane-wave fitting basis and the latter uses a
Gaussian fitting basis. For this reason, the RI approach is often called
“Gaussian density fitting” (GDF) in the context of periodic calculations [56,
57, 58, 59, 60].

The two-electron integral tensor can be further factorized into a product of
five two-index tensors as was done in the tensor hypercontraction (THC) method
of Martinez and coworkers [61, 62, 63]. Factorizations of this form are most
useful for correlated methods where they have the potential to lower the
computational scaling. In this work we present a translational symmetry-
adapted form of the tensor hypercontraction for the two-electron integral
tensors of periodic systems.

###  II.2 The second-quantized Hamiltonian

We can express the second-quantized electronic structure Hamiltonian as

$\displaystyle H$ $\displaystyle=H_{1}+H_{2}\,,$ (6) $\displaystyle H_{1}$
$\displaystyle=\sum_{\sigma}\sum_{\mathbf{k}}\sum_{pq}h_{p\mathbf{k},q\mathbf{k}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}\,,$
(7) $\displaystyle h_{p\mathbf{k},q\mathbf{k}}$
$\displaystyle=T_{p\mathbf{k},q\mathbf{k}}-\frac{1}{2}\sum_{r,\mathbf{Q}}V_{p\mathbf{k},r\mathbf{Q},r\mathbf{Q},q\mathbf{k}}\,,$
(8) $\displaystyle H_{2}$
$\displaystyle=\frac{1}{2}\sum_{\sigma,\tau}\sum_{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}\sum_{pqrs}V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}\,.$
(9)

We first introduce summation limits for each symbol as we will commonly use
short hand summation formulas to indicate multiple sums. For each variable
$\\{p,q,r,s\\}$ summation is performed over the range $[0,N/2-1]$ indexing the
spatial orbital or band, $\\{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}\\}$
summation is performed over the Brillouin zone ($\mathcal{B}\\!Z$) at a set
number of $k$-points of which there are $N_{k}$, and $\\{\sigma,\tau\\}$ are
electron spin variables and summed over $\\{\uparrow,\downarrow\\}$. Non-
modular differences of $\mathbf{k}$, $\mathbf{Q}$, and $\mathbf{k}^{\prime}$
span twice the Brillouin zone. Because $V$ needs to be indexed by values in
the Brillouin zone, we use modular subtraction indicated by ${\ominus}$. That
is, if the number of points in each dimension is $N_{x},N_{y},N_{z}$, we
perform subtraction modulo $N_{x},N_{y},N_{z}$ in each direction,
respectively.

The Hamiltonian is generally complex Hermitian with four-fold symmetry of the
two-electron integrals.111We note that the following generic complex Coulomb
integral symmetries are present
$V_{p\mathbf{k}_{p},q\mathbf{k}_{q},r\mathbf{k}_{r},s\mathbf{k}_{s}}=V_{r\mathbf{k}_{r},s\mathbf{k}_{s},p\mathbf{k}_{p},q\mathbf{k}_{q}}=V_{q\mathbf{k}_{q},p\mathbf{k}_{p},s\mathbf{k}_{s},r\mathbf{k}_{r}}^{*}=V_{s\mathbf{k}_{s},r\mathbf{k}_{r},q\mathbf{k}_{q},p\mathbf{k}_{p}}^{*}$
(10) from integration index relabeling and complex conjugation. In the
following sections we demonstrate how the sparse structure of the two-electron
integral tensor affects the scaling of block encoding the Hamiltonian for
implementation of qubitized quantum walk oracles. The cost of qubitization is
greatly affected by the representational freedom of the the underlying
Hamiltonian expressed as a linear combination of unitaries. We demonstrate how
to construct the sparse, single-factorization (SF), double-factorization (DF),
and tensor-hypercontraction (THC) integral decompositions of Bloch orbital
Hamiltonians and cost out simulations for a variety of materials.

For all algorithms, we will make a comparison to the case of a $\Gamma$-point
calculation on a supercell composed of $N_{k}$ primitive cells in the geometry
described by the $k$-point sampling. This allows us to directly observe the
proposed speedup due to symmetry-adapting. To demonstrate the scaling of
symmetry-adapted block encoding, we estimate quantum simulation resource
requirements for the series of systems listed in Table 1. Range-separated
density fitting [60] is used to construct integrals with Dunning type
correlation-consistent basis sets [64] and the Goedecker-Teter-Hutter (GTH)
family of pseudopotentials for Hartree-Fock [65]. For each Hamiltonian,
cutoffs for the factorization are selected so that the Møller-Plesset second
order perturbation theory (MP2) error in the total energy is below one
milliHartree per cell or formula unit depending on the system. While prior
works used coupled-cluster theory, MP2 is used here for computational
efficiency.

System | Structure | Atoms in Cell | Lattice Parameters | spin-orbitals cc-pVDZ | spin-orbitals cc-pVTZ
---|---|---|---|---|---
C | diamond | 2 | 3.567 [66] | 52 | 116
Si | diamond | 2 | 5.43 [66] | 52 | 116
BN | zinc blende | 2 | 3.616 [66] | 52 | 116
LiCl | rocksalt | 2 | 5.106 [67] | 52 | 98
AlN | wurzite | 4 | (a) 3.11 (c) 4.981 [66] | 104 | 220
Li | bcc | 2 | 3.51 [68] | 52 | 80
Al | fcc | 2 | 4.0479 [69] | 52 | 104
Table 1: Crystal structures lattice parameters used for the systems studied in
this work. The lattice parameters were chosen to be at or near their
experimental equilibrium values.

##  III Qubitization of materials Hamiltonians

Similar to fault tolerant resource estimates for molecular systems represented
in second quantization [36, 32, 70, 71, 37], we compare the number of logical
qubits and number of Toffoli gates required to implement phase estimation on
unitaries that use block encoding [72] and qubitization [33] to encode the
Hamiltonian spectrum in a Szegedy walk operator [73] for various linear
combination of unitaries (LCU) [74] representations of the Hamiltonian. All
LCUs represent the Hamiltonian as

$\displaystyle H=\sum_{\ell=1}^{L}\omega_{\ell}U_{\ell}$ (11)

where $\omega_{\ell}\in\mathbb{R}$, $\omega_{\ell}\geq 0$, and $U_{\ell}$ is a
unitary operator. One can then construct the operators

$\displaystyle\mathrm{PREPARE}|0\rangle^{\otimes\log(L)}\mapsto\sum_{\ell=1}^{L}\sqrt{\frac{\omega_{l}}{\lambda}}|\ell\rangle\equiv|\mathcal{L}\rangle$
(12) $\displaystyle\mathrm{SELECT}|\ell\rangle|\psi\rangle\mapsto|\ell\rangle
U_{\ell}|\psi\rangle$ (13)
$\displaystyle\lambda=\sum_{\ell=1}^{L}\omega_{\ell}$ (14)

where $|\psi\rangle$ is the system register, and $|\ell\rangle$ is an ancilla
register used to index each term in the LCU. The walk operator constructed
from select and a reflection operator built from prepare,
$R=2|\mathcal{L}\rangle\langle\mathcal{L}|\otimes\mathbb{1}-\mathbb{1}$, has
eigenvalues proportional $e^{\pm i\arccos{E_{n}/\lambda}}$ where $E_{n}$ is an
eigenvalue of the Hamiltonian in Eq. (11).

It was shown in References [70] and [32] when ensuring that select is self-
inverse, only the reflection operator $R$ needs to be controlled on the
ancilla for phase estimation and not select. Therefore, the Toffoli cost of
phase estimating the walk operator scales as

$\displaystyle\left\lceil{\frac{\pi\lambda}{2\epsilon_{\mathrm{PEA}}}}\right\rceil\left(C_{S}+C_{P}+C_{P^{\dagger}}+\log(L)\right)$
(15)

where $C_{S}$ is the cost for implementing the select oracle and $C_{P}$ is
the cost for implementing the prepare oracle, $C_{P^{\dagger}}$ is the cost
for the inverse prepare oracle, and $\epsilon_{\mathrm{PEA}}$ is the target
precision for phase estimation. Thus the main costs for sampling from the
eigenspectrum of a second quantized operator are the costs to implement
select, prepare, and prepare†. These costs need to be multiplied by a factor
proportional to $\lambda/\epsilon_{\rm PEA}$ for the number of walk steps
needed for phase estimation. Note that when computing intensive quantities,
such as the the energy per cell, the $\lambda$ factor is scaled by $1/N_{k}$.

The particular choice of LCU changes all of these costs. Prior works have
investigated the resource requirements for simulating molecules with four
different LCUs. While all these methods can be used without modification in
supercell calculations at the $\Gamma$-point, the construction of molecular
select and prepare do not exploit any symmetries and are not applicable away
from the $\Gamma$-point – e.g. at the Baldereschi point [75].

The leading costs in constructing select and prepare for second quantized
Hamiltonians is the circuit primitive that functions similar to a read-only-
memory (ROM) called QROM. The QROM primitive is a gadget that takes a memory
address, potentially in superposition, and outputs data, also potentially in
superposition. There are currently two variations of QROM that have different
costs; traditional QROM that has linear Toffoli complexity when outputting $L$
items with any amount of data associated with each item, and advanced QROM
(called QROAM) with reduced non-Clifford complexity [76]. It uses a select-
swap circuit construction with Toffoli cost

$\displaystyle\left\lceil{\frac{L}{k}}\right\rceil+m(k-1)$ (16)

where $k$ is a power of $2$ for outputting $L$ items of data where each item
of data is $m$ bits long. The notation $k$ here for an integer should not be
confused with $\mathbf{k}$ for the crystal momentum vector. It needs $m(k-1)$
ancillas, so increases the logical ancilla count in exchange for reduced
Toffoli complexity. When $L>m$ this function is minimized by selecting
$k\approx\sqrt{L/m}$ and thus the Toffoli and ancilla cost generically go as
$\mathcal{O}(\sqrt{Lm})$. It is also possible to adjust $k$ to reduce the
ancilla count while increasing the Toffoli count. Having QROAM output the
minimal amount of information to represent the Hamiltonian is at the core of
the $\sqrt{N_{k}}$ improvements we derive in many of the block encodings. We
will also demonstrate that for all LCUs the lowest scaling can be linear in
the Bloch orbital basis size $\mathcal{O}(N_{k}N)$ due to the requirement to
perform unary iteration at least once over the entire basis.

Another primitive that becomes the dominant cost in constructing symmetry-
adapted select is the multiplexed-controlled swap between two registers. The
controlled swap between two registers uses unary iteration [70] on $L$ items
to swap $M$ elements between two registers at the cost of $\mathcal{O}(LM)$
Toffolis. For simulating materials this primitive is commonly encountered when
swapping all band indices with a particular irreducible representation label,
or $k$-point, into a working register at a cost of $\mathcal{O}(N_{k}N)$. The
necessity of coherently moving data thus puts a limit on the total savings one
can achieve by leveraging Abelian symmetries. The cost of moving data must be
weighed against the benefits, which we describe in each section below. In
Table 2 we summarize the space complexity, in terms of logical qubits, and
time complexity, in terms of Toffolis of the four LCUs when considering
translational symmetry on the primitive cell and without (denoted as SC for
supercell).

Table 2: Generically, qubitized quantum walks scale as ${\mathcal{O}}(\sqrt{\Gamma})$ in space and ${\mathcal{O}}(\lambda\sqrt{\Gamma}/\epsilon)$ in time where $\Gamma$ is the amount of information required to specify the Hamiltonian within a particular representation. For double factorization, $\Xi$ is the sum of the average rank of the second factorization which is expected to scale as $\mathcal{O}(N_{k}N)$, which is the number of orbitals in the primitive cell or bands. $\tilde{\Xi}$ is the average rank of the second factorization in the supercell calculation and is also expected to scale as $\mathcal{O}(N_{k}N)$. The tilde on $\mathcal{O}$ is used to account for logarithmic factors, and can include variables not explicitly given in the scaling. $\lambda$ for each LCU is different and is denoted as a subscript indicating the LCU type and if it $\lambda$ for the supercell version. Representation | Qubits | Toffoli Complexity | SC Qubits | SC Toffoli
---|---|---|---|---
sparse | $\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{2})$ | $\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{2}\lambda_{\mathrm{sparse}}/\epsilon)$ | $\widetilde{\mathcal{O}}(N_{k}^{2}N^{2})$ | $\widetilde{\mathcal{O}}(N_{k}^{2}N^{2}\lambda_{\mathrm{sparse},\mathrm{SC}}/\epsilon)$
SF | $\widetilde{\mathcal{O}}(N_{k}N^{3/2})$ | $\widetilde{\mathcal{O}}(N_{k}N^{3/2}\lambda_{\mathrm{SF}}/\epsilon)$ | $\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{3/2})$ | $\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{3/2}\lambda_{\mathrm{SF},\mathrm{SC}}/\epsilon)$
DF | $\widetilde{\mathcal{O}}(\sqrt{N_{k}}N\sqrt{\Xi})$ | $\widetilde{\mathcal{O}}(\sqrt{N_{k}}N\sqrt{\Xi}\lambda_{\mathrm{DF}}/\epsilon)$ | $\widetilde{\mathcal{O}}(N_{k}N\sqrt{\tilde{\Xi}})$ | $\widetilde{\mathcal{O}}(N_{k}N\sqrt{\tilde{\Xi}}\lambda_{\mathrm{DF},\mathrm{SC}}/\epsilon)$
THC | $\widetilde{\mathcal{O}}(N_{k}N)$ | $\widetilde{\mathcal{O}}(N_{k}N\lambda_{\mathrm{THC}}/\epsilon)$ | $\widetilde{\mathcal{O}}(N_{k}N)$ | $\widetilde{\mathcal{O}}(N_{k}N\lambda_{\mathrm{THC},\mathrm{SC}}/\epsilon)$

For the sparse LCU, exploiting primitive cell translational symmetry reduces
the amount of symmetry unique information in the Hamiltonian by a factor of
$N_{k}$, which translates to a reduction of $\sqrt{N_{k}}$ savings in Toffoli
complexity and ancilla complexity. For sparse prepare, the square root savings
originates from the QROAM cost of outputting “alt” and “keep” values for the
coherent alias sampling component of the state preparation. For sparse select
controlled application of all Pauli terms has linear cost in the basis size
$\mathcal{O}(N_{k}N)$ and is not the dominant cost. The supercell calculation
does not exploit the $k$-point symmetry unique non-zero coefficients of the
Hamiltonian and thus has worse scaling.

The single factorization LCU leverages the fact that the Coulomb integral
tensor is positive semidefinite and can be written in a quadratic form. For
molecular systems without symmetry–i.e. $C1$ symmetry, the factorization
results in a three-tensor where there are two orbital indices and one
auxiliary index that scales as the number of orbitals in the system [77]. For
simulations where orbitals now have point group symmetry labels, such as
$k$-points, each three-tensor factor can now be arranged into a five-tensor;
two symmetry labels (irreducible representation labels), two-band labels
(orbital labels), and one auxiliary index which still scales with the number
of bands due to density fitting of the cell periodic part of the density [78,
60]. Thus the origin of the $\sqrt{N_{k}}$ improvement for the symmetry-
adapted block encoding with a single-factorization LCU lies in the fact that
the auxiliary index has $N_{k}$ lower scaling in comparison to a supercell
variation where the Cholesky factorization or density fitting is performed on
the entire supercell two-electron integral tensor. The single factorization
algorithm is also dominated by the QROM cost of prepare–of which there are two
state preparations. The inner state preparation [32] for the $k$-point
symmetry-adapted algorithm requires outputting $\mathcal{O}(N_{k}^{2}N^{3})$
to be used in the state preparation leading to $\mathcal{O}(N_{k}N^{3/2})$
Toffoli and qubit complexity. Contrasting this to the supercell calculation,
we see a $\sqrt{N_{k}}$ savings due to the fact that the inner state
preparation requires only $N_{k}^{2}$ information and not $N_{k}^{3}$
information. We elaborate on this point further in Section III.2. select is
implemented in a similar fashion to sparse, scaling as $\mathcal{O}(N_{k}N)$,
and is not a dominant cost.

The double factorization LCU represents the Hamiltonian in a series of non-
orthogonal bases and leverages that a linear combination of ladder operators
can be constructed by a similarity transform of a single fermionic ladder
operator, or Majorana operator, by a unitary generated by a quadratic
fermionic Hamiltonian. In the molecular case, the dominant cost for these
algorithms is the QROM to output the rotations for the similarity transform
and implementing the basis rotations with the programmable gate array circuit
primitive [36, 32] for select. When taking advantage of primitive cell
symmetry we reduce the amount of data needed to be output by QROM by $N_{k}$,
which results in a $\sqrt{N_{k}}$ savings in the Toffoli complexity. Because
we are using advanced QROM, this output size advantage is also observed in the
logical qubit requirements. As mentioned previously, computing the total
number of Toffolis requires scaling the walk operator cost by a linear
function of $\lambda$. We find that using a canonical orbtial basis set, the
total Toffoli cost is higher than the commensurate supercell total Toffoli
cost because $\lambda$ for the symmetry-adapted case increases. The origin of
the increase is related to a reduced variational freedom when selecting non-
orthogonal bases and is further discussed in Section III.3.

Finally, in the THC LCU there is no asymptotic speedup because the molecular
algorithm had the lowest possible scaling for second quantized algorithms.
This stems from the fact that even iterating over the basis once with unary
iteration to apply an operator indexed by basis element has a Toffoli cost of
$\mathcal{O}(N_{k}N)$. As we will discuss in Section III.4 our symmetry-
adapted algorithm offers other benefits such as enabling the classical
precomputation of the THC factors by exploiting symmetry and lowering the
number of controlled rotations.

We now describe the Hamiltonian factorization used in each LCU, the
calculation of $\lambda$ associated with each Hamiltonian factorization, and
outline the construction of the qubitization oracles. Detailed compilations
are provided for each LCU in the Appendices. In each section we provide
numerical evidence that the symmetry-adapted oracles have the reported scaling
by plotting the Toffoli requirements to synthesize select $+$ prepare $+$
prepare-1 compared against the number of $k$-points sampled ($N_{k}$).

###  III.1 The sparse Hamiltonian representation

In the “sparse” method the Hamiltonian described in Eq. (8) and Eq. (9) is
directly translated to Pauli operators which form the LCU. Under the Jordan-
Wigner transformation, we take

$\displaystyle a_{p\mathbf{k}\sigma}$
$\displaystyle\mapsto\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})/2,$
(17) $\displaystyle a^{\dagger}_{p\mathbf{k}\sigma}$
$\displaystyle\mapsto\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})/2,$
(18)

where the notation $\vec{Z}$ is being used to indicate that there is a string
of $Z$ operators on qubits up to (not including) that on which
$X_{p\mathbf{k}\sigma}$ or $Y_{p\mathbf{k}\sigma}$ acts upon. This requires a
choice of ordering for the qubits indexed by $p$, $\mathbf{k}$, and $\sigma$.
We need only apply the string of $Z$ operators for the same value of $\sigma$,
because we always have matching annihilation and creation operators for the
same spin $\sigma$ (so any $Z$ gates on the other spin would cancel). We also
adopt a convention that the ordering of qubits for the Jordan-Wigner
transformation takes $\mathbf{k}$ as the more significant bits, with qubits
for all $p$ with a given $\mathbf{k}$ grouped together. For most of the
discussion we will not need to explicitly consider this ordering.

With the Jordan-Wigner transform the one-body component of the Hamiltonian
takes on the form

$\displaystyle H_{1}$
$\displaystyle=\frac{i}{4}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}\sum_{p,q=1}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}\right\\}$
$\displaystyle\quad+\frac{i}{4}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}\sum_{p,q=1}{\rm
Im}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}+\sum_{\mathbf{k}}\sum_{p=1}h_{p\mathbf{k},p\mathbf{k}}\openone.$
(19)

We provide the full derivation for this expression in Appendix A.1. To derive
the two-body operator LCU we use only complex conjugation symmetry in contrast
to the molecular derivation that used eight-fold symmetry. The two-body
Hamiltonian can be written as

$\displaystyle
H_{2}=\frac{1}{4}\sum_{\sigma,\tau\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}^{N_{k}}\sum_{p,q,r,s=1}^{N/2}$
$\displaystyle\left[V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}\right.$
$\displaystyle\left.+V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\right],$
(20)

where ${\ominus}$ indicates modular subtraction as defined above. In the case
where $\mathbf{Q}\neq 0$ or $p\neq q$ and $r\neq s$, we can move the creation
and annihilation operators using the fermionic anticommutation relations to
give the term on the second line as

$V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}\,.$
(21)

The Jordan-Wigner representation then gives the expression in square brackets
in Eq. (20) as

$\displaystyle\frac{1}{16}\left\\{V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+iY_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})][\vec{Z}(X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}-iY_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau})][\vec{Z}(X_{s\mathbf{k}^{\prime}\tau}+iY_{s\mathbf{k}^{\prime}\tau})]\right.$
$\displaystyle\left.+V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}[\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-iY_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})][\vec{Z}(X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}+iY_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau})][\vec{Z}(X_{s\mathbf{k}^{\prime}\tau}-iY_{s\mathbf{k}^{\prime}\tau})]\right\\}.$
(22)

Then we can separate Eq. (III.1) into real and imaginary components as

$\displaystyle\frac{1}{8}\left\\{{\rm
Re}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})[\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p,\mathbf{k},\sigma}\vec{Z}Y_{q,\mathbf{k}{\ominus}\mathbf{Q},\sigma})][\vec{Z}X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}X_{s\mathbf{k}^{\prime}\tau}+\vec{Z}Y_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}Y_{s\mathbf{k}^{\prime}\tau}]\right.$
$\displaystyle-{\rm
Re}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})[\vec{Z}Y_{p,\mathbf{k},\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})][\vec{Z}Y_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}X_{s\mathbf{k}^{\prime}\tau}-\vec{Z}X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}Y_{s\mathbf{k}^{\prime}\tau}]$
$\displaystyle+{\rm
Im}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})[\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})][\vec{Z}X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}X_{s\mathbf{k}^{\prime}\tau}+\vec{Z}Y_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}Y_{s\mathbf{k}^{\prime}\tau}]$
$\displaystyle\left.+{\rm
Im}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})[\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})][\vec{Z}Y_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}X_{s\mathbf{k}^{\prime}\tau}-\vec{Z}X_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}\vec{Z}Y_{s\mathbf{k}^{\prime}\tau}]\right\\}.$
(23)

In accounting for cases where $\mathbf{Q}=0$ with $p=q$ or $r=s$, the same
expression is obtained, but there are also one-body terms obtained. These
result in a total one-body operator

$\displaystyle\tilde{H}_{1}=H_{1}+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p,q=1}^{N/2}\left(\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}\,.$
(24)

A full derivation of this expression can be found in Appendix A.2.

Using the representation of the one-body and two-body operators as Pauli
operators, we have a linear combination of unitaries form. The $\lambda$
associated with this LCU is

$\displaystyle\lambda$ $\displaystyle=\lambda_{\tilde{H}_{1}}+\lambda_{H_{2}}$
(25) $\displaystyle\lambda_{\tilde{H}_{1}}$
$\displaystyle=\sum_{\mathbf{k}}\sum_{pq}\left\\{\left|\mathrm{Re}[h_{p\mathbf{k},q\mathbf{k}}]+\mathrm{Re}\left[\sum_{\mathbf{k}^{\prime},r}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right]\right|+\left|\mathrm{Im}[h_{p\mathbf{k},q\mathbf{k}}]+\mathrm{Im}\left[\sum_{\mathbf{k}^{\prime},r}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right]\right|\right\\}$
(26) $\displaystyle\lambda_{H_{2}}$
$\displaystyle=\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}}\sum_{pqrs}\left\\{\left|{\rm
Re}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})\right|+\left|{\rm
Im}(V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}})\right|\right\\}.$
(27)

In determining $\lambda_{\tilde{H}_{1}}$ there is a factor of $2$ due to the
summation over spin $\sigma$ and then a factor of $2$ accounting for the fact
that each expression in braces in Eq. (III.1) is the sum of two different
Pauli strings. As a result these factors have cancelled the original $1/4$
prefactor. In the expression for $\lambda_{\tilde{H}_{1}}$ we have also summed
over the native one-body terms and the contributions from the two-body terms.
For $\lambda_{H_{2}}$ we had a factor of $1/8$ in Eq. (III.1), which is
multiplied by the factor of $1/4$ in Eq. (20). The two sums over spin $\sigma$
and $\tau$ give a factor of $4$. Then for each of the real and imaginary parts
in Eq. (III.1) there were sums over 8 Pauli strings, giving a factor of 8. As
a result these factors have also cancelled in the expression for
$\lambda_{H_{2}}$. Note that there is a factor of 2 between this expression
and that in [32], even when we just consider $V$ that is real. The reason is
that in Ref. [32] there was eight-fold symmetry, where here we only have four-
fold symmetry. That is, here we have symmetry when simultaneously swapping the
pairs $p,q$ and $r,s$, whereas in Ref. [32] there are two symmetries from
swapping $p,q$ or $r,s$ on their own. That meant it was possible to express
the Hamiltonian as in Eq. (A2) of that work, then in Eq. (A3) of that work the
Jordan-Wigner mapping was used in the form

$a_{p\sigma}^{\dagger}a_{q\sigma}+a_{q\sigma}^{\dagger}a_{p\sigma}\mapsto\frac{X_{p\sigma}\vec{Z}X_{q\sigma}+Y_{p\sigma}\vec{Z}Y_{q\sigma}}{2}.$
(28)

In this mapping there has been a cancellation of half the Pauli strings, which
results in $\lambda$ being reduced by a factor of 2. Here we only have four-
fold symmetry, so the value of $\lambda$ for the two-body term is a factor of
2 larger than that in [32].

In order to implement the Hamiltonian as a linear combination of unitaries,
the first step is to perform a state preparation on
$\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime},p,q,r,s$. This state preparation
corresponds to the sum, then we will perform controlled operations for each of
the operators in Eq. (III.1). The state preparation is applied using coherent
alias sampling as described in [71]. Because there are multiple variables that
the state needs to be prepared over, it is convenient to use the QROM to
output “ind” values as well as “alt” and “keep” values. Both “ind” and “alt”
give values of all variables
$\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime},p,q,r,s$. Then an inequality test
is performed between keep and an equal superposition state, and the result is
used to control a swap between ind and alt.

There are both real and imaginary values of
$V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}$,
so we also include a qubit to distinguish between these values in the state
preparation. We also do not prepare all values of
$\\{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}\\}$.
There is symmetry in swapping $p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q})$
with $r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}$, or
simultaneously $p\mathbf{k}$ with $q(\mathbf{k}{\ominus}\mathbf{Q})$ and
$r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})$ with $s\mathbf{k}^{\prime}$. Only
those values of
$\\{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}\\}$
that give unique values of $V$ will be prepared. Then the full range can be
obtained by using qubits to control these swaps. There is also a complex
conjugate needed in the symmetry, which can be applied with a Clifford gate.

The dominant complexity in the preparation comes from the QROM. The number of
items of data is $\mathcal{O}(N_{k}^{3}N^{4})$, and by using the advanced form
of QROM the complexity can be made approximately the square root of the total
amount of data (number of items of data times the size of each). The size of
each item of data is logarithmic in $N_{k}$ and $N$ as well as the allowable
error. Therefore, the scaling of the complexity can be given ignoring these
logarithmic parts as $\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{2})$.

To describe the controlled operations needed in order to implement the
operation as in Eq. (III.1), we need to account for the fact that there are
two lines for each of the real and imaginary components of $V$. In addition,
for each line in Eq. (III.1) there is a product of two factors, each of which
is a sum of two terms. To describe the linear combination of unitaries we
therefore introduce three more qubits.

  * •

The first is used to distinguish between the two lines for each of the real
and imaginary parts in Eq. (III.1).

  * •

The second distinguishes between the two terms in the first set of square
brackets.

  * •

The third distinguishes between the two terms in the second set of square
brackets.

When implementing the controlled operations, we perform four operations of the
form of $\vec{Z}X$ or $\vec{Z}Y$, with $X$ or $Y$ being applied on target
qubits indexed by $p\mathbf{k}$, $q(\mathbf{k}{\ominus}\mathbf{Q})$ and so
forth. These Pauli strings are applied using the approach of [32], but in this
case there is the additional complication that we need to select between $X$
or $Y$. This selection can be performed simply by performing the controlled
Pauli string twice, once for $X$ and once for $Y$. The complexity is
proportional to $N_{k}N$, which is trivial compared to the complexity of the
state preparation. The choice of whether $X$ or $Y$ is performed depends on
the value of the three qubits selecting between the terms, as well as the
qubit selecting between the real and imaginary parts. The processing of these
qubits to determine the appropriate choice of $X$ or $Y$ can be performed with
a trivial number of gates.

The last part to consider is how the implementation of the one-body part of
the Hamiltonian is integrated with the implementation of the two-body part. In
the state preparation, amplitudes corresponding to the real and imaginary
parts of $h_{p\mathbf{k},q\mathbf{k}}$ will be produced, as well as a qubit
selecting between the one- and two-body parts. That qubit will be used to also
select between the choice of $X$ and $Y$. For the one-body part there is a
product of only two of the Pauli strings, so the other two will not be applied
at all for the one-body part. See Appendix A.3 for a more detailed description
of the implementation.

In Figure 1 we plot the Toffoli complexity to implement select $+$ prepare $+$
prepare-1 for simulating the aforementioned sample systems using a symmetry-
adapted select and prepare at different Monkhorst-Pack grids and different
number of bands (cc-pVDZ and cc-pVTZ). We compare the symmetry-adapted
calculations to supercell calculations using the same select and prepare. The
supercell calculations do not explicitly take into account the symmetry of the
primitive cell in the full simulation cell. To demonstrate the symmetry-
adapted Hartree-Fock orbitals do not appreciably change the overall scaling
with respect to a supercell calculation we plot the total $\lambda$ for the
supercell calculation (which reruns Hartree-Fock on the supercell) and the
symmetry-adpated version.

(a)(b) Figure 1: (a) Sparse Toffoli step complexity versus the number of
$k$-points for systems in Table 1 using the cc-pVDZ and cc-pVTZ basis set
$\Gamma$-centered Monkhorst-Pack grids of size [1, 1, 1] to [3, 3, 3]. Each
point is a single system described at a particular basis set and $k$-mesh
where the threshold for zeroing each nonzero two-electron integral coefficient
is determined by MP2 as described earlier. The scaling for implementing the
block encoding is shown in the legend. To isolate the $N_{k}$ scaling behavior
we divide the Toffoli step complexity by the square of the number of basis
functions $(N^{2})$. For supercell we expect a scaling going as
$\mathcal{O}(N_{k}^{2})$ and for symmetry-adapted block encodings we expect a
scaling going as $\mathcal{O}(N_{k}^{1.5})$. The ideal symmetry-adapted
scaling is not reached due to finite size effects which are further discussed
in Figure 2. (b) Total $\lambda$ for the symmetry-adapted version (denoted UC)
and the supercell calculation without explicit primitive cell symmetry
(denoted SC) demonstrating no deterioration of $\lambda$ by symmetry-adapting.
Similar to (a) all points are a particular system from the benchmark set in a
fixed basis and $k$-mesh.

Though Figure 1 indicates some computational advantage for the symmetry-
adapted case the expected $\sqrt{N_{k}}$ improvement over the supercell case
is not easily observed. The cost of the sparse method largely depends on the
number of nonzero elements of $V$, which is generically expected to go as
$\mathcal{O}(N_{k}^{3}N^{4})$ for the symmetry-adapted case and
$\mathcal{O}(N_{k}^{4}N^{4})$ for the supercell case. The scaling ultimately
depends on the number of nonzero elements in each block of integrals (indexed
by three-momentum indices), which we expect to be independent of supercell
size $N_{k}$. For Diamond we plot this dependence in Figure 2 and demonstrate
that convergence is slow and there is a strong $N_{k}$ dependence in the
number of nonzero elements in each two-electron integral block. This
dependence makes observing the improvement in Toffoli cost for symmetry-
adapted oracles difficult in the low $N_{k}$ regime.

Figure 2: $N_{k}$ dependence on the number of non-zero elements in each two-
electron integral block indexed by irrep. labels for Diamond in a single-zeta-
valence basis with a $k$-mesh shifted to (1/8, 1/8, 1/8) of the simulation
cell. In the thermodynamic limit this value should be independent of $N_{k}$
and thus both supercell (SC) and symmetry-adapted (UC) should have no
correlation with $N_{k}$–i.e. a slope of zero. The $N_{k}$ dependence for
small $N_{k}$ makes it difficult to observe the asymptotic improvements from
symmetry-adapting the sparse qubitization oracles.

###  III.2 The single-factorization Hamiltonian representation

For the “single-factorization” method, the Cholesky decomposition of the
2-electron integral tensor can be applied iteratively or the factorized forms
can be directly recovered from a density fitted representation of the atomic
orbital integral. The quadratic representation of the two-electron integral
tensor is

$V_{p\mathbf{k}_{p},q\mathbf{k}_{q},r\mathbf{k}_{r},s\mathbf{k}_{s}}=\sum_{n}L_{p\mathbf{k}_{p}q\mathbf{k}_{q},n}L_{s\mathbf{k}_{s}r\mathbf{k}_{r},n}^{*}$
(29)

where $\mathbf{k}_{p}+\mathbf{k}_{r}=\mathbf{k}_{q}+\mathbf{k}_{s}$ modulo a
reciprocal lattice vector $\mathbf{G}$, or
$\mathbf{k}_{p}-\mathbf{k}_{q}-(\mathbf{k}_{s}-\mathbf{k}_{r})=\mathbf{G}$. We
can identify
$\mathbf{k}_{p}-\mathbf{k}_{q}=\mathbf{Q}+\mathbf{G}=\mathbf{k}_{s}-\mathbf{k}_{r}$.
Thus the two-body interaction operator can be written as

$\displaystyle\hat{H}_{2}^{\prime}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{Q}}^{N_{k}}\sum_{n}^{M}\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\left(\sum_{\tau\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}^{\prime}}^{N_{k}}\sum_{rs}^{N/2}L^{*}_{s\mathbf{k}^{\prime}r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),n}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}\right).$
(30)

Due to the reduced symmetry of the complex valued two-electron integral tensor
we take additional steps to form Hermitian operators which can be expressed as
Pauli operators under the Jordan-Wigner transform. We express each one-body
operator in the product of particle-conserving one-body operators forming the
two-electron operator as

$\hat{\rho}_{n}(\mathbf{Q})=\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right),\qquad\hat{\rho}^{\dagger}_{n}(\mathbf{Q})=\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}L^{*}_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right).$
(31)

We now take a linear combination of $\hat{\rho}_{n}(\mathbf{Q})$ to form
Hermitian operators and represent our two-electron integral operator as a sum
of squares of Hermitian operators that are amenable to the approach for the
qubitization of one-body sparse operators via a linear combination of
unitaries. These operators are denoted $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$ and are defined as

$\displaystyle\hat{A}_{n}(\mathbf{Q})=\frac{1}{2}(\hat{\rho}_{n}(\mathbf{Q})+\hat{\rho}^{\dagger}_{n}(\mathbf{Q})),$
(32)
$\displaystyle\hat{B}_{n}(\mathbf{Q})=\frac{i}{2}(\hat{\rho}_{n}(\mathbf{Q})-\hat{\rho}^{\dagger}_{n}(\mathbf{Q})),$
(33)

to give

$\hat{H}^{\prime}_{2}=\frac{1}{2}\sum_{\mathbf{Q}}^{N_{k}}\sum_{n}^{M}\left(\hat{A}^{2}_{n}(\mathbf{Q})+\hat{B}^{2}_{n}(\mathbf{Q})\right).$
(34)

We have taken advantage of the translational symmetry by performing the sum
over $\mathbf{Q}$ outside the squares of $\hat{A}$ and $\hat{B}$, which
reduces the amount of information needed in the representation. In the case
$\mathbf{Q}\neq 0$ we can write $\hat{A}_{n}$ as

$\displaystyle\hat{A}_{n}(\mathbf{Q}\neq 0)$
$\displaystyle=\frac{1}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}\left(L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}^{*}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right)$
$\displaystyle=\frac{1}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]\left(a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right)$
$\displaystyle\quad+\frac{i}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]\left(a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right).$
(35)

Applying the Jordan-Wigner representation then gives

$\displaystyle\hat{A}_{n}(\mathbf{Q}\neq
0)=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}$
$\displaystyle\left(\frac{i{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right.$
$\displaystyle+\left.\frac{i{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right).$
(36)

The same reasoning can be performed for $\hat{B}_{n}(\mathbf{Q}\neq 0)$, which
gives the plus and minus signs between
$a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}$
and $a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}}$ in
Eq. (35) reversed, so the roles of the real and imaginary parts are reversed.
As a result we obtain

$\displaystyle\hat{B}_{n}(\mathbf{Q}\neq
0)=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}$
$\displaystyle\left(\frac{i{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right.$
$\displaystyle+\left.\frac{i{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right).$
(37)

Accounting for the cases with $\mathbf{Q}=0$, we may use the same expressions
with an extra identity, which yields a one-body correction when squaring. We
show in Appendix B.1 that this results in the total one-body operator

$\displaystyle\tilde{H}_{1}$
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p,q=1}^{N/2}\left(h_{p\mathbf{k},q\mathbf{k}}+\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}$
(38)

as before. Therefore the associated $\lambda$ is again
$\lambda_{\tilde{H}_{1}}$ as given in Eq. (25). The $\lambda$ for the two-body
term is then

$\displaystyle\lambda_{V}=\frac{1}{2}\sum_{\mathbf{Q}}\sum_{n}^{M}\left(\sum_{\mathbf{k},pq}(|{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|+|{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|)\right)^{2}.$ (39)

This expression can be obtained by first summing the absolute values of the
weights in the linear combination of unitaries for $\hat{A}$ and $\hat{B}$ to
give

$\sum_{\mathbf{k},pq}(|{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|+|{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|).$ (40)

This is obtained by noting that the sum over the spin gives a factor of 2, and
there are two unitary operators for each of the real and imaginary parts;
together these cancel the factor of 4. Then this expression is squared for
each of $\hat{A}$ and $\hat{B}$, and there is a sum over $\mathbf{Q}$ and $n$
in Eq. (34). A further factor of $1/2$ is obtained because we use amplitude
amplification on each operator as described in [32], thus giving our
expression for $\lambda_{V}$.

Next we describe the method to block encode the Hamiltonian in this single-
factorized representation. The key idea is to perform a state preparation over
$\mathbf{Q}$ and $n$, then block encode the squares of
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$ using a single step of
oblivious amplitude amplification (which saves a factor of 2 for the value of
$\lambda$). That is, we perform block encodings of $\hat{A}_{n}(\mathbf{Q})$
and $\hat{B}_{n}(\mathbf{Q})$, reflect on an ancilla register, then apply the
block encodings again.

For the initial state preparation on $\mathbf{Q},n$, the number of items of
data is $MN_{k}+1$, where the $+1$ is for the one-body part of the
Hamiltonian. This state preparation is via coherent alias sampling, so the
dominant cost is from the QROM needed to output ind and alt values. That has
complexity scaling as $\widetilde{\mathcal{O}}(\sqrt{MN_{k}})$, where the
tilde accounts for the size of the items of data.

For both $\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$ we have
weightings according to the real and imaginary parts of
$L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}$, but the difference is in
what operations are performed in the sum. Therefore, for an LCU block
encoding, the state preparation step may be identical between
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$. For each value of
$\mathbf{Q},n$, the number of unique values of $\mathbf{k},p,q$ to consider is
$N_{k}N^{2}/4$. Unlike the supercell case we cannot take advantage of symmetry
between $p$ and $q$, because we have $p\mathbf{k}$ and
$q(\mathbf{k}{\ominus}\mathbf{Q})$. The relation between $\mathbf{k}$ and
$\mathbf{k}{\ominus}\mathbf{Q}$ is governed by the value of $\mathbf{Q}$ which
is given in the outer sum, and so we cannot exchange $p$ and $q$. There is a
further factor of $2$ for the number of items of data, because both real and
imaginary parts are needed.

Accounting for the values of $\mathbf{Q},n$, the total number of items of data
that must be output by the QROM used in the state preparation is
$(MN_{k}+1)N_{k}N^{2}/2=\mathcal{O}(N_{k}^{2}N^{3})$, given that $M$ scales as
$\mathcal{O}(N)$. Again because the size of the items of data is logarithmic,
this gives a complexity $\widetilde{\mathcal{O}}(N_{k}N^{3/2})$. In contrast,
in the supercell calculation, each $\hat{A}_{n}$ and $\hat{B}_{n}$ would have
$\mathcal{O}(N_{k}^{2}N^{2})$ entries, and the rank would be
$\mathcal{O}(N_{k}N)$, for a total number of items of data
$\mathcal{O}(N_{k}^{3}N^{3})$. That would give a complexity
$\widetilde{\mathcal{O}}(N_{k}^{3/2}N^{3/2})$, so there is a factor of
$\sqrt{N_{k}}$ improvement obtained by taking advantage of the symmetry.

In the state preparation we only prepare $p,q$ for $p\leq q$, and the full
range of values should be produced using a swap controlled by an ancilla
register. A further subtlety in the implementation as compared to prior work
is that the complex conjugate is needed as well. This may be implemented using
a sign flip on the qubit indicating the imaginary part, so it is just a
Clifford gate.

A major difference is in the selection of operations for
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$. We see that there are
two steps where we need to apply an operation of the form $\vec{Z}X$ or
$\vec{Z}Y$, and the choice of $X$ or $Y$. The selection of where the $X$ or
$Y$ is applied (indicated by the subscript) can be implemented in the standard
way. The choice of whether $X$ or $Y$ is applied depends on four qubits.

  1. 1.

The qubit selecting between the one- and two-body parts.

  2. 2.

A qubit selecting between $A$ and $B$, which can simply be prepared in an
equal superposition using a Hadamard because there are equal weightings
between these operators.

  3. 3.

A qubit selecting between the real and imaginary parts of
$L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}$, which was prepared in the
state preparation.

  4. 4.

A qubit selecting between the two terms shown above in each line of the
expressions for $\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$. This
qubit can also be prepared using a Hadamard.

Using a trivial number of operations on these qubits we can determine whether
it is $X$ or $Y$ that needs to be performed. The cost of the controlled
unitary is doubled because we apply a controlled $\vec{Z}X$ and a controlled
$\vec{Z}Y$, but this cost is trivial compared to the state preparation cost so
has little effect on the overall complexity.

A further subtlety in the implementation is that in the second implementation
of $\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$, we simply use the
qubit flagging the one-body part to control whether the Pauli string
$\vec{Z}X$ or $\vec{Z}Y$ is applied at all. This ensures that the square is
not obtained for the one-body part. For a more in-depth explanation of the
implementation, see the circuit diagram in Figure 3 and the explanation in
Appendix B.2.

Figure 3: The circuit for performing the state preparation and controlled
operations for the single factorization approach. The register labelled $\ell$
is a contiguous register for $\mathbf{Q},n$, with $\mathbf{Q}$ also output in
the state preparation. The inner state preparation uses $\ell$ as a control.
Then the minus on $\mathbf{Q}$ controlled by $\mathbf{k}$ is to compute
$\mathbf{k}{\ominus}\mathbf{Q}$, with that register reset to $\mathbf{Q}$
later with a controlled addition. A qubit is used to swap $p$ and $q$ to
generate that symmetry, and another is used to swap the spin up and spin down
components of the system so the selection only need act on the spin down
component. The qubits labelled Re/Im, $A/B$, and “term” are the qubits
selecting the real versus imaginary parts, $A$ versus $B$, and the two terms
in each line of the Hamiltonian. These correspond to $b_{1},b_{2},b_{3}$ in
Appendix B.2, and the Toffoli and CNOT gates are used so that the “term” qubit
can be used to select whether the Paul string with $X$ or $Y$ is applied. The
selection is performed twice, once for each of the Pauli strings and so is
controlled by $\mathbf{k}{\ominus}\mathbf{Q}$ and $q$ the first time, then
$\mathbf{k}$ and $p$ the second time. A controlled phase between Re/Im and
“term” (also controlled by the success flag qubits) is used to generate the
correct sign for the term. The block encoding of $A/B$ is performed twice,
with the reflection on the ancilla qubits in the middle generating the step of
oblivious amplitude amplification. The third register flags that we have the
two-body part of the Hamiltonian, and is used to control the block encoding of
$A/B$ the second time to ensure it is not performed for the one-body part.
(a)(b) Figure 4: (a) Number of $k$-points verses Toffoli cost to implement the
block encoding for the single factorization LCU evaluated for the benchmark
systems listed in Table 1 described using the cc-pVDZ and cc-pVTZ basis sets
and $\Gamma$-centered Monkhorst-Pack grids of size [1, 1, 1] to [3, 3, 3].
Each point is a single system described at a particular basis set and $k$-mesh
where the range of the auxiliary index of the Cholesky factorization is
selected to produce two-electron integrals corresponding to an MP2 error of
one 1 milliHartree per unit cell with respect to an untruncated auxiliary
index range. We divide the Toffoli complexity for implementing select \+
prepare \+ prepare-1 by $N^{3/2}$, which is the shared scaling in the number
of bands. The different scaling in number of $k$-points becomes clear: $N_{k}$
for symmetry-adapted block encodings and $N_{k}^{3/2}$ for supercell non-
symmetry-adapted block encodings. We observe similar behavior for qubit count
and for plotting oracle Toffoli complexity versus the number of bands. (b) The
value of $\lambda$ per unit cell ($\lambda/N_{k}$) as a function of the total
system size $NN_{k}$ for the same systems described with the same cutoffs used
in (a).

Figure 4 demonstrates the $\sqrt{N_{k}}$ improvement in constructing the walk
operator by symmetry-adapting. Even for small $N_{k}$ there is a clear
separation between the cost of supercell (SC) and symmetry-adapted oracles
that agrees with the theoretical scalings of 1.5 and 1.0, respectively.

###  III.3 The double-factorization Hamiltonian representation

In the sparse and SF LCU approaches we have found that there is a factor of
$\sqrt{N_{k}}$ savings in Toffoli costs and logical qubit costs for symmetry-
adapted block encoding constructions over their non-symmetry-adapted
counterparts (supercell calculations). The double-factorization (DF)
representation continues this trend, though the origin of the speedup is
different. In the double-factorization circuits each unitary of the LCU is a
rank-one one-body operator that can be thought of as the outer product of two
vectors of ladder operators, where each vector of ladder operators is obtained
by a Givens rotation with multiqubit control based on other indices. First
notice that for SF there is $\mathcal{O}(N_{k}^{2}N^{3})$ data to output to
specify the Hamiltonian via the Cholesky factors. The factors come from two
momentum indices $\mathbf{k},\mathbf{Q}$, two band indices $p,q$, and one
auxiliary index $n$. In this section we demonstrate that by using a workspace
register to apply Givens rotations to pairs of band indices,
$\\{\mathbf{k},\mathbf{k}{\ominus}\mathbf{Q}\\}$, the complexity of the DF LCU
can also be improved by a factor of $\sqrt{N_{k}}$ over supercell
calculations.

To construct the DF LCU, we will separate $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$ out into sums over $\mathbf{k}$. To express this,
instead of having $\rho_{n}(\mathbf{Q})$, we define
$\rho_{n}(\mathbf{Q},\mathbf{k})$

$\hat{\rho}_{n}(\mathbf{Q},\mathbf{k})=\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{pq}^{N/2}L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right),\qquad\hat{\rho}^{\dagger}_{n}(\mathbf{Q},\mathbf{k})=\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{pq}^{N/2}L^{*}_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right)$
(41)

so then the Hermitian one-body operators that are squared to form the two body
part of the Hamiltonian are

$\displaystyle\hat{A}_{n}(\mathbf{Q})=\sum_{\mathbf{k}}\frac{1}{2}(\hat{\rho}_{n}(\mathbf{Q},\mathbf{k})+\hat{\rho}^{\dagger}_{n}(\mathbf{Q},\mathbf{k}))\,,$
(42)
$\displaystyle\hat{B}_{n}(\mathbf{Q})=\sum_{\mathbf{k}}\frac{i}{2}(\hat{\rho}_{n}(\mathbf{Q},\mathbf{k})-\hat{\rho}^{\dagger}_{n}(\mathbf{Q},\mathbf{k}))\,.$
(43)

Just as in the single factorization case, we have the two-body part of the
Hamiltonian

$\hat{H}^{\prime}_{2}=\frac{1}{2}\sum_{\mathbf{Q}}^{N_{k}}\sum_{n}^{M}\left(\hat{A}^{2}_{n}(\mathbf{Q})+\hat{B}^{2}_{n}(\mathbf{Q})\right).$
(44)

We can write $\hat{A}_{n}(\mathbf{Q})$ as

$\displaystyle\hat{A}_{n}(\mathbf{Q})=\sum_{\mathbf{k}}\left[U^{A}_{n}(\mathbf{Q},\mathbf{k})\left(\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{p}(\mathbf{Q},n,\mathbf{k})n_{p\mathbf{k}\sigma}\right)U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right]$
(45)

where the basis rotation unitary $U_{n}(\mathbf{Q},\mathbf{k})$ acts on
orbitals indexed by $\mathbf{k}$ and $\mathbf{k}{\ominus}\mathbf{Q}$,
$\Xi_{\mathbf{Q},n,\mathbf{k},A}$ corresponds to a rank cutoff for $A$, and
$f_{p}^{A}(\mathbf{Q},n,\mathbf{k})$ is the eigenvalue of the one body
operator that is diagonalized by $U_{n}(\mathbf{Q},\mathbf{k})$. The
expression for $\hat{B}_{n}(\mathbf{Q})$ is similar, and we use
$\Xi_{\mathbf{Q},n,\mathbf{k},B}$ to denote the rank cutoff.

In practice, for the implementation we would apply a different basis rotation
for each individual value of $p$. As explained by [32], when doing that the
number of Givens rotations needed only corresponds to the number of orbitals
it is acting upon, instead of the square. Here we have two momentum modes
$\\{\mathbf{k},\mathbf{k}{\ominus}\mathbf{Q}\\}$ with $N$ orbitals for each,
suggesting there should be $2N$. However, there is no mixture between the
different spin states indexed by $\sigma$, so that gives the number of
orbitals as $N$.

To quantify the amount of information needed to specify the rotations for the
Hamiltonian, there is a $\mathbf{Q}$ summation, $n$ summation, $\mathbf{k}$
summation, $p$ summation, and we need to specify $N$ Givens rotations for
each. In turn, each Givens rotation needs two angles. The total data here
therefore scales as

$\widetilde{\mathcal{O}}\left(N\sum_{\mathbf{Q},n,\mathbf{k}}(\Xi_{\mathbf{Q},n,k,A}+\Xi_{\mathbf{Q},n,k,B})\right),$
(46)

where a factor of $N$ comes from the number of Givens rotations, and the tilde
accounts for the bits of precision given for the rotations. By analogy with
the supercell case, it is convenient to define an average rank

$\Xi\coloneqq\frac{1}{2N_{k}M}\sum_{\mathbf{Q},n,\mathbf{k}}(\Xi_{\mathbf{Q},n,k,A}+\Xi_{\mathbf{Q},n,k,B}).$
(47)

This has division by $N_{k}$ for the $\mathbf{Q}$ sum and $M$ for the $n$ sum,
but no division by a factor accounting for $\mathbf{k}$. That is, it is the
average rank for each value of $\mathbf{Q}$ and $n$, with the sum over
$\mathbf{k}$ regarded as part of the rank. Then it is most closely analogous
to the rank in the supercell case, and it is found that it similarly scales as
$\mathcal{O}(N_{k}N)$. In terms of $\Xi$, the scaling of the amount of data
can be given as $\widetilde{\mathcal{O}}(N_{k}N^{2}\Xi)$, using
$M=\mathcal{O}(N)$.

Next we describe in general terms how to perform the block encoding of the
linear combination of unitaries, with the full explanation in Appendix C. As
in the case of single factorisation, the general principle is to perform state
preparation over $\mathbf{Q}$ and $n$, then block encode the squares of
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$ using oblivious
amplitude amplification. The difference is that $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$ are now block encoded in a factorized form. In more
detail, the key parts are as follows.

  1. 1.

Perform a state preparation over $\mathbf{Q}$ and $n$, as well as a qubit
distinguishing between $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$. Using the advanced QROM, the complexity of this
state preparation scales approximately as the square root of the number of
items of data, so as $\widetilde{O}(\sqrt{N_{k}N})$. The tilde accounts for
logarithmic factors from the size of the output. For convenience here we use a
contiguous register for combined values of $\mathbf{Q}$ and $n$.

  2. 2.

Apply a QROM which outputs the value of $\mathbf{Q}$, as well as an offset
needed for the contiguous register needed in the state preparation for
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$.

  3. 3.

Perform the inner state preparation over $\mathbf{k}$ and $p$. Here the number
of items of data is ${\mathcal{O}}(N_{k}N\Xi)$, accounting for the sums over
$\mathbf{Q},n,\mathbf{k},p$. The complexity via advanced QROM is approximately
the square root of this quantity, $\widetilde{\mathcal{O}}(\sqrt{N_{k}N\Xi})$.

  4. 4.

Apply the QROM again to output the rotation angles for the Givens rotations
needed for the basis rotation. This time the size of the output scales as
$\mathcal{O}(N)$, so the complexity of the QROM scales as
$\widetilde{\mathcal{O}}(N\sqrt{N_{k}\Xi})$. This is the dominating term in
the complexity.

  5. 5.

Use control qubits to swap the system registers into the correct location.
This is done first controlled by a qubit labelling the spin, $\sigma$, which
is similar to what was done in prior work. The new feature here is that
registers containing $\mathbf{k}$ and $\mathbf{k}{\ominus}\mathbf{Q}$ are also
used to swap system registers into $N$ target qubits.

  6. 6.

Apply the Givens rotations on these $N$ target qubits. The complexity here
only scales as $\widetilde{\mathcal{O}}(N)$, so is smaller than in the other
steps.

  7. 7.

Apply a controlled $Z$ for part of the number operator. This comes from
representing the number operator as $(\openone-Z)/2$ and combining the
identity with the one-body part of the Hamiltonian.

  8. 8.

Invert the Givens rotations, controlled swaps, QROM for the Givens rotations,
and state preparation over $\mathbf{k},p$. The complexities here are similar
to those in the previous steps, but the complexities for QROM erasure are
reduced. This completes the block encoding of $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$.

  9. 9.

Perform a reflection on the ancilla qubits used for the state preparation on
$\mathbf{k},p$. This is needed for the oblivious amplitude amplification.

  10. 10.

Perform steps 3 to 8 again for a second block encoding of
$\hat{A}_{n}(\mathbf{Q})$ and $\hat{B}_{n}(\mathbf{Q})$. This together with
the reflection gives a step of oblivious amplitude amplification, and
therefore the squares of $\hat{A}_{n}(\mathbf{Q})$ and
$\hat{B}_{n}(\mathbf{Q})$.

  11. 11.

Invert the QROM from step 2. This has reduced complexity because it is an
erasure.

  12. 12.

Invert the state preparation from step 1.

A quantum circuit for the procedure is shown in Figure 5. This is similar to
Figure 16 in [32], except it is including the extra parts needed in order to
account for the momentum $\mathbf{k}$ used here. In particular, $\ell$ shown
here is a contiguous register for $\mathbf{Q},n$, and $p$ shown in the diagram
is actually a contiguous register for $\mathbf{k},p$. The values of
$\mathbf{Q}$ and $\mathbf{k}$ need to be output via QROM after the state
preparations. Then $\mathbf{k}{\ominus}\mathbf{Q}$ is computed and
$\mathbf{k}$ and $\mathbf{k}{\ominus}\mathbf{Q}$ are used to swap the required
part of the system register into $N$ target qubits where we apply the Givens
rotations.

Figure 5: The circuit for performing the state preparation and controlled
operations for the double factorization approach. The register labelled $\ell$
is a contiguous register for preparing $\mathbf{Q},n$, with $\mathbf{Q}$
output next in the QROM. The register labelled $p$ is actually a contiguous
register for preparing both $\mathbf{k}$ and $p$. The value of $\mathbf{k}$ is
output in the next step together with the rotations. Then the minus on
$\mathbf{Q}$ controlled by $\mathbf{k}$ is to compute
$\mathbf{k}{\ominus}\mathbf{Q}$. The “$\times$” on
$\mathinner{|{\psi_{\downarrow}}\rangle}$ controlled by the $\mathbf{Q}$ and
$\mathbf{k}$ registers indicates that these registers are used to swap the
required part of the system register into $N$ target qubits that the Givens
rotations $R$ act upon.

The lambda value for the Hamiltonian can be calculated by determining the
total L1-norm of the coefficients of the unitaries used to represent the
Hamiltonian. To determine this norm, note first that the number operator is
replaced with $(\openone-Z)/2$, and the identity is combined with the one-body
part of the Hamiltonian. For $\hat{A}_{n}(\mathbf{Q})$, what is implemented
therefore corresponds to

$\displaystyle-\frac{1}{2}\sum_{\mathbf{k}}\left[U^{A}_{n}(\mathbf{Q},\mathbf{k})\left(\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{p}(\mathbf{Q},n,\mathbf{k})Z_{p\mathbf{k}\sigma}\right)U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right].$
(48)

Summing the absolute values of coefficients here gives

$\displaystyle\sum_{\mathbf{k}}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}|f^{A}_{p}(\mathbf{Q},n,\mathbf{k})|,$
(49)

where the sum over the spin $\sigma$ has given a factor of 2 which canceled
the factor of $1/2$. In implementing the square of $\hat{A}_{n}(\mathbf{Q})$
we use oblivious amplitude amplification, which provides a factor of $1/2$ to
$\lambda$. Combining this with the $1/2$ in the definition of
$\hat{H}^{\prime}_{2}$ gives $1/4$, and combining with the contribution from
$\hat{B}_{n}(\mathbf{Q})$ then gives

$\displaystyle\lambda_{\mathrm{DF},2}=\frac{1}{4}\sum_{\mathbf{Q},n}\left[\left(\sum_{\mathbf{k},p}^{N_{k}\Xi_{\mathbf{Q},n,\mathbf{k},A}}|f^{A}_{n}(p,\mathbf{Q},\mathbf{k})|\right)^{2}+\left(\sum_{\mathbf{k},p}^{N_{k}\Xi_{\mathbf{Q},n,\mathbf{k},B}}|f^{B}_{n}(p,\mathbf{Q},\mathbf{k})|\right)^{2}\right],$
(50)

where the superscript $B$ on $f$ indicates the corresponding quantity for
$\hat{B}_{n}(\mathbf{Q})$.

The one-body Hamiltonian is adjusted by the one-body term arising from the
identity in the representation of the number operator in the two-body
Hamiltonian. This yields an effective one-body Hamiltonian (see Appendix C)

$\displaystyle
H_{1}^{\prime}=\sum_{\mathbf{k},p,q,\sigma}\left(h_{p\mathbf{k},q\mathbf{k}}+\sum_{\mathbf{k}^{\prime},r}V_{r\mathbf{k}^{\prime},r\mathbf{k}^{\prime},q\mathbf{k},p\mathbf{k}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(51)

We can rewrite this as

$\displaystyle
H_{1}^{\prime}=\sum_{\mathbf{k},\sigma}\left[U^{C}(\mathbf{k})\left(\sum_{p}^{N/2}\lambda_{\mathbf{k},p}n_{p\mathbf{k}\sigma}\right)U^{C}(\mathbf{k})^{\dagger}\right]$
(52)

where $\lambda_{\mathbf{k},p}$ are eigenvalues of the matrix indexed by $p,q$
in the the brackets in Eq. (51). Thus the L1-norm of $H_{1}^{\prime}$ is the
sum

$\displaystyle\lambda_{\mathrm{DF},1}=\sum_{\mathbf{k}}\sum_{p}|\lambda_{\mathbf{k},p}|.$
(53)

(a)(b) Figure 6: (a) The number of $k$-points verses Toffoli cost to implement
the block encoding for the double factorization LCU evaluated for the
benchmark systems listed in Table 1 described using the cc-pVDZ and cc-pVTZ
basis sets and $\Gamma$-centered Monkhorst-Pack grids of size [1, 1, 1] to [3,
3, 3]. Each point is a single system described at a particular basis set and
$k$-mesh where the threshold to keep eigenvalues and vectors of the second
factorization is selected to produce two-electron integrals corresponding to
an MP2 error of one 1 milliHartree with respect to an untruncated double
factorization. On average this corresponds to a threshold value of $1\times
10^{-4}$ for the benchmark systems. The expected $\mathcal{O}(\sqrt{N_{k}})$
scaling improvement for symmetry-adapted walk operators is demonstrated. (b)
The value of $\lambda$ per unit cell as a function of the total system size
$NN_{k}$ for the same systems described with the same cutoffs used in (a). The
reduced variational freedom in compression of the two-electron integral
tensors for the symmetry-adapted walk operator construction translates to an
increased value of $\lambda$ at all system sizes.

Figure 6 demonstrates the improved $\sqrt{N_{k}}$ scaling of the block
encodings coming from reducing the number of controlled rotations by $N_{k}$.
Unlike the SF case $\lambda$ for DF has worse scaling in the symmetry-adapted
setting compared to the supercell case. This is rationalized by the fact that
there is a larger degree of variational freedom in the second factorization
for supercell calculations (and thus more compression) compared to the
symmetry-adapted case. The $\lambda$ value is basis set dependent and can
potentially be reduced by orbital optimization [79].

###  III.4 The tensor hypercontraction Hamiltonian representation

In the tensor hypercontraction (THC) LCU representation the fact that the two-
electron integrals can be represented in a symmetric Canonical Polyadic like
decomposition is used to define a set of non-orthogonal basis function in
which to represent the Hamiltonian, and we use a similar infrastructure to the
DF algorithm to implement each term in the factorization (which is in a
different non-orthogonal basis) sequentially. In the following section we
describe the Bloch orbital version (symmetry-adapted) of the THC decomposition
and the resulting LCU, $\lambda$ calculation, and qubitization complexities.
First we review the salient features of tensor hypercontraction for the
molecular case before introducing symmetry labels. Recall that in the
molecular THC approach we expand density like terms over a grid of $M$ points
(labeled $\mu$) and weight each grid point with a function $\xi_{\mu}(r)$

$\phi_{p}(\mathbf{r})\phi_{q}(\mathbf{r})\approx\sum_{\mu}\xi_{\mu}(\mathbf{r})\phi_{p}(\mathbf{r}_{\mu})\phi_{q}(\mathbf{r}_{\mu})$
(54)

which allows us to write the two-electron integral tensor as

$V_{pqrs}=\sum_{\mu\nu}\chi_{p}^{(\mu)}\chi_{q}^{(\mu)}\zeta_{\mu\nu}\chi_{r}^{(\nu)}\chi_{s}^{(\nu)}$
(55)

where the central tensor is defined as

$\zeta_{\mu\nu}=\int d\mathbf{r}\,\int
d\mathbf{r}^{\prime}\,\frac{\xi_{\mu}(\mathbf{r})\xi_{\nu}(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}.$
(56)

In order to incorporate translational symmetry into the THC factorization the
decomposition of the density is performed on the _cell periodic_ part of the
Bloch orbitals as [80, 81]

$u_{p\mathbf{k}_{p}}^{*}(\mathbf{r})u_{q\mathbf{k}_{q}}(\mathbf{r})\approx\sum_{\mu}\xi_{\mu}(\mathbf{r})u_{p\mathbf{k}_{p}}^{*}(\mathbf{r}_{\mu})u_{q\mathbf{k}_{q}}(\mathbf{r}_{\mu}),$
(57)

where
$u_{p\mathbf{k}_{p}}(\mathbf{r})=e^{-i\mathbf{k}_{p}\mathbf{r}}\phi_{p\mathbf{k}_{p}}(\mathbf{r})$.
Then the two-electron integral tensor has the form

$\displaystyle
V_{p\mathbf{k}_{p},q\mathbf{k}_{q},r\mathbf{k}_{r},s\mathbf{k}_{s}}$
$\displaystyle=\int d\mathbf{r}\int
d\mathbf{r}^{\prime}\phi_{p\mathbf{k}_{p}}^{*}\phi_{q\mathbf{k}_{q}}V(\mathbf{r},\mathbf{r}^{\prime})\phi_{r\mathbf{k}_{r}}^{*}\phi_{s\mathbf{k}_{s}}$
$\displaystyle=\sum_{\mu\nu}u_{p\mathbf{k}_{p}}^{*}(\mathbf{r}_{\mu})u_{q\mathbf{k}_{q}}(\mathbf{r}_{\mu})\zeta_{\mu\nu}^{\mathbf{k}_{p},\mathbf{k}_{q},\mathbf{k}_{r},\mathbf{k}_{s}}u_{r\mathbf{k}_{r}}^{*}(\mathbf{r}_{\nu})u_{s\mathbf{k}_{s}}(\mathbf{r}_{\nu})$
$\displaystyle=\sum_{\mu\nu}\chi_{p\mathbf{k}_{p}}^{(\mu)*}\chi_{q\mathbf{k}_{q}}^{(\mu)}\zeta_{\mu\nu}^{\mathbf{k}_{p},\mathbf{k}_{q},\mathbf{k}_{r},\mathbf{k}_{s}}\chi_{r\mathbf{k}_{r}}^{(\nu)*}\chi_{s\mathbf{k}_{s}}^{(\nu)}$
(58)

where $\chi_{q\mathbf{k}_{q}}^{(\mu)}=u_{q\mathbf{k}_{q}}(\mathbf{r}_{\mu})$,
$V(\mathbf{r},\mathbf{r}^{\prime})=|\mathbf{r}-\mathbf{r}^{\prime}|^{-1}$, and

$\displaystyle\zeta_{\mu\nu}^{\mathbf{k}_{p},\mathbf{k}_{q},\mathbf{k}_{r},\mathbf{k}_{s}}$
$\displaystyle=\int d\mathbf{r}\int
d\mathbf{r}^{\prime}e^{-i(\mathbf{k}_{p}-\mathbf{k}_{q})\cdot\mathbf{r}}\xi_{\mu}(\mathbf{r})V(\mathbf{r},\mathbf{r}^{\prime})\xi_{\nu}(\mathbf{r}^{\prime})e^{i(\mathbf{k}_{s}-\mathbf{k}_{r})\cdot\mathbf{r}^{\prime}}.$
(59)

Some care needs to be taken when bringing this into a form similar to Eq. 5.
First recall that we have
$\mathbf{k}_{p}-\mathbf{k}_{q}+\mathbf{k}_{r}-\mathbf{k}_{s}=\mathbf{G}_{pqrs}$,
where $\mathbf{G}_{pqrs}$ is a reciprocal lattice vector, and we are working
with a uniform $\Gamma$-point centered momentum grid with dimensions
$\mathbf{N}=[N_{x},N_{y},N_{z}]$ and $N_{k}=N_{x}N_{y}N_{z}$ . To eliminate
one of the four momentum modes, we identify
$\mathbf{Q}=\mathbf{k}_{p}{\ominus}\mathbf{k}_{q}$, and
$\mathbf{Q}=\mathbf{k}_{s}{\ominus}\mathbf{k}_{r}$, and set
$\mathbf{k}_{p}=\mathbf{k}$, $\mathbf{k}_{q}=\mathbf{k}{\ominus}\mathbf{Q}$,
$\mathbf{k}_{s}=\mathbf{k}^{\prime}$ and
$\mathbf{k}_{r}=(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})$. To evaluate the
$\zeta$ tensor we still need to know the values of
$\mathbf{k}_{p}-\mathbf{k}_{q}$ in absolute terms given a value for
$\mathbf{Q}$ and $\mathbf{k}$. We note that mapping the difference
$\mathbf{k}_{p}-\mathbf{k}_{q}$ back into our $k$-point mesh amounts to adding
a specific reciprocal lattice vector
$\mathbf{G}^{\mathbf{Q}}_{pq}=(\mathbf{k}_{p}-\mathbf{k}_{q})-\mathbf{Q}=(\mathbf{k}-(\mathbf{k}{\ominus}\mathbf{Q}))-\mathbf{Q}\equiv\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}$,
with a similar expression for $\mathbf{k}^{\prime}$ (the subtraction here is
not modular). Thus, given a $\mathbf{Q}$ and $\mathbf{k}$ we can determine
$\mathbf{k}-\mathbf{Q}$ and $\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}$.
With these replacements we can write

$\displaystyle
V_{p\mathbf{k}_{p},q\mathbf{k}_{q},r\mathbf{k}_{r},s\mathbf{k}_{s}}\rightarrow
V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}$
$\displaystyle=\sum_{\mu\nu}\chi_{p\mathbf{k}}^{(\mu)*}\chi_{q\mathbf{k}{\ominus}\mathbf{Q}}^{(\mu)}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})}^{(\nu)*}\chi_{s\mathbf{k}^{\prime}}^{(\nu)}$
$\displaystyle=\sum_{\mu\nu}\chi_{p\mathbf{k}}^{(\mu)*}\chi_{q\mathbf{k}{\ominus}\mathbf{Q}}^{(\mu)}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}},\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})}^{(\nu)*}\chi_{s\mathbf{k}^{\prime}}^{(\nu)},$
(60)

where we have used

$\displaystyle\zeta_{\mu\nu}^{\mathbf{k}_{p},\mathbf{k}_{q},\mathbf{k}_{r},\mathbf{k}_{s}}\rightarrow\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}$
$\displaystyle=\int d\mathbf{r}\int
d\mathbf{r}^{\prime}e^{-i(\mathbf{Q}+\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}})\cdot\mathbf{r}}\xi_{\mu}(\mathbf{r})V(\mathbf{r},\mathbf{r}^{\prime})\xi_{\nu}(\mathbf{r}^{\prime})e^{i(\mathbf{Q}+\mathbf{G}_{\mathbf{k}^{\prime}\mathbf{k}^{\prime}-\mathbf{Q}})\cdot\mathbf{r}^{\prime})}$
$\displaystyle=\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}},\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}}.$
(61)

In practice there are at most 8 values of $\mathbf{G}$, so we only need to
classically determine at most $8^{2}N_{k}$ values of $\zeta$, as opposed to
$N_{k}^{3}$.

We can then write

$\displaystyle H_{2}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}\sum_{pqrs}\sum_{\sigma\tau}V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{Q},\mathbf{k},\mathbf{k}^{\prime}}\sum_{pqrs}\sum_{\sigma\tau}\sum_{\mu\nu}\chi_{p\mathbf{k},\mu}^{*}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\mu}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}},\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\nu}^{*}\chi_{s\mathbf{k}^{\prime},\nu}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}\sum_{\mu\nu}\sum_{\sigma\tau}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}$
$\displaystyle\quad\times\left(\sum_{\mathbf{k}|\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}=\mathbf{G}_{1}}\sum_{pq}\chi_{p\mathbf{k},\mu}^{*}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\mu}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\left(\sum_{\mathbf{k}^{\prime}|\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}=\mathbf{G}_{2}}\sum_{rs}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\nu}^{*}\chi_{s\mathbf{k}^{\prime},\nu}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}\right),$
(62)

where in going from the second to the third line of Eq. 62 we have rewritten
the sum over $\mathbf{k}$ and $\mathbf{k}^{\prime}$ as a double sum over all
$8^{2}$ values of $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$, and a restricted sum
on $\mathbf{k}$ such that for a given $\mathbf{G}_{1}$ and $\mathbf{Q}$ we
only sum over those $\mathbf{k}$ which satisfy
$\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}=\mathbf{G}_{1}$. Here the
notation $\mathbf{G}_{\mathbf{k}_{p},\mathbf{k}_{q}}$ is used as equivalent to
$\mathbf{G}_{pq}$ above. The fourfold symmetry of the two-electron integrals
carries over to analogous symmetries in $\zeta$, which are listed in Appendix
D.

We will then define $\tilde{\chi}$ which are individually normalized for each
$\mathbf{k}$ and $\mu$ so
$\sum_{p}\tilde{\chi}_{p\mathbf{k},\mu}^{*}\tilde{\chi}_{p\mathbf{k},\mu}=1$
and

$\displaystyle\mathcal{N}_{\mathbf{k},\mu}\tilde{\chi}_{p\mathbf{k},\mu}=\chi_{p\mathbf{k},\mu}$
(63)

with
$\mathcal{N}_{\mathbf{k},\mu}:=\sqrt{\sum_{p}|\chi_{p\mathbf{k},\mu}|^{2}}$.
We then use these normalized $\tilde{\chi}$ to give transformed annihilation
and creation operators

$\displaystyle
c_{\mu\mathbf{k}\sigma}=\sum_{p}\tilde{\chi}_{p\mathbf{k},\mu}a_{p\mathbf{k}\sigma},\qquad
c_{\mu\mathbf{k}\sigma}^{\dagger}=\sum_{p}\tilde{\chi}_{p\mathbf{k},\mu}^{*}a_{p\mathbf{k}\sigma}^{\dagger}\,.$
(64)

We can then write the two-body Hamiltonian as

$\displaystyle\hat{H}_{2}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}\sum_{\mu,\nu}\sum_{\sigma\tau}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}$
$\displaystyle\quad\times\left(\sum_{\mathbf{k}|\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}=\mathbf{G}_{1}}\mathcal{N}_{\mathbf{k},\mu}\mathcal{N}_{\mathbf{k}{\ominus}\mathbf{Q},\mu}c_{\mu\mathbf{k}\sigma}^{\dagger}c_{\mu(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\left(\sum_{\mathbf{k}^{\prime}|\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}=\mathbf{G}_{2}}\mathcal{N}_{\mathbf{k}^{\prime}{\ominus}\mathbf{Q},\nu}\mathcal{N}_{\mathbf{k}^{\prime},\nu}c_{\nu(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}c_{\nu\mathbf{k}^{\prime}\tau}\right).$
(65)

A complication for the implementation is that we would like to be able to
choose the relative weighting between $\zeta$ and $\chi$ such that

$\sum_{\mathbf{k}|\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}=\mathbf{G}}\mathcal{N}_{\mathbf{k},\mu}\mathcal{N}_{\mathbf{k}{\ominus}\mathbf{Q},\mu}=1.$
(66)

The difficulty here is that the values of $\mathcal{N}_{\mathbf{k},\mu}$ only
depend on $\mathbf{k},\mu$, because they are based on
$\chi_{p\mathbf{k},\mu}$. This sum is also dependent on $\mathbf{Q}$ and
$\mathbf{G}$, so for this normalization condition to hold it would mean we
need to have $\chi_{p\mathbf{k},\mu}$ also dependent on $\mathbf{Q}$ and
$\mathbf{G}$ in a multiplicative factor (so a non-$\mu$-dependent way). That
will leave the normalized $\tilde{\chi}_{p\mathbf{k},\mu}$ unaffected, but
means that the values of $\mathcal{N}_{\mathbf{k},\mu}$ need to have
dependence on $\mathbf{Q},\mathbf{G}$, which will need to be taken account of
in the state preparation.

The form in Eq. (III.4) then gives us a recipe for block encoding the
Hamiltonian as a linear combination of unitaries.

  1. 1.

First prepare a superposition state proportional to

$\sum_{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2},\mu,\nu}\sqrt{|\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}|}\mathinner{|{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2},\mu,\nu}\rangle}.$
(67)

This state may be prepared via the coherent alias sampling approach with a
complexity dominated by the complexity of the QROM. Accounting for symmetry
the dimension is about $32N_{k}M^{2}$ and the size of the QROM output is
approximately the log of that plus the number of bits for the keep
probability. That gives a Toffoli complexity scaling as

$\sqrt{32N_{k}M^{2}[\log(32N_{k}M^{2})+\aleph]}\,.$ (68)

  2. 2.

For each of the two expressions in brackets in Eq. (III.4), a preparation over
$\mathbf{k}$ or $\mathbf{k}^{\prime}$ is needed to give a state of the form

$\sum_{\mathbf{k}|\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}=\mathbf{G}}\sqrt{\mathcal{N}_{\mathbf{k},\mu}\mathcal{N}_{\mathbf{k}{\ominus}\mathbf{Q},\mu}}\mathinner{|{\mathbf{k}}\rangle}.$
(69)

As explained above, the values of $\mathcal{N}_{\mathbf{k},\mu}$ need to be
chosen with (implicit) dependence on $\mathbf{Q},\mathbf{G}$ for this to be a
normalised state. This means that the amplitudes here need to be indexed by
$\mathbf{k}$, $\mathbf{Q}$, $\mathbf{G}_{1}$ and $\mu$. The restricted range
of values in the sum over $\mathbf{k}$ means that the indexing over
$\mathbf{k},\mathbf{Q},\mathbf{G}_{1}$ gives $N_{k}^{2}$ items of data, which
is multiplied by $M$ for the indexing over $\mu$. So there are $N_{k}^{2}M$
items of data needed, which is smaller than that in the first step, because it
is missing the factor of 32 and typically $N_{k}<M$. Given that the output
size is approximately $\log(N_{k})+\aleph$, the Toffoli complexity is
approximately

$\sqrt{N_{k}^{2}M[\log(N_{k})+\aleph]}.$ (70)

This cost is incurred twice, once for each of the factors in brackets in Eq.
(III.4).

  3. 3.

For each of the $c$ annihilation and creation operators we perform a rotation
of the basis from $a$. This is done in the following way.

    1. (a)

First use the spin $\sigma$ or $\tau$ to control a swap of the system
registers. This is done once and inverted for each of the two $c^{\dagger}c$
factors. Each of these 4 swaps has cost $N_{k}N/2$.

    2. (b)

Then use $\mathbf{k}$ or $\mathbf{k}{\ominus}\mathbf{Q}$ to control the swap
of the registers we wish to act on into working registers. The value of
$\mathbf{k}{\ominus}\mathbf{Q}$ is used for
$c_{\mu(\mathbf{k}{\ominus}\mathbf{Q})\sigma}$, and needs to be computed to
use as a control. Each of these eight swaps may be done with a Toffoli
complexity approximately as half the number of system registers $N_{k}N/2$.

    3. (c)

Next $\mathbf{k}$ (or $\mathbf{k}{\ominus}\mathbf{Q}$) and $\mu$ (or $\nu$)
are used as a control for a QROM to output the angles for Givens rotations.
There are two angles for each of $N/2$ Givens rotations, so if they have
$\beth$ each the size of the output is $N\beth$. Then the QROM complexity is
about

$\sqrt{N_{k}N^{2}\beth}.$ (71)

This must be done 4 times (and has a smaller erasure cost).

    4. (d)

The sequence of $N/2$ Givens rotations is performed, each with 4 individual
rotations on $\beth$, for a cost of $2N\beth$. This cost is incurred 8 times,
twice for each of the annihilation and creation operators.

  4. 4.

After the rotation of the basis, we simply need to perform the linear
combination of $\vec{Z}X$ and $\vec{Z}Y$ for $c^{\dagger}$ and $c$. The $X$ or
$Y$ is applied in a fixed location, but the $\vec{Z}$ needs to be applied on a
range of qubits chosen by $\mathbf{k}$ or $\mathbf{k}{\ominus}\mathbf{Q}$. We
therefore have approximately $N_{k}$ for the unary iteration for each
$\vec{Z}$ for a total cost of about $4N_{k}$.

Lastly we would perform reflections on control ancillas as usual to construct
a qubitised quantum walk from the block encoding. This cost is trivial
compared to that in the other steps. For a more detailed explanation, see the
circuit diagram in Figure 7 and the discussion in Appendix D.

Figure 7: The quantum circuit for the block encoding of the THC
representation, split into two parts with the right half at the bottom. The
top shows the portion of the circuit for the first part controlled by
$\mathbf{k}^{\prime}$ and $\nu$, and the bottom shows the (right) part of the
circuit where it is controlled by $\mathbf{k}$ and $\mu$. The dotted
rectangles show the regions for implementing the $c$ and $c^{\dagger}$
operators together with the Givens rotations needed to change the basis. The
swaps controlled by the $\mathbf{k}^{\prime}$ and $\mathbf{k}$ registers are
to move the appropriate qubits into target registers in order to apply the
Givens rotations. The $c$ and $c^{\dagger}$ are applied using a superposition
of $X$ and $iY$ applied using an ancilla qubit (not shown for simplicity),
together with a string of $Z$ gates for the Jordan-Wigner representation. The
preparation at the beginning includes an inequality test between $\mu$ and
$\nu$ to give a qubit flagging whether the real or imaginary part is produced.
To make the implementation self-inverse, the $\mu,\nu$ and
$\mathbf{k},\mathbf{k}^{\prime}$ pairs of registers are swapped in the middle
(the left of the lower half). Also, an $X$ gate is applied to the qubit that
controls the swaps at the beginning and end. Figure 8: Violin plot of absolute
errors in the $k$-THC-MP2 energy per cell for the benchmark set in Table 1.
Here we compare the MP2 errors as a function of the THC rank parameter
$c_{\mathrm{THC}}$ using ISDF or subsequent reoptimization to generate the THC
factors. (a)(b) Figure 9: (a) The number of $k$-points verses Toffoli cost to
implement the block encoding for the THC factorization LCU evaluated for the
benchmark systems listed in Table I described using the cc-pVDZ and cc-pVTZ
basis sets and $\Gamma$-centered Monkhorst- Pack grids of size [1, 1, 1] to
[3, 3, 3]. Each point is a single system described at a particular basis set
and k-mesh where the range of the auxiliary index of the THC factorization is
selected to produce two-electron integrals corresponding to an MP2 error of
one 1 milliHartree with respect to an untruncated auxiliary index range. This
corresponds to an auxiliary index that is eight times the number of orbitals
in the primitive cell for symmetry adapted THC, and eight times the total
number of orbitals ($N_{k}N$) for supercell THC. We divide the Toffoli
complexity for implementing SELECT + PREPARE + PREPARE-1 by N, which is the
shared scaling in the number of bands. While we observe a $\sqrt{N_{k}}$
scaling improvement for symmetry-adapted walk operations, we believe this is a
finite size effect and both methods should scale linear with $N_{k}$ for
sufficiently large $N_{k}$ The value of $\lambda$ as a function of the total
system size $NN_{k}$ for the same systems described with the same cutoffs used
in (a). The reduced variational freedom in compression of the two-electron
integral tensors for the symmetry-adapted walk operator construction
translates to an increased value of $\lambda$ at all system sizes. Figure 10:
Synthetic data for the number of Toffolis required to implement the
qubitization oracles with the $k$-THC factorization demonstrating the
challenge of extracting the correct asymptotic scaling with limited finite
size data. To generate the data we used the system parameters of carbon
diamond in the cc-pVDZ basis set ($N=52$, $M=208$)

The $\lambda_{\rm THC}$ value has a one-body component and two-body component.
Unlike molecular THC where the two-body component is reduced because we evolve
by number operators in the non-orthogonal basis, in this version of the THC
algorithm we will evolve by ladder operators in a non-orthogonal basis, and
thus there is no one-body part to remove. The one-body contribution to
$\lambda_{\rm THC}$, $\lambda_{{\rm THC},1}$, is computed in a similar way as
for the double factorization algorithm but noting that the extra factor of
$1/2$ coming from the $Z$ operator is no-longer present to cancel the factor
of two from spin summing. The one-body contribution to $\lambda_{\rm THC}$ is

$\displaystyle\lambda_{{\rm THC},1}$
$\displaystyle=2\sum_{\mathbf{k}}\sum_{p}|\lambda_{p,\mathbf{k}}|.$ (72)

The two-body contribution to $\lambda_{\rm THC}$, $\lambda_{{\rm THC},2}$, is
determined by summing over all unitaries in the LCU. This summation can be
rewritten in the form

$\displaystyle\lambda_{{\rm THC},2}$
$\displaystyle=2\sum_{\mathbf{Q}}\sum_{\mu,\nu}\sum_{\mathbf{G}_{1},\mathbf{G}_{2}}\left(|\mathrm{Re}[\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}]|+|\mathrm{Im}[\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}]|\right)$
$\displaystyle\quad\times\left(\sum_{\mathbf{k}|\mathbf{G}_{\mathbf{k},\mathbf{k}{\ominus}\mathbf{Q}}=\mathbf{G}_{1}}\mathcal{N}_{\mathbf{k},\mu}\mathcal{N}_{\mathbf{k}{\ominus}\mathbf{Q},\mu}\right)\left(\sum_{\mathbf{k}^{\prime}|\mathbf{G}_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}=\mathbf{G}_{2}}\mathcal{N}_{\mathbf{k}^{\prime}{\ominus}\mathbf{Q},\nu}\mathcal{N}_{\mathbf{k}^{\prime},\nu}\right)$
(73)

using the expression for $\zeta$ described in Eq. (III.4).

To obtain resource estimates for THC with $k$-points we follow a similar
procedure to previous molecular work [32] and first compress the rank of the
THC factors ($M=c_{\mathrm{THC}}N/2$, where $c_{\mathrm{THC}}$ is the THC rank
parameter). In particular, we use the interpolative separable density fitting
(ISDF) approach [80, 82, 83] as a starting point before subsequently
reoptimizing these factors in order to compress the THC rank while
regularizing $\lambda$ [32, 37], which we will call $k$-THC. Further details
of this procedure are provided in Appendix G. In Fig. 8 we demonstrate that a
$c_{\mathrm{THC}}=8$ is sufficient to obtain MP2 correlation energies within
approximately $0.1$ mHa/Cell for a subset of the systems considered in the
benchmark set. We note that the equivalent ISDF rank may be on the order of
10-15 for comparable accuracy, which would correspond to a much larger value
for $\lambda$.

Fig. 9 (a) demonstrates a $\sqrt{N_{k}}$ scaling improvement of the block
encodings in the symmetry adapted case. Note that this $\sqrt{N_{k}}$ speedup
for the block encodings is partially a finite size effect. In Fig. 10 we plot
the Toffoli complexity per step as a function of $N_{k}$ using artificially
generated data to explore the large $N_{k}$ behavior. We see that depending on
the fitting range employed the extracted asymptotic scaling trends towards
linear. While ultimately both the symmetry-adapted and supercell encodings
should scale linearly with the system size due to the cost of unary iteration
over all basis states, there are several factors that yield a $\sqrt{N_{k}}$
saving in the symmetry-adapted case, and the relative size of the prefactors
becomes important. Similar to DF, we find from Fig. 9 (b) that $\lambda$ in
the symmetry-adapted setting exhibits slightly worse scaling than for
supercell calculations. This worsening of $\lambda$ in the symmetry-adapted
case can be understood again as a reduction in variational freedom in the
symmetry adapated case, leading to smaller compression. Note that while Eq. 73
nominally scales cubicly with $N_{k}$, we expect each individual matrix
element to decay like $N_{k}^{-1}$, which yields the expected quadratic
dependence of $\lambda$, or a linear dependence of $\lambda$ when targeting
the total energy per cell. In the supercell case, there are simply
$M^{2}=(N_{k}N)^{2}$ elements in the central tensor, which in turn controls
the scaling of $\lambda$. From Table 2 and Fig. 9 we can conclude that there
is asymptotically no advantage to incorporating symmetry in the THC
factorization for the Toffoli complexity, with both the supercell and
symmetry-adapted methods exhibiting approximately quadratic scaling with
system size for a fixed target accuracy of the total energy per cell.

##  IV Scaling comparison and runtimes for diamond

We now compare runtimes and estimate total physical requirements to simulate
Diamond as a representative material. In Figure 11 we plot the total Toffoli
complexity for the sparse, SF, DF, and THC LCUs using symmetry-adapted block
encodings and supercell calculations for Diamond with cc-pVDZ and cc-pVTZ
basis sets at various Monkhorst-Pack samplings. In sparse and SF there is a
clear asymptotic separation between supercell and symmetry-adapted Toffoli
counts. This is expected from the fact that both block encoding constructions
are asymptotically improved and $\lambda$ does not increase. For the DF case,
total Toffoli complexity for supercell and symmetry-adapted cases is similar
due to the larger $\lambda$ for the symmetry-adapted algorithm. For THC, the
total Toffoli complexity is similar in the supercell and symmetry adapted
case, but the asymptotic scaling is identical for the supercell and symmetry-
adapted algorithms. This is due to the increase in $\lambda$ for the symmetry-
adapted algorithm.

In Table 3 we tabulate the quantum resource requirements and estimated
runtimes after compiling into a surface code using physical qubits with error
rates of $0.01$% and a $1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ cycle time.
We assume four Toffoli factories similar to References [32] and [37] and
observe that for systems with 52-1404 spin-orbitals the quantum resource
estimates are roughly in line with extrapolated estimates from the molecular
algorithms.

It is important to note that while the THC resource requirements look
competitive for these small systems, in its current form it is not a practical
way to simulate materials at scale. This is due to the prohibitive cost of
reoptimizing the THC factors which significantly limits the system sizes that
can be simulated. Moreover, as discussed in Section III.4, we caution that the
THC trend lines are only valid within the fitting range, and we expect that
asymptotic THC Toffoli count will trend more towards $\mathcal{O}(N_{k}^{2})$
in the thermodynamic limit.

(a)(b) Figure 11: (a) Total Toffoli requirements for Diamond in a cc-pVDZ basis at various Monkhorst-Pack samplings of the Brillouin zone with $\Gamma$-point centered grids of size [1,1,1] to [3, 3, 3]. Dashed lines are fits to the supercell data that is not plotted. Solid lines are fits to the symmetry-adapted data shown as data points. (b) Total logical qubits for symmetry-adapted oracles and supercell (dotted lines). All values are estimated from 0.1 mHa per unit cell thresholds on the MP2 energy. In the case of THC we only plot the symmetry adapted data due limited THC data arising from difficulty in optimizing the supercell THC factors.  LCU |  $k$-mesh | Toffolis | Logical Qubits | Physical Qubits[M] | Surface Code Runtime [days]
---|---|---|---|---|---
sparse | $[1,1,1]$ | $4.84\times 10^{9}$ | 2478 | $2.20$ | $9.10\times 10^{-1}$
| $[2,2,2]$ | $2.66\times 10^{12}$ | 75287 | $90.57$ | $5.77\times 10^{2}$
| $[3,3,3]$ | $1.06\times 10^{14}$ | 374274 | $543.76$ | $2.61\times 10^{4}$
SF | $[1,1,1]$ | $3.20\times 10^{9}$ | 2283 | $2.05$ | $6.02\times 10^{-1}$
| $[2,2,2]$ | $3.27\times 10^{12}$ | 20567 | $24.91$ | $7.11\times 10^{2}$
| $[3,3,3]$ | $1.13\times 10^{15}$ | 47665 | $69.52$ | $3.10\times 10^{5}$
DF | $[1,1,1]$ | $9.61\times 10^{8}$ | 2396 | $1.55$ | $1.81\times 10^{-1}$
| $[2,2,2]$ | $6.74\times 10^{10}$ | 18693 | $18.47$ | $1.27\times 10^{1}$
| $[3,3,3]$ | $1.09\times 10^{12}$ | 68470 | $82.39$ | $2.37\times 10^{2}$
THC | $[1,1,1]$ | $1.67\times 10^{10}$ | 18095 | $14.20$ | $3.14$
| $[2,2,2]$ | $4.85\times 10^{11}$ | 36393 | $35.60$ | $1.05\times 10^{2}$
Table 3: Diamond represented in a cc-pVDZ basis (52 spin-orbitals in the
primitive cell) at various $k$-mesh sizes and the associated quantum resource
requirements to compute the total energy per cell to within 1 kcal / mol. The
surface code runtime is estimated using four T-factories, a physical error
rate of $0.01$%, and a cycle time of $1\text{\,}\mathrm{\SIUnitSymbolMicro
s}$. The physical qubit count is given in millions.

##  V Classical and quantum simulations of LNO

In this section, we compare modern classical computational methods with
quantum resource estimates in the context of a challenging problem of
industrial interest: the ground state of LiNiO2.

###  V.1 LNO background

Layered oxides have been the most popular cathode active materials for Li-ion
batteries since their commercialization in the early ‘90s. While LiCoO2 is
still the material of choice in the electronics industry, the increasing
human, environmental and financial cost of cobalt spells out the need for
cobalt-free cathode active materials, especially for automotive
applications[84, 85].

The isostructural compound LiNiO2 (LNO) had been identified as an ideal
replacement for LiCoO2 already in the ‘90s, due to its comparably high
theoretical capacity at a lower cost [86, 87]. Despite its numerous drawbacks,
LNO still serves as the perfect model system for many derivative compounds
such as lithium nickel-cobalt-manganese (NCM) and lithium nickel-cobalt-
aluminum oxides (NCA) that are nowadays the gold standard in the automotive
industry [38]. Moreover, the constant demand for better performing materials
pushes the amount of substituted Ni to the dilute regime and the research
trend is approaching the asymptotic LiNiO2 limit, making LiNiO2 a system of
interest in battery research [38].

Even the nature of the ground state of LNO is still under debate. The
universally observed rhombohedral $\mathrm{R\bar{3}m}$ symmetry [38], with Ni
being octahedrally coordinated to six oxygen atoms through six equivalent Ni-O
bonds conflicts with the renowned Jahn-Teller (JT) activity of low-spin
trivalent Ni, which has been experimentally proven on a local scale [38]. In a
recent DFT study [39], we argued that this apparent discrepancy might be
resolved by the dynamics and low spatial correlation of Jahn-Teller
distortions. In that work, the energy distance between Jahn-Teller distorted
and non-distorted candidates (Figure 12) compared to zero-point vibrational
energies makes a strong argument in favor of the dynamic Jahn-Teller effect. A
non-JT distorted structure resulting from the disproportionation of Ni3+ has
also been reported as a ground state candidate [40] despite the 1:1 ratio
between long and short Ni-O bonds, which conflicts with the experimentally
determined 2:1 ratio. In the original study, the stability of this structure
has been found to depend heavily on the value of the on-site Hubbard
correction applied to the PBE functional. With the SCAN-rVV10 functional (with
and without on-site Hubbard correction) [39], this candidate is consistently
less stable than the JT-distorted models; it is also worth mentioning that the
on-site Hubbard correction considerably increases the stability of the JT-
distorted models. The dependence of Jahn-Teller stabilization energies on the
functional had already been observed by Radin [88] and is ascribed to the
difficulty to adequately describe the doubly degenerate high-symmetry,
undistorted state.

In light of previous studies, we will focus on four candidate structures for
the LNO ground state. These structures are shown in Figure 12. We will
furthermore focus only on the energetics of the problem. The goal is to
compute the relative energies of these different crystal structures without
the uncertainty of DFT.

Figure 12:  The four known LiNiO2 polymorphs: high-symmetry
$\mathrm{R\bar{3m}}$, collinear JT-distorted $\mathrm{C2/m}$, zig-zag JT-
distorted $\mathrm{P2{{}_{1}}/c}$, and disproportionated $\mathrm{P2/c}$.
Green spheres represent Li, gray polyhedra are NiO6 octahedra, and elongated
Ni-O bonds are depicted as bold blue arrows.

###  V.2 Correlated $k$-point calculations

Local-basis quantum chemistry methods for electron correlation have been
increasingly applied to periodic systems as an alternative to DFT with more
controllable accuracy. Here we apply two such methods, second order Møller-
Plesset perturbation theory (MP2)[89, 90] and coupled cluster singles and
doubles (CCSD)[91, 92], to the three distorted structures of LNO (Figure 12).
Local basis methods like these can be directly compared to quantum algorithms
described in this work, since both are formulated within the same framework of
a crystalline Gaussian one-particle basis. While these methods cannot be
easily applied to the symmetric structure, which is metallic at the mean-field
level, they should provide accurate results for the distorted structures
provided that the finite-size and finite-basis errors can be controlled. All
mean field, MP2 and CCSD calculations were performed with the PySCF program
package [93, 94]. QMCPACK [95, 96] was used to perform the ph-AFQMC
calculations, where we used at least 600 walkers and a timestep of 0.005 Ha-1.
The population control bias was found to be negligible. In all calculations,
we use separable, norm-conserving GTH pseudopotentials [97, 65] that have been
recently optimized for Hartree-Fock [98]. In all calculations on LNO we use
the GTH basis sets [41, 99] (GTH-SZV and GTH-DZVP specifically) that are
distributed with the CP2K [42] and PySCF [94] packages. In Figure 13 we show
the convergence of the minimal-basis MP2 energy as a function of effective
cell size for increasingly large $k$-point calculations. This demonstrates the
essential difficulty in converging to the bulk limit for correlated
calculations: the finite-size error will converge with $n_{k}^{-1/3}$.
Shifting the $k$-point grid to (1/8, 1/8, 1/8) and/or twist averaging (TA)
does not change the asymptotic behavior of the energy. In all other LNO
calculations, we use $\Gamma$-centered $k$-point grids. In all calculations,
the density of $k$-points along each reciprocal lattice vector was chosen so
that the density of $k$-points is as close to constant as possible.

Figure 13:  Convergence of the total MP2 energy for the C2/m structure using
the minimal basis set as a function of $n_{k}^{-1/3}$ where $n_{k}$ is the
number of $k$-points. $\Gamma$-centered and shifted $k$-point grids are
compared with and without twist averaging (TA).

While a minimal basis is useful for a qualitative understanding of the finite-
size error, it is does not provide sufficient accuracy to resolve the
different LNO structures examined in this work. The double-zeta basis set
(GTH-DZVP) is large enough to provide qualitative accuracy, but converging the
result to the bulk limit is prohibitively expensive for the systems considered
here. We can nonetheless provide some estimates of the ground-state CCSD and
MP2 (DZVP) energies as shown in Table 4 and 5. Since there is no evidence of
particularly ”strong correlation” in any of these systems (see Appendix E for
a more detailed discussion), MP2 and CCSD should provide qualitatively correct
estimates of the ground state energy. The unusually large MP2 correlation
energy for the P2/c structure suggests it may not be as reliable for this
structure, and this suspicion is confirmed by the CCSD and ph-AFQMC
calculations. For CCSD and ph-AFQMC, the P21/c structure is lowest in energy
which agrees qualitatively with the DFT calculations in Ref. [39]. However,
this prediction carries with it a great deal of uncertainty due to the small
simulation size, small one-particle basis set, and error in the MP2/CCSD/ph-
AFQMC approximations.

structure |  $k$-points | ROHF | MP2 | CCSD | ph-AFQMC
---|---|---|---|---|---
C2/m | 2x2x1 | -206.557491 | -0.750524 | -0.767350 | -0.7997(5)
P21/c | 1x2x1 | -206.567049 | -0.747717 | -0.765445 | -0.7982(5)
P2/c | 1x1x1 | -206.551767 | -0.811386 | -0.780580 | -0.8078(4)
Table 4:  ROHF total energy and MP2, CCSD, and ph-AFQMC correlation energies computed with a double-zeta basis set and a small k-mesh, the equivalent of 4 primitive formula units (16 atoms total) for each distorted structure. All units are Hartrees per formula unit. structure | ROHF | MP2 | CCSD | ph-AFQMC
---|---|---|---|---
C2/m | 260 | 184 | 208 | 218(18)
P2/c | 416 | -1317 | 4 | 155(17)
Table 5:  Energies in meV relative to the energy of the P21/c structure for
each method. The P21/c structure is lowest in energy for all methods except
MP2.

###  V.3 Single shot density matrix embedding theory

Another way to apply high-level correlated methods to periodic solids is
through quantum embedding methods in which a local impurity is treated with a
high-level method and the remainder of the system, the bath, is treated at a
lower level of theory. For periodic solids, dynamical mean-field theory (DMFT)
is perhaps the most widely successful such method [100, 101, 102, 103, 20].
Density matrix embedding theory (DMET) is an efficient quantum embedding
method for the ground state of quantum systems [104, 105], and it has recently
been applied to periodic solids with a fully ab initio Hamiltonian [20].

Though very large impurities are necessary to converge to the bulk limit of
the correlated method used for the impurity, a fixed impurity size provides a
local, systematically improvable approximation to the correlation energy. This
is particularly useful in cases where a local treatment of correlation is
sufficient for a qualitatively correct solution. Here we apply DMET with a
CCSD impurity solver to the distorted structures in a minimal basis set (GTH-
SZV). The libdmet code [19] was used for the DMET calculations with the PySCF
program package [93, 94] used in the impurity solver.

Figure 14 shows the minimal-basis DMET results with for each of the three
distorted structures.

Figure 14:  Convergence of the total DMET energy for a 4 formula unit (16
atom) impurity with respect to effective size of the mean-field calculation
for the different distorted structures.

In the context of DMET, we can effectively converge the mean-field part of the
problem. Unfortunately, the small one-particle basis set and modest impurity
size make it unable to meaningfully resolve these three structure of LiNiO2.
Quantum simulation can potentially overcome some of these limits by acting as
a lower scaling, unbiased impurity solver [106].

###  V.4 Quantum resource estimates for LNO

Quantum resource estimates for LNO using the SF and DF LCUs are reported in
Table 6. THC is not reported due to the difficulty of re-optimizing the THC
tensors to have low L1-norm as discussed in References [32, 37, 79]. For the
sparse LCU, a threshold of $1\times 10^{-4}$ was determined by averaging the
thresholds for the systems in Table 1 required to achieve 1 m$E_{\mathrm{H}}$
per unit cell. For the SF LCU the truncation of the auxiliary basis was set to
eight times the number of molecular orbitals which was determined by requiring
the error in the MP2 energy for the smallest C2/m system to be less than 1
m$E_{\mathrm{h}}$ per formula unit. For DF, the same requirement was used to
determine a cutoff for the second factorization of $1\times 10^{-3}$. The
trends are consistent with what was observed in Section IV: DF is consistently
more efficient than either sparse or SF LCUs.

For the smaller systems, these calculations are anticipated to be useful for
benchmarking faster classical methods. For the larger systems, the estimated
run times are daunting, but we are optimistic that further algorithmic
improvements can make calculations like these feasible in the future.

System | LCU |  $k$-mesh | $\lambda$ | Num. Spin-Orbs. | Toffolis | Logical Qubits | Physical Qubits [M] | run time [days]
---|---|---|---|---|---|---|---|---
$\mathrm{R\bar{3}m}$ | Sparse | $[2,2,2]$ | $120382.037$ | 116 | 6.16$\times 10^{13}$ | 166946 | $242.72$ | 1.51$\times 10^{4}$
|  | $[3,3,3]$ | $718377.133$ | 116 | 3.57$\times 10^{15}$ | 1625295 | $2808.82$ | 9.82$\times 10^{5}$
| SF | $[2,2,2]$ | $183778.821$ | 116 | 7.86$\times 10^{13}$ | 89162 | $129.77$ | 1.93$\times 10^{4}$
|  | $[3,3,3]$ | $2966279.293$ | 116 | 4.60$\times 10^{15}$ | 404723 | $699.68$ | 1.27$\times 10^{6}$
| DF | $[2,2,2]$ | $10730.422$ | 116 | 4.97$\times 10^{12}$ | 149939 | $180.16$ | 1.08$\times 10^{3}$
|  | $[3,3,3]$ | $44794.803$ | 116 | 7.28$\times 10^{13}$ | 598286 | $869.02$ | 1.79$\times 10^{4}$
$\mathrm{C2/m}$ | Sparse | $[2,2,1]$ | $58422.522$ | 116 | 1.03$\times 10^{13}$ | 83532 | $100.47$ | 2.53$\times 10^{3}$
|  | $[4,4,2]$ | $893339.394$ | 116 | 5.37$\times 10^{15}$ | 3051285 | $5272.93$ | 1.48$\times 10^{6}$
| SF | $[2,2,1]$ | $95803.204$ | 116 | 2.05$\times 10^{13}$ | 44657 | $53.90$ | 5.05$\times 10^{3}$
|  | $[4,4,2]$ | $2899609.300$ | 116 | 5.23$\times 10^{15}$ | 405310 | $700.69$ | 1.44$\times 10^{6}$
| DF | $[2,2,1]$ | $4873.648$ | 116 | 1.18$\times 10^{12}$ | 75178 | $90.44$ | 2.56$\times 10^{2}$
|  | $[4,4,2]$ | $51416.281$ | 116 | 9.82$\times 10^{13}$ | 598736 | $869.68$ | 2.41$\times 10^{4}$
$\mathrm{P2/c}$ | Sparse | $[1,1,1]$ | $84977.359$ | 464 | 2.06$\times 10^{13}$ | 99918 | $120.21$ | 5.07$\times 10^{3}$
|  | $[2,2,2]$ | $1627121.892$ | 464 | 1.67$\times 10^{16}$ | 3182362 | $6454.14$ | 4.59$\times 10^{6}$
| SF | $[1,1,1]$ | $201894.726$ | 464 | 8.74$\times 10^{13}$ | 92786 | $135.04$ | 2.15$\times 10^{4}$
|  | $[2,2,2]$ | $5666363.179$ | 464 | 2.07$\times 10^{16}$ | 839487 | $1450.95$ | 5.68$\times 10^{6}$
| DF | $[1,1,1]$ | $2753.901$ | 464 | 9.72$\times 10^{11}$ | 75834 | $91.23$ | 2.11$\times 10^{2}$
|  | $[2,2,2]$ | $40788.113$ | 464 | 1.40$\times 10^{14}$ | 1192900 | $1732.40$ | 3.44$\times 10^{4}$
$\mathrm{P2_{1}/c}$ | Sparse | $[1,2,1]$ | $105584.297$ | 232 | 3.39$\times 10^{13}$ | 182864 | $265.83$ | 8.34$\times 10^{3}$
|  | $[2,4,2]$ | $1714723.913$ | 232 | 1.50$\times 10^{16}$ | 3116825 | $6321.24$ | 4.12$\times 10^{6}$
| SF | $[1,2,1]$ | $271178.934$ | 232 | 8.92$\times 10^{13}$ | 96882 | $140.98$ | 2.19$\times 10^{4}$
|  | $[2,4,2]$ | $7798992.981$ | 232 | 2.13$\times 10^{16}$ | 438080 | $757.32$ | 5.85$\times 10^{6}$
| DF | $[1,2,1]$ | $3958.111$ | 232 | 1.27$\times 10^{12}$ | 75383 | $90.69$ | 2.76$\times 10^{2}$
|  | $[2,4,2]$ | $46189.645$ | 232 | 1.23$\times 10^{14}$ | 1192758 | $1732.20$ | 3.02$\times 10^{4}$
Table 6: Quantum Resource estimates for all four LNO structures normalized by
the number of formula units represented in each simulation cell.
$\mathrm{R\bar{3}m}$ and $\mathrm{C2/m}$ are both one formula unit while
$\mathrm{P2/c}$ is four formula units and $\mathrm{P2_{1}/c}$ is two formula
units. The sparse threshold is selected to be $1.0\times 10^{-4}$, the SF the
auxiliary index is truncated at eight times the number of molecular orbitals,
and the DF the second factorization is truncated at $1.0\times 10^{-4}$.

##  VI Conclusion

In this work we developed the theory of symmetry-adapted block encodings for
extended system simulation using four different representations of the
Hamiltonian as LCUs in order to improve quantum resource costs for reaching
the thermodynamic limit when simulating solids. In order to realize an
asymptotic speedup due to symmetry, we substantially modify the block
encodings compared with their molecular counterparts. To demonstrate these
asymptotic improvements we compiled constant factors for all four LCUs and
compared their performance on a suite of benchmark systems and a realistic
problem in materials simulation. We find that despite a clear asymptotic
speedup for walk operator construction there are competing factors (such as
lower compression in Hamiltonian tensor factorizations) that make it difficult
to observe a large speedup using symmetry. It was recently shown that
variationally constructing tensor compressions for Hamiltonian simulation can
improve quantum resource requirements [107, 79, 37] and thus we believe the
compressions can be improved to ultimately demonstrate a speedup for these
types of simulations.

For the sparse and SF LCUs we derive a $\mathcal{O}(\sqrt{N_{k}})$ speedup in
constructing select and prepare by ensuring only the minimal amount of
symmetry unique information is accessed by the quantum circuit through QROM.
In both cases a speedup is observable, though it is much clearer in the SF
case. Observing the sparse LCU speedup is more challenging due to the
difficulty of converging the $N_{k}$ and $N$ dependence of the two-electron
integrals. Compared with the molecular case where sparse was competitive with
the DF and THC algorithms [71, 32], largely due to the simplicity of select,
we find that sparse is not viable for converging to the thermodynamic limit of
solids.

The DF and THC tensor factorizations yield LCUs as unitaries in non-orthogonal
bases and lead to much higher compression than sparse and SF LCUs. In the DF
case we derive an asymptotic $\mathcal{O}(\sqrt{N_{k}})$ improvement in