# Optimized Quantum Simulation Algorithms for Scalar Quantum Field Theories

Andrew Hardy<EMAIL_ADDRESS>Department of Physics, University
of Toronto, Toronto, ON M5S 1A7 , Canada  Priyanka Mukhopadhyay
<EMAIL_ADDRESS>Department of Computer Science, University
of Toronto, Toronto, ON, M5S 2E4, Canada  These authors made equivalent
contributions to this work.  M. Sohaib Alam<EMAIL_ADDRESS>Quantum Artificial
Intelligence Laboratory (QuAIL), NASA Ames Research Center, Moffett Field, CA,
94035, USA  USRA Research Institute for Advanced Computer Science (RIACS),
Mountain View, CA, 94043, USA  These authors made equivalent contributions to
this work.  Robert Konik<EMAIL_ADDRESS>Condensed Matter and Materials Science
Division, Brookhaven National Laboratory, Upton, NY 11973, USA  Layla Hormozi
Computational Science Initiative, Brookhaven National Laboratory, Upton, NY
11973, USA  Eleanor Rieffel<EMAIL_ADDRESS>Quantum Artificial
Intelligence Laboratory (QuAIL), NASA Ames Research Center, Moffett Field, CA,
94035, USA  Stuart Hadfield  Quantum Artificial Intelligence Laboratory
(QuAIL), NASA Ames Research Center, Moffett Field, CA, 94035, USA  USRA
Research Institute for Advanced Computer Science (RIACS), Mountain View, CA,
94043, USA  João Barata  Physics Department, Brookhaven National Laboratory,
Upton, NY 11973, USA  Raju Venugopalan  Physics Department, Brookhaven
National Laboratory, Upton, NY 11973, USA  Dmitri E. Kharzeev  Center for
Nuclear Theory, Department of Physics and Astronomy, Stony Brook University,
Stony Brook, NY 11794-3800, USA  Physics Department, Brookhaven National
Laboratory, Upton, NY 11973, USA  Nathan Wiebe<EMAIL_ADDRESS>Department of Computer Science, University of Toronto, Toronto, ON, M5S 2E4,
Canada  Pacific Northwest National Laboratory, Richland, WA, 99354, USA
Canadian Institute for Advanced Studies, Toronto, ON, M5G 1M1, Canada

###### Abstract

We provide practical simulation methods for scalar field theories on a quantum
computer that yield improved asymptotics as well as concrete gate estimates
for the simulation and physical qubit estimates using the surface code. We
achieve these improvements through two optimizations. First, we consider a
different approach for estimating the elements of the S-matrix. This approach
is appropriate in general for 1+1D and for certain low-energy elastic
collisions in higher dimensions. Second, we implement our approach using a
series of different fault-tolerant simulation algorithms for Hamiltonians
formulated both in the field occupation basis and field amplitude basis. Our
algorithms are based on either second-order Trotterization or qubitization.
The cost of Trotterization in occupation basis scales as
$\widetilde{O}(\lambda N^{7}|\Omega|^{3}/(M^{5/2}\epsilon^{3/2})$ where
$\lambda$ is the coupling strength, $N$ is the occupation cutoff $|\Omega|$ is
the volume of the spatial lattice, $M$ is the mass of the particles and
$\epsilon$ is the uncertainty in the energy calculation used for the
$S$-matrix determination. Qubitization in the field basis scales as
$\widetilde{O}(|\Omega|^{2}(k^{2}\Lambda+kM^{2})/\epsilon)$ where $k$ is the
cutoff in the field and $\Lambda$ is a scaled coupling constant. We find in
both cases that the bounds suggest physically meaningful simulations can be
performed using on the order of $4\times 10^{6}$ physical qubits and $10^{12}$
$T$-gates which corresponds to roughly one day on a superconducting quantum
computer with surface code and a cycle time of 100 ns, placing simulation of
scalar field theory within striking distance of the gate counts for the best
available chemistry simulation results.

††preprint: APS/123-QED

##  I Introduction

Quantum simulation has been one of the great success stories of quantum
computing [1, 2, 3, 4, 5]. It has led to the realization that non-relativistic
quantum mechanics can, under most physically reasonable assumptions, be
simulated in polynomial time on a quantum computer. Further, recent work in
chemistry has shown that physically meaningful problems can be simulated using
fewer than a million physical qubits and a few hours worth of compute time
[6]. In contrast, there have been numerous classical methods developed to
simulate quantum field theories (QFT), with the state of the art approaches
summarized in [7]. However, even sign-problem free approaches relying on
Hamiltonian representations [8, 9], require an exponential number amount of
resources [10]. This opens up the possibility that quantum computing may
provide the only means possible for us to numerically simulate some of the
most challenging problems in quantum field theory [11, 12, 13, 14, 15].
Despite these advances, we do not yet know whether all physically meaningful
quantum field theories can be efficiently simulated, nor are there at present
detailed estimates of the memory and number of quantum operations required to
perform quantum simulations of field theories.

The seminal work of Jordan, Lee and Preskill (JLP) [16, 17] provided a major
step towards addressing the question of whether quantum computers can
efficiently simulate quantum field theory algorithms that are capable of
simulating scattering events for scalar $\phi^{4}$ theory. In addition to
describing aspects of the dynamics and couplings of the Higgs boson in the
standard model, $\phi^{4}$ theory is an excellent model field field theory
that captures non-trivial features of the dynamics of the other fundamental
quantum fields in the standard model. Even in 1+1D, classical $\phi^{4}$
theory has non-trivial bound states described by kinks and anti-kinks. The
theory also contains stationary topological solitons that are periodic in
time, known as “breather”-modes [18], providing insight into the nonlinear
dynamics of domains walls in fields ranging from condensed matter to cosmology
[18]. Further, at the very high energies now accessible at colliders, as
exemplified in the parton model of Feynman and Bjorken [19, 20], scattering
processes involving fermions or bosons become alike. Thus, the study of the
$\phi^{4}$ scalar field theory, at high energies, provides some insight into
the dynamics of more phenomenologically rich theories, which are harder to
directly simulate.

The JLP [16, 17] proposal focused on computing the elements of the Scattering
Matrix, or $S$-matrix, that describes the entire information content of a
scattering process in quantum mechanics and, more generally, in quantum field
theory. Specifically, the $S$-matrix describes how free fields of a theory at
asymptotic time $T=-\infty$ evolve through the scattering process into either
the same or other asymptotic free states of the theory at time $T=\infty$. The
$S$ matrix thus allows one to predict scattering cross-sections, which
represent the probability of a scattering event at a given impact parameter,
which is then be integrated over a range of impact parameters. The computation
of in-vacuum to the out-vacuum transitions embedded in the $S$ matrix is, for
$1+1$ dimensions, promise- complete. This means that if an arbitrary quantum
computation can be thought of as such a scattering experiment, and if a
classical computer could efficiently compute these matrix elements, then
quantum computers would be no more powerful than probabilistic classical
computers. In the language of complexity theory, = .

The computation of the entire $S$-matrix of a QFT is exponentially difficult,
although specific matrix elements can be computed within additive error in
polynomial time on a quantum computer. In [16] the authors provide a quantum
algorithm for computing such matrix elements by directly performing a
simulation and computing the overlaps between the outgoing distribution and
the incoming one. They approximate the solution by preparing well separated
wave packets, allowing them to evolve in time and scatter, and then measuring
the resulting products of the scattering.

Since the original work of Jordan, Lee and Preskill, quantum simulation
algorithms have undergone a mini revolution. The central problem in quantum
simulation algorithms is to translate the unitary time dynamics of a closed
quantum system into a sequence of quantum gates that, up to small
approximation error, implements the unitary. At the time of the original JLP
result, the dominant form of quantum simulation were product formulas [16, 3,
2] which approximate a complicated time-evolution operator as a product of
simple time-evolution operators (which are often so simple that they can be
synthesized from scratch into elementary gates. Since then a host of new
methods have been developed for performing simulation; amongst these,
qubitization [5, 21] has emerged as a leading non-Trotter based simulation
approach. Neither of the two methods is believed to fully advantageous
relative to the other. Qubitization provides provably optimal scaling in worst
case scenarios but is inflexible and is not suited to take advantage of
commutation relations between terms. Product (Trotter) formulas are the
opposite in this regard: they take advantage of small commutators but fail to
achieve optimal scaling in the worst case scenario.

Subsequent work has provided a number of optimized methods for simulating
various aspects of $\phi^{4}$ theory [11, 22]. There have also been
considerable progress for various other field theories [23, 24, 25].However,
at present we still do not know whether there exist practical simulation
algorithms for these methods in the sense that we do not yet know whether it
is possible to perform such simulations using a reasonable number of gate
operations and physical qubits on surface-code based quantum computers.

###  I.1 Our Contributions

We provide a cost-analysis of fault tolerant calculation of scattering matrix
elements using fixed volume methods for scalar $\phi^{4}$ theory in this paper
that is appropriate in general for $1+1$D and for certain low-energy elastic
collisions in higher dimensions. We achieve this through two main
contributions. First, we propose an alternative approach to JLP [16] by which
we are able to readily extract values of the $S$-matrix through the
measurement of multi-particle energies in finite volume. These techniques can
also be applied in certain cases to elastic collisions in higher dimensions.
Specifically, we discuss below (specifically in Claim 1) that direct
calculation of the elements of the $S$-matrix is not often needed and the low-
lying eigenvalues of the field Hamiltonian can be used to estimate, under
specific circumstances, values of the $S$-matrix that would otherwise be
impractical to compute directly as per [16, 17].

We implement our approach using a series of different fault tolerant
simulation algorithms, which form the second major contribution of this paper.
We consider primarily Hamiltonians formulated both in the field occupation
basis and field amplitude basis. The occupation basis Hamiltonian is simulated
with a Trotterization algorithm, while for the amplitude basis one we describe
four simulation procedures - one with Trotterization, and the three others
with modern qubitization methods. These qubitization approaches are optimized
through our introduction of new unitary block encodings (implemented through
LCU or linear combinations of unitaries approaches [3]) of $\phi$, $\phi^{2}$
and $\phi^{4}$ operators and in order to reduce the resource cost. Further,
the implemented circuits have been optimized with recent optimization
techniques [26, 27]. The Trotterization algorithm in the occupation basis
performs optimally in the non-interacting limit, but is outperformed by
qubitization methods (in the amplitude basis) with increasing interaction
strength. However, the main appeal of the occupation basis algorithm is that
the basis provides a natural extension for state preparation and measurement
of direct scattering calculations in higher dimensions. In Table 1 we have
compared the cost of implementing the various algorithms in terms of number of
qubits and T-gates used. Here, we mention that we have implemented our
algorithms with the Clifford+T gate set, which is the most popular, fault-
tolerant, universal gate set considered for quantum simulation algorithms [28,
29, 30, 15] and affords a wide variety of exact and approximate unitary
synthesis methods [31, 32, 33, 34, 35, 36]. At times, it can exactly
synthesize useful gates that could at best be approximately synthesized in
other gate libraries [37, 38, 39]. The resource overhead in fault-tolerant
implementation of the non-Clifford T gate is the highest in most error-
correction schemes, including the surface code. A bound on the T-count can be
can be used as a marker to reflect the complexity of fault-tolerant
implementation of a quantum algorithm[40].

We then consider implementing the Trotter and qubitization based simulation
algorithms on top of the surface code and consider the overheads of magic
state distillation in the simulation algorithm. The space-time volume needed
to implement a $T$ gate is expected to dwarf the costs of all other gates by
factors of $100$ or greater. Our calculations for surface code costs are based
on canonical approaches [40]. This facilitates the translation for others
towards any preferred modern methods, (including those that may exist in
future) in the ever-expanding forefront of quantum error correction [41, 42].

Algorithm Qubit $\\#$ T gate $\mathbf{\\#}$ Reference Occ. Trotter
$O\left(N|\Omega|+\log{\left(\frac{\lambda N}{M\epsilon_{E}}\right)}\right)$
$O\left(\frac{\lambda
N^{7}|\Omega|^{3}}{M^{5/2}\epsilon_{E}^{3/2}}\log(\lambda\frac{N|\Omega|}{M\epsilon_{E}})\right)$
Theorem 12 Amp. Trotter
$O\left(|\Omega|\log(k)+\log_{2}\left(\frac{|\Omega|\Lambda k\left(\Lambda
k+M^{2}\right)}{\epsilon_{E}}\right)\right)$
$\tilde{O}\left(\frac{|\Omega|^{3/2}\sqrt{\Lambda^{2}k^{5}+\Lambda
M^{2}k^{4}}\log^{4}(k)}{\epsilon_{E}^{3/2}}\right)$ Theorem 21, Lemma 210 Amp.
Qubitized I
$O\left(|\Omega|\log{k}+\log^{2}{k}+\log{\left(\frac{|\Omega|\left[k^{2}\Lambda+kM^{2}\right]}{\epsilon_{E}}\right)}\right)$
$O\left(\frac{|\Omega|^{2}}{\epsilon_{E}}\left[k^{2}\Lambda+kM^{2}\right]\log^{2}{k}\right)$
Theorem 34 Amp. Qubitized IIIa
$O\left(|\Omega|\log{k}+\log{\left[\frac{|\Omega|\left(k^{2}\Lambda+kM^{2}\right)}{\epsilon_{E}}\right]}\right)$
$O\left(\frac{|\Omega|^{2}}{\epsilon_{E}}\left[k^{2}\Lambda+kM^{2}\right]\left[\log^{4}{k}\right]\right)$
Theorem 34 Amp. Qubitized IIIb
$O\left(|\Omega|\log{k}+\log{\left[\frac{|\Omega|\left(k^{2}\Lambda+kM^{2}\right)}{\epsilon_{E}}\right]}\right)$
$O\left(\frac{|\Omega|^{2}}{\epsilon_{E}}\left[k^{2}\Lambda+kM^{2}\right]\left[\log^{2}{k}\right]\right)^{*}$
Proposition 35

Table 1: Qubit counts and T gate cost scalings for computing $\phi^{4}$
scattering matrix elements as a function of the lattice size $|\Omega|$, the
(rescaled) coupling constant $\Lambda$ and mass $M$, which hide dependence on
the lattice size, the amplitude cutoff $k$ in the field amplitude basis in
units of the bin size, the momentum cutoff $N$ in the occupation basis, and
the target precision in the energy estimate $\epsilon_{E}$, which dictates the
precision of the S-matrix through Eq. (38). The $\tilde{O}$ notation hides
further multiplicative factors at most logarithmic in all variables. (*) For
Algorithm IIIb the estimate on the T-count has been obtained assuming
Conjecture 28.

The organization of paper is as follows. In Section II we provide a brief
review of the scalar $\phi^{4}$ QFT. We specifically provide derivations of
the Hamiltonian for $\phi^{4}$ theory in both amplitude and occupation bases.
We then discuss how phase estimation can be used to estimate the scattering
amplitude in Section III. We use these techniques in Section IV to provide
quantum algorithms for estimating the elements of the $S$-matrix (under the
assumption of elastic collisions) in both the amplitude and the occupation
basis. This section gives not only the asymptotic scalings required but also
the constant factor analysis needed to estimate the energy within fixed error.
Section V contains our analysis of the implementation of the quantum
algorithms using the surface code along with the space-time needs for the
simulation. We finally conclude in Section VI and discuss future applications.

##  II Review of $\phi^{4}$ Scalar Field Theory

In this section we introduce notation for the $\phi^{4}$ scalar field theory.
In order to simulate scattering processes in the theory, we need a discrete
Hamiltonian representation of the continuum Lagrangian. We provide a concise
definition and derivation of the Hamiltonian in both a field and occupation
basis for the field theory. The regimes and limits where each basis may be
ideal are discussed in more detail below.

We begin by considering the free-theory for a scalar bosonic field, denoted by
$\phi(x)$. The first term in our Lagrangian density is the free Klein-Gordon
Lagrangian given by

$\displaystyle\mathcal{L}_{0}$
$\displaystyle=\frac{1}{2}\left(\partial_{\mu}\phi\right)^{2}-\frac{1}{2}m^{2}\phi^{2}=\frac{1}{2}\dot{\phi}^{2}-\frac{1}{2}(\nabla\phi)^{2}-\frac{1}{2}m^{2}\phi^{2}\,$
(1)

where $m$ is the bare mass of the scalar boson field. The corresponding
conjugate field is given by the Euler-Lagrange equations to state

$\pi(\mathbf{x})\equiv\frac{\partial\mathcal{L}_{0}}{\partial\dot{\phi}(\mathbf{x})}$
(2)

The simplest interaction to consider for this Lagrangian would be a quartic
interaction term given by

$\mathcal{L}_{int}=\lambda\phi^{4}$ (3)

The corresponding Hamiltonian for the theory is given as

$H=\int d^{d}x\mathcal{H}=\int
d^{d}x\left[\frac{1}{2}\pi(\mathbf{x})^{2}+\frac{1}{2}(\nabla\phi(\mathbf{x}))^{2}+\frac{1}{2}m^{2}\phi(\mathbf{x})^{2}+\lambda\phi(\mathbf{x})^{4}\right]$
(4)

Both formalisms considered here would also extend to interaction terms of
arbitrary order $\phi^{\alpha}$. These higher-order terms appear in general
phenomenological Landau theories of phase transitions in non-relativistic
condensed matter and have rich applications in nonlinear dynamics [18]. In
this work however we focus on the $\phi^{4}$ term more in the context of the
high-energy community. This interaction is the simplest renormalizable term
that is also bounded from below.

We consider a theory of scalar $\phi^{4}$ inside a $d$-dimensional space
confined to the $d$-cube of side length $L$, considering the momentum cut-off
= $\Lambda(\epsilon)$, with ${}^{d}\sqrt{|\Omega|}=P$ modes in each dimension.
For simplicity, we will refer to this volume region as $\Omega$ and we assume
that this region is discretized into $|\Omega|=(L/a)^{d}$ points with a
spacing between each of the points of $a$.

###  II.1 Field Amplitude Basis

We begin by considering our Hilbert space for the theory to be of the form

$\mathcal{H}=\bigotimes_{{\bf x}\in\Omega}\mathcal{H}_{\phi},$ (5)

where $\mathcal{H}_{\phi}$ is the Hilbert space that describes a field at each
of the lattice sites. The discretized lattice form of the Hamiltonian is

$H=\sum_{\mathbf{x}\in\Omega}a^{d}\left[\frac{1}{2}\pi(\mathbf{x})^{2}+\frac{1}{2}\left(\nabla_{a}\phi\right)^{2}(\mathbf{x})+\frac{1}{2}m^{2}\phi(\mathbf{x})^{2}+\frac{\lambda}{4!}\phi(\mathbf{x})^{4}\right]$
(6)

where $\Omega$ is the set of all (spatial) lattice sites,
$\mathbf{x}=a\mathbf{n}$, where $a$ is the lattice spacing and
$\mathbf{n}\in\mathbb{Z}^{\times d}$, and $d$ is the number of spatial
dimensions.

We then consider the conjugate momentum operator

$\pi(x)=\mathcal{F}^{\dagger}\phi(x)\mathcal{F},$ (7)

where $F$ is the discrete quantum Fourier transform acting on our truncated
space. Starting with Eq. (6), we identify

$\displaystyle\left(\nabla\phi(\mathbf{x})\right)_{i}\rightarrow\Delta_{i}\phi(\mathbf{x})$
$\displaystyle\equiv$
$\displaystyle\frac{\phi(\mathbf{x}+a\hat{x}_{i})-\phi(\mathbf{x})}{a}$
$\displaystyle\Rightarrow\sum_{i=1}^{d}\left(\Delta_{i}\phi(\mathbf{x})\right)^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{d}\frac{\phi^{2}(\mathbf{x}+a\hat{x}_{i})+\phi^{2}(\mathbf{x})-2\phi(\mathbf{x})\phi(\mathbf{x}+a\hat{x}_{i})}{a^{2}}$
(8)

If we assume periodic boundary conditions, then at each lattice site
$\mathbf{x}$, there is one contribution of $\phi(\mathbf{x})^{2}$ from each of
the $d$ lattice sites at $\mathbf{x}^{(i)}=\mathbf{x}-a\hat{x}_{i}$ (without
periodic boundary conditions, we would need to be careful about handling the
sites at the edges), and a single additional contribution coming from the
$\phi(\mathbf{x})^{2}$ term. Therefore, the Hamiltonian becomes

$H=\sum_{\mathbf{x}\in\Omega}\left[\frac{1}{2}a^{d}\pi^{2}(\mathbf{x})+\frac{1}{2}(d+1+a^{2}m^{2})a^{d-2}\phi^{2}(\mathbf{x})+\frac{\lambda}{4!}a^{d}\phi^{4}(\mathbf{x})-a^{d-2}\sum_{i=1}^{d}\phi(\mathbf{x})\phi(\mathbf{x}+a\hat{x}_{i})\right]$
(9)

Let us now define

$\displaystyle\Phi:=a^{\frac{d}{2}-1}\phi,$
$\displaystyle\Pi:=a^{\frac{d}{2}}\pi,$ $\displaystyle M:=am,\quad$
$\displaystyle\Lambda:=a^{4-d}\lambda$ (10)

Then, in terms of the scaled variables, our Hamiltonian simplifies to

$H_{amp}=\sum_{\mathbf{x}\in\Omega}\left[\frac{1}{2}\Pi^{2}(\mathbf{x})+\frac{1}{2}(M^{2}+d+1)\Phi^{2}(\mathbf{x})+\frac{\Lambda}{4!}\Phi^{4}(\mathbf{x})-\sum_{i=1}^{d}\Phi(\mathbf{x})\Phi(\mathbf{x}+a\hat{x}_{i})\right]$
(11)

Thus, we see that the derivative term is effectively replaced by a nearest
neighbor interaction. We shall make use of this expression, which we denote by
$H_{amp}$ in both the Qubitization-based and Trotterization based
implementation in the field amplitude basis (Section IV.5).

The field amplitude basis is defined to be that in which the field operator is
diagonal

$\Phi\ket{\phi}=\phi\ket{\phi}$ (12)

and, in the continuum, is related to the Fock basis via
$\ket{\phi}=\sum_{n=0}^{\infty}\ket{n}\bra{n}\phi\rangle$, where

$\bra{\phi}\ket{n}=\left(\frac{M}{\pi}\right)^{1/4}\frac{1}{\sqrt{2^{n}n!}}e^{-\frac{M\phi^{2}}{2}}H_{n}\left(\sqrt{M}\phi\right)$
(13)

where $H_{n}(x)$ is the (physicists’) Hermite polynomial of order $n$. This
basis has the advantage that all terms in the $\phi^{4}$ Hamiltonian, except
for the momentum term, are diagonal. As such, the basis should operate better
in large $\lambda$ coupling limits. The momentum operator, in turn, is related
to the field operator through a Fourier transform. We discuss Trotter and
qubitization approaches to simulating the $\phi^{4}$ Hamiltonian in this basis
in Section IV.5.

###  II.2 Field Occupation Basis

We note that the decomposition of the Hilbert in Eq. (5) is not the only way
to latticize the theory. In fact, despite its clear appeal when studying low
energy properties, such a discretization might not be ideal at high energies.
As we mentioned previously, an emergent particle description appears
(“partons”) in high energy experiments, and therefore in this context it may
be more appropriate to discretize the theory in this parton basis. Such a
description starts then by dividing the Hilbert space into single particle
sectors

$\mathcal{H}=\bigotimes_{{l}}^{\infty}\mathcal{H}_{l},$ (14)

where $\mathcal{H}_{l}$ describes the Hilbert of a single particle, which is
span by the local vacuum and all the possible momentum (or position) modes
that can be occupied. One then evolves the theory by freely propagating the
particles with the quadratic part of the Hamiltonian in each particle sector,
while the interactions connect different $l$’s. Such a procedure has been
applied to several problems in high energy physics, see for example [22, 43,
44, 45].

In the same spirit, the typical approach to solving the $\phi^{4}$ theory
begins with diagonalizing the non-interacting Klein-Gordon Hamiltonian. This
is done by introducing bosonic ladder operators
$a^{\dagger}_{\mathbf{p}},a_{\mathbf{p}}$ which correspond to the creation and
annihilation of quanta inside the scalar field. They obey the canonical
bosonic commutation relation

$[a_{\mathbf{p}},a^{\dagger}_{\bf{q}}]=(2\pi)^{D}\delta(p-q)$ (15)

where $D$ is the spatial dimensionality. However these operators on their own
are not Hermitian. As such, we construct Hermitian operators (and their
unitary exponentiation) through considering them in conjunction with their
Hermitian conjugates. As such, we begin with their implementation.

To carry out our analysis in this basis, we require the mode expansion of the
scalar field and its conjugate momentum

$\displaystyle\phi(\mathbf{x})$ $\displaystyle=$
$\displaystyle\int\frac{d^{D}p}{(2\pi)^{D}}\frac{e^{i\mathbf{p}\cdot\mathbf{x}}}{\sqrt{2\omega_{\mathbf{p}}}}\left(a_{\mathbf{p}}+a_{-\mathbf{p}}^{\dagger}\right)$
$\displaystyle\pi(\mathbf{x})$ $\displaystyle=$
$\displaystyle\int\frac{d^{D}p}{(2\pi)^{D}}(-i)e^{i\mathbf{p}\cdot\mathbf{x}}\sqrt{\frac{\omega_{\mathbf{p}}}{2}}\left(a_{\mathbf{p}}-a_{-\mathbf{p}}^{\dagger}\right)$
(16)

where $\omega_{\mathbf{p}}=\sqrt{|\mathbf{p}|^{2}+M^{2}}$, and $D$ denotes the
spatial dimensionality. The non-interacting part of the scalar field
Hamiltonian, or the Klein-Gordon Hamiltonian $H_{0}$, is then given by

$\displaystyle H_{0}$ $\displaystyle=$ $\displaystyle\int
d^{D}x\frac{1}{2}\left[\left(\pi(\mathbf{x})\cdot\pi(\mathbf{x})\right)+\left(\nabla\phi(\mathbf{x})\cdot\nabla\phi(\mathbf{x})\right)+M^{2}\left(\phi(\mathbf{x})\cdot\phi(\mathbf{x})\right)\right]$
(17) $\displaystyle=$ $\displaystyle\int
d^{D}x\int\frac{d^{D}\mathbf{p}d^{D}\mathbf{p}^{\prime}}{(2\pi)^{D}}\frac{e^{i(\mathbf{p}+\mathbf{p^{\prime}})\cdot\mathbf{x}}}{(2\pi)^{D}}\left\\{-\frac{\sqrt{\omega_{\mathbf{p}}\omega_{\mathbf{p}^{\prime}}}}{4}\left(a_{\mathbf{p}}-a_{-\mathbf{p}}^{\dagger}\right)\left(a_{\mathbf{\mathbf{p}^{\prime}}}-a_{-\mathbf{\mathbf{p}^{\prime}}}^{\dagger}\right)\right.$
$\displaystyle\left.+\frac{m^{2}-\mathbf{p}\cdot\mathbf{p^{\prime}}}{4\sqrt{\omega_{\mathbf{p}}\omega_{\mathbf{p^{\prime}}}}}\left(a_{\mathbf{p}}+a_{-\mathbf{p}}^{\dagger}\right)\left(a_{\mathbf{p^{\prime}}}+a_{-\mathbf{p^{\prime}}}^{\dagger}\right)\right\\}$
$\displaystyle=$ $\displaystyle\int\frac{d^{D}p}{(2\pi)^{D}}\
\omega_{\mathbf{p}}\left(a_{\mathbf{p}}^{\dagger}a_{\mathbf{p}}+\frac{1}{2}[a_{\mathbf{p}},a_{\mathbf{p}}^{\dagger}]\right)$

The second term in the parentheses above is proportional to $\delta(0)$, and
represents the zero point energy of all the harmonic oscillator modes in the
free scalar field. This formally leads to an infinite contribution to the
Hamiltonian, which must be removed to compute physical quantities. The
resulting free Hamiltonian we consider is

$H_{0}=\int\frac{d^{D}\pi}{(2\pi)^{D}}\ \omega_{\mathbf{p}}\
a_{\mathbf{p}}^{\dagger}a_{\mathbf{p}}$ (18)

This Hamiltonian is now a normal ordered product of ladder operators, where
all creation operators appear on the left of all annihilation operators. This
reordering is necessary to avoid spurious divergences and is denoted
$H_{0}=:H_{0}:$.

This form has two primary appeals. One is that in the limit of small
$\lambda$, our approach approaches a trivial diagonal matrix. This limit up to
arbitrary precision remains a challenging computational task. The primary
appeal of this approach is that for explicit scattering process that require
single-particle packet preparation, this protocol also has much simpler state
preparation procedure. Next, the interaction part of the Hamiltonian is given
by

$\begin{gathered}H_{\lambda}=\int d^{D}x\frac{\lambda}{4!}\phi^{4}\\\
=\frac{\lambda}{4!}d^{D}x\int\frac{d^{D}\mathbf{p}_{1}d^{D}\mathbf{p}_{2}d^{D}\mathbf{p}_{3}d^{D}\mathbf{p}_{4}}{(2\pi)^{4D}4\sqrt{\omega_{\mathbf{p}_{1}}\omega_{\mathbf{p}_{2}}\omega_{\mathbf{p}_{3}}\omega_{\mathbf{p}_{4}}}}\left(a_{\mathbf{p}_{1}}e^{+i\mathbf{p}_{1}\cdot\mathbf{x}}+a^{\dagger}_{\mathbf{p}_{1}}e^{-i\mathbf{p}_{1}\cdot\mathbf{x}}\right)\left(a_{\mathbf{p}_{2}}e^{+i\mathbf{p}_{2}\cdot\mathbf{x}}+a^{\dagger}_{\mathbf{p}_{2}}e^{-i\mathbf{p}_{2}\cdot\mathbf{x}}\right)\\\
\times\left(a_{\mathbf{p}_{3}}e^{+i\mathbf{p}_{3}\cdot\mathbf{x}}+a^{\dagger}_{\mathbf{p}_{3}}e^{-i\mathbf{p}_{3}\cdot\mathbf{x}}\right)\left(a_{\mathbf{p}_{4}}e^{+i\mathbf{p}_{4}\cdot\mathbf{x}}+a^{\dagger}_{\mathbf{p}_{4}}e^{-i\mathbf{p}_{4}\cdot\mathbf{x}}\right)\\\
=\int\frac{\lambda}{4!}d^{D}x\int\frac{d^{D}\mathbf{p}_{1}d^{D}\mathbf{p}_{2}d^{D}\mathbf{p}_{3}d^{D}\mathbf{p}_{4}}{(2\pi)^{4D}}\frac{e^{i(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3}+\mathbf{p}_{4})\cdot\mathbf{x}}}{4\sqrt{\omega_{\mathbf{p}_{1}}\omega_{\mathbf{p}_{2}}\omega_{\mathbf{p}_{3}}\omega_{\mathbf{p}_{4}}}}\\\
\left(a_{\mathbf{p}_{1}}+a_{-\mathbf{p}_{1}}^{\dagger}\right)\left(a_{\mathbf{p}_{2}}+a_{-\mathbf{p}_{2}}^{\dagger}\right)\left(a_{+}a_{-\mathbf{p}_{3}}^{\dagger}\right)\left(a_{\mathbf{p}_{4}}+a_{-\mathbf{p}_{4}}^{\dagger}\right)\\\
=\frac{\lambda}{4!}\int\frac{d^{D}\mathbf{p}_{1}d^{D}\mathbf{p}_{2}d^{D}\mathbf{p}_{3}}{(2\pi)^{3D}}\frac{1}{4\sqrt{\omega_{\mathbf{p}_{1}}\omega_{\mathbf{p}_{2}}\omega_{\mathbf{p}_{3}}\omega_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}}}\left\\{\left(a_{\mathbf{p}_{1}}+a_{-\mathbf{p}_{1}}^{\dagger}\right)\left(a_{\mathbf{p}_{2}}+a_{-\mathbf{p}_{2}}^{\dagger}\right)\right.\\\
\times\left.\left(a_{\mathbf{p}_{3}}+a_{-\mathbf{p}_{3}}^{\dagger}\right)\left(a_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}+a_{(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}^{\dagger}\right)\right\\}\end{gathered}$
(19)

For the algorithm, we consider the normal ordered form $:H_{\phi}:$, discussed
in 240

We consider the scalar field on a finite sized, discrete lattice, borrowing
the notation from above which requires us to replace the continuous integrals
with discrete Riemann sums as

$\int_{0}^{2\pi/a}\frac{d^{D}p}{2\pi}\rightarrow\sum_{\alpha=1}^{|\Omega|}\frac{(2\pi)}{(2\pi)|\Omega|(a)^{D}}=\frac{1}{(a)^{D}}\sum_{\alpha=1}^{|\Omega|}\frac{1}{|\Omega|}$
(20)

Now the Hamiltonians reads as follows

$:H_{0}:=\frac{1}{|\Omega|}\sum_{\mathbf{p}}\omega_{\mathbf{p}}a^{\dagger}_{\mathbf{p}}a_{\mathbf{p}}$
(21)

We note that $(a)^{D}|\Omega|=L$, but we retain this notation for if
$|\Omega|$ increases with fixed $L$, i.e, taking the $\alpha\rightarrow 0$
continuum limit this notation makes the scaling explicit. Furthermore, we take
proper discretization of our operators
$a_{\mathbf{p}}\rightarrow\sqrt{a}a_{\mathbf{p}}$ so that these are
dimensionless. [46]. An advantage of the momentum space-representation is that
it allows for the retention of the continuum dispersion as prepared to in
position space [22].

The $\lambda$ term is equivalently given by

$\begin{gathered}H_{\lambda}=\frac{\lambda}{4!|\Omega|^{3}}\sum_{\\{\mathbf{p}_{i}\\}}\frac{1}{4\sqrt{\omega_{\mathbf{p}_{1}}\omega_{\mathbf{p}_{2}}\omega_{\mathbf{p}_{3}}\omega_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}}}\left\\{\left(a_{\mathbf{p}_{1}}+a_{-\mathbf{p}_{1}}^{\dagger}\right)\left(a_{\mathbf{p}_{2}}+a_{-\mathbf{p}_{2}}^{\dagger}\right)\left(a_{-\mathbf{p}_{3}}+a_{\mathbf{p}_{3}}^{\dagger}\right)\left(a_{(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}+a_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}^{\dagger}\right)\right\\}\end{gathered}$
(22)

We leave the normal ordering until we consider the cases of momentum
$\mathbf{p}_{i}$ overlaps.

As a final note, due to the relativistic dispersion relationship the frequency
$\omega_{\bf p}$ depends on the momentum as well as the mass via $\omega_{\bf
p}=\sqrt{M^{2}+|{\bf p}|^{2}}$. For this reason the frequencies
$\omega_{\bf{p}}$ subtley contain the mass. If we are interested in the low
momentum regime then we can take the maximum and minumum relevant frequencies
to be $\omega_{min}=M$ and $\omega_{max}=\sqrt{M^{2}+P_{max}^{2}}\sim M$ which
gives us the required mass dependence of the Hamiltonian.

##  III Scattering Matrix Estimation

The cost of any end-to-end quantum simulation depends on the observables that
we wish to measure in addition to the processing required to prepare the
initial state and simulate the dynamics of the quantum system. Unlike previous
work, which directly examines scattering in field theories [16, 17], energy
estimation is computationally simple and opens up the possibility of practical
simulations of quantum dynamics.

We will now go into detail about our alternative approach to that of [16, 17],
for the computation of values of the scattering matrix or $S$-matrix. The
scattering matrix describes the relationship between the incoming and outgoing
momenta for a quantum system. In quantum field theory in particular, it is
useful to define the $S$-matrix in terms of its dynamics in an interaction
frame. In particular, if we define the Hamiltonian for a field to be

$H=H_{\rm free}+H_{\rm int}$ (23)

where $H_{\rm free}$ and $H_{\rm int}$ are the free and interaction terms
where $H_{\rm int}$ vanishes at infinity. In order to unambiguously define the
dynamical phase of the non-interacting system as $T\rightarrow\infty$ one can
go into an interaction frame of $H_{\rm free}$ to obtain the interaction
picture Hamiltonian

$H_{I}(t)=e^{iH_{\rm free}t}H_{\rm int}e^{iH_{\rm free}t}.$ (24)

With this definition in mind we can express the $S$-matrix through its matrix
elements as

$\displaystyle S_{\vec{p},\vec{k}}=_{\text{ out }}\left\langle
p_{1}p_{2}\cdots\mid k_{\mathcal{A}}k_{\mathcal{B}}\right\rangle_{\text{in }}$
$\displaystyle=\lim_{T\rightarrow\infty}\langle\underbrace{p_{1}p_{2}\cdots}_{T}\mid\underbrace{k_{\mathcal{A}}k_{\mathcal{B}}}_{-T}\rangle$
(25) $\displaystyle=\lim_{T\rightarrow\infty}\left\langle
p_{1}p_{2}\cdots\left|\mathcal{T}e^{-i\int_{-\infty}^{\infty}dtH_{I}(t)}\right|k_{\mathcal{A}}k_{\mathcal{B}}\right\rangle.$

More specifically, by isolating the interaction terms $S=1+iT$, you can
identify

$\displaystyle\left\langle p_{1}\cdots
p_{n}|iT|p_{\mathcal{A}}p_{\mathcal{B}}\right\rangle$ (26)
$\displaystyle=\lim_{\mathcal{T}\rightarrow\infty(1-i\epsilon)}\left({}_{0}\left\langle
p_{1}\cdots
p_{n}\left|T\left(\exp\left[-i\int_{-T}^{T}dtH_{I}(t)\right]\right)\right|p_{\mathcal{A}}p_{\mathcal{B}}\right\rangle_{0}\right)$

The most important element of $T$ can be given through the matrix element
($\mathcal{M}$), which can be though of as the scattering amplitude and
includes the important dynamics.

$\left\langle
p_{1}p_{2}\cdots|iT|k_{\mathcal{A}}k_{\mathcal{B}}\right\rangle=(2\pi)^{4}\delta^{(4)}\left(k_{\mathcal{A}}+k_{\mathcal{B}}-\sum
p_{f}\right)\cdot i\mathcal{M}\left(k_{\mathcal{A}},k_{\mathcal{B}}\rightarrow
p_{f}\right)$ (27)

This $\mathcal{M}$ will give all the necessary calculations for the cross-
section if necessary. For reference: To one-loop order, we can compute
$\mathcal{M}$ in $\varphi^{4}$ as

${}_{0}\left\langle p_{1}p_{2}\left|T\left(-i\frac{\lambda}{4!}\int
d^{4}x\varphi_{I}^{4}(x)\right)\right|p_{\mathcal{A}}p_{\mathcal{B}}\right\rangle_{0}=-i\lambda(2\pi)^{4}\delta^{(4)}\left(p_{\mathcal{A}}+p_{\mathcal{B}}-p_{1}-p_{2}\right)$
(28)

which identifies $\mathcal{M}=-\lambda$.

Computing the entire $S$-matrix is exponentially difficult since the number of
output momenta and input momenta scales exponentially with the number of
particles permitted. But we can readily extract at least the elastic parts of
the S-matrix through the measurement of multi-particle energies in finite
volume. We will also comment on the ability to obtain inelastic information.
Here we will argue that the measurement of energies in finite volume leads
directly to information on the S-matrix elements. The simplest case where this
occurs is for a non-interacting $1+1D$ theory and where both particles have
mass $m$. In this case, we can go into a center of mass frame where without
loss of generality both particles have momenta $p$ and $-p$ such that
$p_{1}=-p_{2}$ and by conservation of momentum the output momenta also must be
in $\\{-p,p\\}$ In finite volume, the momenta, $p_{i},i=1,2$ of the particles
are

$p_{i}=\frac{2\pi n_{i}}{L},~{}~{}n_{i}\in\mathbb{Z},$ (29)

where $n_{i}$ are integers indexing the two particles and $L$ is the length of
the system. The energy, $E$, of the two particle state is

$E=\sum_{i=1,2}\left(m^{2}+\frac{4\pi^{2}n_{i}^{2}}{L^{2}}\right)^{1/2},$ (30)

where we have taken $c=\hbar=1$. Now what happens when interactions are turned
on? Provided the energy of the two-particle state is below threshold for
particle production, the quantization condition for the two momenta is altered
to

$\displaystyle e^{ip_{1}L}S(p_{1},p_{2})$ $\displaystyle=$ $\displaystyle
1=e^{2\pi n_{1}/L};$ (31) $\displaystyle e^{ip_{2}L}S(p_{2},p_{1})$
$\displaystyle=$ $\displaystyle 1=e^{2\pi n_{2}/L},$ (33)

where $S(p_{1},p_{2})=e^{i\delta(p_{1},p_{2})}$ is the two-body scattering
phase, which is a a phase because the scattering is elastic in $1+1$
dimensions.

By taking logarithms, the quantization conditions can be written in the form

$\displaystyle\frac{2\pi n_{1}}{L}$ $\displaystyle=p_{1}-\frac{i}{L}\log
S(p_{1},p_{2});$ (34) $\displaystyle\frac{2\pi n_{2}}{L}$
$\displaystyle=p_{2}-\frac{i}{L}\log S(p_{2},p_{1}).$ (36)

If one solves the quantization condition for $p_{1}$ and $p_{2}$ one
immediately knows the energy of the two-particle state via

$E=\sum_{i=1,2}(m^{2}+p_{i}^{2})^{1/2},$ (37)

where $p_{i}$ necessarily deviate from their free values. Now if instead we
start with knowledge of $E$ via measurement, we can reverse the process and
infer $S(p_{1},p_{2})$. And by measuring $E$ for different $L$, we can
determine $S(p_{1},p_{2})$ over a range of $p_{1}$ and $p_{2}$ because while
$p_{i}$ are not free, they still will behave as $1/L$ over a wide range of
volumes.

In any determination of the energy, $E$, of a two particle state en route to
the determination of an S-matrix element, there will be some uncertainty
associated with its determination, $\updelta E$. If we work in the center-of-
mass frame and the particles have equal and opposite momentum, i.e.,
$p_{1}=-p_{2}$, the associated uncertainty in the scattering phase
$\delta(p,-p)$ is given by

$\updelta\delta(p,-p)=-\frac{LE}{8p}\updelta E$ (38)

We can see that uncertainty in the determination of the energy affects most
dramatically the determination of the scattering phase at small momentum.

As an illustration of this in the $\varphi^{4}$ theory, we will consider an
example from the ordered phase of the model (i.e., $m^{2}<0$). The spectrum of
the theory in this regime consists of kinks and bound states of kinks (meson-
like states). As $\lambda$ is increased from 0, the bound states disappear at
a $\lambda_{0}$ from the spectrum and only kink-like states exist. At a large
enough $\lambda=\lambda_{c}$, the model becomes disordered. In our example, we
will consider the regime $\lambda_{0}<\lambda<\lambda_{c}$.

In Fig. 1, we plot the low-lying energy levels with zero total momenta as a
function of system size, $L$, as computed in Ref.[47] using Hamiltonian
truncation methods. The ground state energy has been subtracted. The lowest
lying level is the state whose energy is nearly degenerate to the ground state
(we expect such a near-degeneracy in finite volume in the broken phase).
Beyond this state are energies corresponding to two-kink and four-kink states.
We can use the two-kink energies to back out the scattering phase as described
above. This phase is plotted in Fig. 2 as a function of energy, $E$. Here the
energy of the two-particle state is parameterized by $\theta$ via :

$E=2M_{kink}\cosh(\theta)$

In this parameterization the momenta of the two kinks are $\pm
M_{kink}\sinh(\theta)$.

Figure 1: The low lying spectrum of the $\varphi^{4}$ in its broken phase
($m^{2}=-0.25$, $\lambda=0.15$). The ground state energy has been subtracted.
From Fig. 10 of [47]. Figure 2: The scattering phase inferred from the
energies of the two kink states presented in Fig. 1. From Fig. 11 of [47].

What we have discussed so far concerns the elastic part of the $2-2$
scattering matrix. We can also access inelastic information from the
measurement of energies. We can parameterize the S-matrix $S(p_{1},p_{2})$
solely in terms of $\theta$ ($\theta$ is related to the Mandelstam variable,
$s$, via $s=4m^{2}\cosh^{2}(\theta/2)$). The unitarity condition on the
S-matrix then reads

$S(\theta+i\epsilon)S(-\theta+i\epsilon)=f(\theta),$ (39)

where $f(\theta)$ is real positive on the real line. If $\theta_{0}$ marks the
threshold beyond which inelastic processes are possible,
$f(|\theta|<\theta_{0})=1$, while $0<f(|\theta|>\theta_{0})<1$. Using the
analytic structure of the S-matrix in the complex-$\theta$ plane, we can write
$S(\theta)$ in the following form:

$S(\theta)=\pm\prod_{j}\frac{\sinh(\theta)+i\sin(\beta_{j})}{\sinh(\theta)-i\sin(\beta_{j})}\exp[-\int^{\infty}_{-\infty}\frac{d\theta^{\prime}}{2\pi
i}\frac{\log
f(\theta^{\prime})}{\sinh(\theta-\theta^{\prime}+i\epsilon)}\bigg{]}.$ (40)

The first part of the parameterization involving the angles $\beta_{j}$ are
so-called Castillejo-Dalitz-Dyson (CDD) factors. The second part of the
parameterization contains information about the inelastic part of the 2-2
scattering inasmuch as it depends on the region in $\theta$ where $f(\theta)$
differs from 1. Notice that because of analyticity, the inelastic part of the
S-matrix influences the elastic region $|\theta|<\theta_{0}$. With
sufficiently accurate measurements of the energies, one can use this
parameterization to determine the number and values of the CDD angles,
$\beta_{j}$, as well as the function $f(\theta)$. In practice, extracting the
function $f(\theta)$ is difficult as it involves the inversion of an ill-posed
integral equation. However the integrated quantity,

$\int^{\infty}_{-\infty}\frac{d\theta^{\prime}}{2\pi i}\frac{\log
f(\theta^{\prime})}{\sinh(\theta-\theta^{\prime}+i\epsilon)},$ (41)

has been successfully extracted from numerical data, for example, in the case
of the non-integrable Ising field theory [48].

While we have focused our discussion here on scattering in $1+1d$, our
approach can be extended to scattering in higher dimensions, i.e., the elastic
scattering phases of 2-2 particle scattering, $\delta_{l}$, in different
angular momentum channels can be connected to measurements of the two-particle
energy in finite volume [49]. While originally formulated in [49] for two
spinless identical particles, extensions have been made to identical particles
with spin, asymmetric volumes, and amplitudes containing external current -
see [50] for a review. In certain cases, inelastic information involving 3
particle scattering is available in higher dimensions [51, 52, 53, 54].

Note that the last approach differs from that taken by Jordan Lee and Preskill
[16, 17]. Their approach involved preparing a Gaussian wave packet in the
interacting eigenbasis rather than an eigenvector of the interacting
Hamiltonian via phase estimation. Their approach has the advantage that it can
be applied in circumstances where the gap is large and no prior knowledge of
the eigenvectors is given. Directly preparing the eigenstates as we discuss
above can be more computationally efficient than the approach given by JLP,
but requires efficient approximations to the eigenstates which may not be
available in all settings.

To summarize, we see from this discussion that the excited state energy is a
meaningful quantity for scattering theory that can be used to estimate, under
specific circumstances, values of the $S$-matrix that would otherwise be
impractical to compute. The circumstances under which the energy can be used
to provide information about the scattering is given in the following claim
with increasing levels of generalization:

######  Claim 1.

Let $\tilde{S}$ be a finite dimensional version of the $S$-matrix for a
Hamiltonian operator $H(t)$. The value of the excited state energies can be
used to infer certain elements of $\tilde{S}$ if one of the following occurs

  1. 1.

We are interested in $2\mapsto 2$ elastic scattering of particles in one
spatial and one temporal dimension ($1+1\text{D}$) where the two-particle
final state is below the threshold for particle production.

  2. 2.

We are interested in inelastic processes in $1+1\text{D}$ in integrated form
(as encoded in Eq.(40).

  3. 3.

We are interested in computing the scattering phases of elastic $2\mapsto 2$
scattering in higher dimensions [49].

  4. 4.

We are interested in certain relatively simple scattering processes that
include particle decay in dimensions higher than $1+1\text{D}$. This is
focused primarily on scattering involving 3-particles, an active topic of
current research [51, 52, 53, 54].

Barring the cases mentioned in the above claim, there are other reasons why
the energies of both the ground state and low lying excited states is
independently interesting for the theory. It can allow us to understand phase
transitions in the strongly interacting theory [55, 56]. For all the above
reasons, we focus our attention on the computation of the energy of the low
lying excited states which is also a well studied quantity to estimate using
quantum computers [28, 30].

##  IV Simulation Algorithms

Modern simulation algorithms have converged in recent years to the point that
there is no optimal single quantum simulation algorithm. Rather, different
algorithms tend to have advantages and disadvantages in different regimes [15,
57, 58]. In this spirit, we provide five simulation algorithms for scalar
field theory in the occupation as well as field amplitude basis. These
algorithms employ Trotter decompositions or qubitization and use optimized
circuits to construct the operator exponentials or oracle calls that both
methods require. We further observe that the methods, as expected, have
advantages and disadvantages with respect to each other. Most obviously, the
Trotter algorithm for simulations in the occupation basis is by far the most
efficient for weak coupling ($\lambda\ll 1$) but the qubitized field amplitude
basis is likely to scale better in the strong coupling regime.

The material in this section is laid out as follows. In Section IV.1 we
describe an algorithm for the Hamiltonian formulated in the field occupation
basis, while in Section IV.5 we describe several algorithms for the
Hamiltonian in the field amplitude basis. The techniques and concepts in these
two sections can be followed independently of the other.

###  IV.1 Field Occupation Basis

The interaction Hamiltonian can be decomposed into four cases based on
matching momentum indices. We consider the normal ordered Hamiltonian so that
the unitary phases corresponding to eigenvalues that are used to construct the
$S-$matrix are computed with reference to the vacuum energy. First we state
some essential results that are required to map the Hamiltonian from the
bosonic to qubit space. Then we synthesize quantum circuits implementing the
exponentiated Hamiltonian.

We discretize the Hilbert space of occupation into $N$ distinct momentum
states. We have a register of $(N+1)V$ qubits to store information about the
occupation states and we use the following one-hot unary encoding to map an
occupation state to a qubit state. We index the qubits by a pair of integers,
such as $(p,n)$, where $p$ corresponds to a momentum mode and $n$ to a
momentum state. For each such pair $(p,n)$, we have a quantum state on
$(N+1)V$ qubits, in which each qubit is $\ket{0}$, except the
$(p,n)^{\text{th}}$ one, which is $\ket{1}$. We denote this state by
$\ket{p,n}$, which is

$\displaystyle\ket{p,n}$ $\displaystyle=$
$\displaystyle\ket{0_{1,0},\ldots,0_{p-1,N};0_{p,0},\ldots,0_{p,n-1},1_{p,n},0_{p,n+1},\ldots,0_{p,N};0_{p+1,0},\ldots
0_{V,N}}$ (42) $\displaystyle=$
$\displaystyle\left(\otimes_{j=1}^{p-1}\ket{0_{j,0}\ldots,0_{j,N}}\right)\bigotimes\ket{0_{p,0},\ldots
1_{p,n}\ldots,0_{p,N}}\bigotimes\left(\otimes_{j=p+1}^{|\Omega|}\ket{0_{j,0}\ldots,0_{j,N}}\right).$

We emphasis that each Hilbert subspace $\mathcal{H}\mathbf{p}$ is spanned only
by vectors of Hamming weight 1 in order to ensure the unary encoding. We now
consider a construction of the Hamiltonian that preserves this Hamming weight
without restricting to total number conservation or conservation for each
$p-$mode.

For convenience, we denote an operator $A_{pn}$ acting on $pn^{th}$ qubit by
$A_{p,n}$ or $\left(A_{n}\right)\mathbf{p}$. The qubit mapping for the
creation and annihilation operators is as follows.

$\displaystyle a_{\mathbf{p}}^{\dagger}$ $\displaystyle=$
$\displaystyle\sum_{n}\sqrt{n+1}\left(\sigma_{n}^{-}\sigma_{n+1}^{+}\right)_{\mathbf{p}}$
$\displaystyle a_{\mathbf{p}}$ $\displaystyle=$
$\displaystyle\sum_{n}\sqrt{n+1}\left(\sigma_{n}^{+}\sigma_{n+1}^{-}\right)_{\mathbf{p}},$
(43)

where

$\sigma^{+}=\frac{1}{2}(X-iY),\quad\sigma_{-}=\frac{1}{2}(X+iY).$ (44)

and therefore

$\displaystyle
a_{\mathbf{p}}^{\dagger}\ket{p,n}=\sqrt{n+1}\ket{p,n+1}\qquad\text{and}\qquad
a_{\mathbf{p}}\ket{p,n}=\sqrt{n}\ket{p,n-1}$ (45)

Now, considering Hermitian pairing of operators, we have

$\displaystyle a_{\mathbf{p}}+a_{\mathbf{p}}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n}\sqrt{n+1}\left(X_{n}X_{n+1}+Y_{n}Y_{n+1}\right)_{\mathbf{p}}.$
(46)

In theory the number of momentum states range till infinity, but for our
simulation, we truncate the Hilbert space and have $N$ momentum states, thus
$n$ varies from $0$ to $N$ in the above summations. This truncation in the
bosonic occupation is proportional to the maximum energy expected to be
simulated in semi-elastic collisions $N\propto\frac{E}{\omega_{\mathbf{p}}}$.
It is expected that the error in this truncation is exponentially small with
respect to the cutoff [59, 60, 61, 7]. However, careful numerical analysis of
these (and other finite-system size effects) will be required in
implementation [62, 63].

The details of the proofs for the following lemmas have been provided in
Appendix A

A similar approach was laid out in [22]. We note several distinctions between
the current implementation and previous work. The previous work has two
conceptual appeals. One is the use of a binary representation for the
occupation basis of $a_{\mathbf{p}}+a^{\dagger}_{\mathbf{p}}$. Further work
may examine efficient fault-tolerant arithmetic primitives required for such
an encoding, which we avoid here. The second is a diagonalization of the
position space field amplitude

$\phi_{\mathbf{n}}^{(i)}\equiv\frac{a_{\mathbf{n}}^{(i)}+a_{\mathbf{n}}^{(i)\dagger}}{\sqrt{2}}$
(47)

where

$a_{\mathbf{n}}\equiv\frac{1}{\sqrt{\mathcal{V}}}\sum_{\mathbf{q}}A_{\mathbf{q}}e^{i2\pi\mathbf{n}\cdot\mathbf{q}/N_{s}}$
(48)

and

$A_{\mathbf{q}}\equiv\frac{1}{2}\left[\omega_{\mathbf{q}}^{-\frac{1}{2}}+\omega_{\mathbf{q}}^{\frac{1}{2}}\right]a_{\mathbf{q}}+\frac{1}{2}\left[\omega_{\mathbf{q}}^{-\frac{1}{2}}-\omega_{\mathbf{q}}^{\frac{1}{2}}\right]a_{-\mathbf{q}}^{\dagger}$
(49)

This is implemented with a squeezing operator framework. In a fault-tolerate
framework, the phase factors of the squeezing operators introduce an
additional error without a clear mapping to the final error. These squeezing
factors also require a large number of controlled $R_{Z}$ gates Extending such
an elegant implementation into the fault-tolerant regime would be a intriguing
further work.

######  Lemma 2.

If $\hat{n}_{p}$ is the number operator on momentum mode $p$ then for any
integer $r\geq 1$ we have,

$\left(\hat{n}_{p}\right)^{r}=\sum_{n}\frac{n^{r}}{2}(I_{n}-Z_{n})_{p}.$ (50)

######  Lemma 3.

If $m\geq 1$ and $r\geq 0$ are integers then we have,

$(a_{p}^{\dagger})^{m}(n_{p})^{r}+(n_{p})^{r}(a_{p})^{m}=\frac{1}{2}\sum_{n}\sqrt{\frac{(n+m)!}{n!}}n^{r}(X{{}_{n}}X_{n+m}+Y_{n}Y_{n+m})_{p}.$
(51)

We first consider the non-interacting term in Eq.17 and derive the circuit
required to simulate the unitary exponential, with the following result.

######  Lemma 4 (Complexity of non-interacting Hamiltonian Simulation).

We require at most $N|\Omega|$ number of $R_{z}$ gates to simulate
$e^{-iH_{0}^{\prime}t}$, where $H_{0}^{\prime}=:H_{0}:$.

###### Proof.

The normal ordered $H_{0}$ i.e. $:H_{0}:$ is,

$:H_{0}:=\sum_{p}\omega_{p}a^{\dagger}_{p}a_{p}=\sum_{p}\omega_{p}\hat{n}_{p}=\sum_{n,p}\frac{n\omega_{p}}{2}(I_{n}-Z_{n})_{p},$
(52)

using Lemma 50. The above expression, up to some global phase is equal to a
sum of Z operators. And the exponential of each of them can be implemented
with a $R_{z}$ gate, whose angle depends on the coefficient. Thus the lemma
follows. ∎

##### Circuits for simulating the interacting Hamiltonian :

We now consider the interacting part of the Hamiltonian i.e. $H_{\varphi}$,
which we are re-writing as follows, for convenience. Let
$S_{4\mathbf{p}}=\\{\mathbf{p}=(p_{1},p_{2},p_{3},p_{4}):p_{i}\in\Gamma;i=1,2,3,4\\}$,
be a set consisting of ordered 4-tuples of the momentum mode, that respects
the conservation of momentum constraint ($p_{1}+p_{2}=\pm(p_{3}+p_{4})$). We
parameterize this constraint by working in the basis such that such that $\pm
p_{3}=p_{1}+k$ and $\pm p_{4}=p_{2}-k$.

Now, we can divide the terms in the above sum into 4 groups, based on equality
of the momentum modes in $\mathbf{p}$, i.e.
$H_{\varphi}:=H_{1\varphi}+H_{2\varphi}+H_{3\varphi}+H_{4\varphi}$. For each
such group, we map the resulting bosonic expression into the qubit space using
Equation 46, Lemma 50, Lemma 3, and obtain the Pauli expression. From this, we
derive the quantum circuit for the exponentiated sum in each group. We follow
the methods described in [26]. In short, we first derive an eigenbasis for a
set of mutually commuting Pauli terms and calculate the number of distinct
non-zero eigenvalues (ignoring sign) for the required (weighted) sum of these
Paulis. This number is equal to the number of (controlled)-rotations we
require. After that, we design the explicit quantum circuit, with appropriate
optimizations. Each circuit is implemented with
Clifford+T+(controlled)-$R_{z}$ gates.

We have delegated the detailed description of the circuit design to Appendix A
and in this section we only mention the relevant results required for
estimating the gate counts. As indicated, we have separated the T-gates
arising from the approximately implementable (controlled)-rotation gates and
the exactly implementable Toffoli or multi-controlled-X gates. So we report
the number of (controlled)-$R_{z}$ required and the number of T gates in Lemma
5-8 do not include the T-gates from these rotations. The total T-count
estimate is obtained by multiplying the number of (controlled)-$R_{z}$ gates
with the T-count of individual rotation gates and adding the product to the
T-count arising from the exactly implementable parts.

In the following discussions let
$S_{4\vec{n}}=\\{\vec{n}=(n_{1},n_{2},n_{3},n_{4}):n_{i}=0,1,\ldots,N;i=1,2,3,4\\}$
be an ordered 4-tuple of momentum states.We reserved $\mathbf{p}$ for vectors
of unbolded $p,k,q$ which themselves may be $D$-dimensional vectors.

###  IV.2 Group I: All distinct momentum modes : $\mathbf{p_{1}\neq
p_{2},k\neq 0}$.

We take Eq.22 with Hermitian terms grouped as

$\begin{gathered}H_{\lambda}=\frac{\lambda}{4!|\Omega|^{3}}\sum_{\\{\mathbf{p}_{i}\\}}\frac{1}{4\sqrt{\omega_{\mathbf{p}_{1}}\omega_{\mathbf{p}_{2}}\omega_{\mathbf{p}_{3}}\omega_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}}}\left\\{\left(a_{\mathbf{p}_{1}}+a_{\mathbf{p}_{1}}^{\dagger}\right)\left(a_{\mathbf{p}_{2}}+a_{\mathbf{p}_{2}}^{\dagger}\right)\left(a_{\mathbf{p}_{3}}+a_{\mathbf{p}_{3}}^{\dagger}\right)\left(a_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}+a_{-(\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{3})}^{\dagger}\right)\right\\}\end{gathered}$
(53)

This state divides into a tensor product of terms of the following form.

$\displaystyle:H_{\theta}:=\frac{\lambda}{24|\Omega|}\sum_{\mathbf{p}\in
S_{4\mathbf{p}}}\prod_{p_{i}\in\mathbf{p}}\frac{1}{\sqrt{2w_{p_{i}}}}\left(a_{p_{i}}+a_{p_{i}}^{\dagger}\right)$
(54)

We denote the sum of the terms having distinct momentum modes by
$H_{1\varphi}$. We note that this Hamiltonian is trivially normally order
because distinct momentum operators commute. After the qubit mapping using Eq.
46, we get the following.

$\displaystyle H_{1\varphi}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{24|\Omega|^{3}}\sum_{\mathbf{p}\in
S_{4\mathbf{p}}}\prod_{p_{i}\in\mathbf{p}}\frac{1}{\sqrt{2w_{p_{i}}}}\left(\sum_{n}\frac{1}{2}\sqrt{n+1}\left(X_{p_{i},n}X_{p_{i},n+1}+Y_{p_{i},n}Y_{p_{i},n+1}\right)\right)$
(55) $\displaystyle=$ $\displaystyle\frac{\lambda}{96\cdot
16|\Omega|^{3}}\sum_{\mathbf{p}\in S_{4\mathbf{p}}}\sum_{\vec{n}\in
S_{4\vec{n}}}\prod_{(p_{j},n_{j})\in(\mathbf{p},\vec{n})}\sqrt{\frac{n_{j}+1}{w_{p_{j}}}}\left(X_{p_{j},n_{j}}X_{p_{j},n_{j}+1}+Y_{p_{j},n_{j}}Y_{p_{j},n_{j}+1}\right)$

######  Lemma 5.

It is possible to implement $e^{-iH_{1\varphi}t}$ using at most
$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{48}$ $cR_{z}$,
$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{4}$ T,
$\frac{11N^{4}|\Omega|^{2}(|\Omega|-1)}{16}$ CNOT and
$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{24}$ H gates per time step $t$.

The proof of the above lemma, including the detail of the circuit design has
been given in Appendix A.

### Group II : Two distinct momentum modes : $\mathbf{p_{1}=p_{2},k\neq 0}$.

Let $H_{2\varphi}$ is the sum of the terms with momentum modes satisfying the
given constraint. Then,

$\displaystyle:H_{2\varphi}:$ $\displaystyle=$
$\displaystyle\frac{\lambda}{24|\Omega|^{2}}\sum_{p,k}\frac{1}{4\omega_{p}\sqrt{\omega_{p+k}\omega_{p-k}}}\left((a_{p}+a^{\dagger}_{p})^{2}(a_{p+k}+a^{\dagger}_{p+k})(a_{p-k}+a^{\dagger}_{p-k})\right)$
$\displaystyle=$
$\displaystyle\frac{\lambda}{96}\sum_{p,k}\frac{1}{\omega_{p}\sqrt{\omega_{p+k}\omega_{p-k}}}\left((a^{\dagger}_{p})^{2}+(a_{p})^{2}+2n_{p}\right)\left(a_{p+k}+a^{\dagger}_{p+k}\right)\left(a_{p-k}+a^{\dagger}_{p-k}\right),$

and after using after using Eq. 46, Lemma 50, Lemma 3 we get the following.

$\displaystyle:H_{2\varphi}:$ $\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|^{2}}\sum_{p,k}\frac{1}{\omega_{p}\sqrt{\omega_{p+k}\omega_{p-k}}}\left(\left(\sum_{n}\frac{\sqrt{(n+2)(n+1)}}{2}(X_{n}X_{n+2}+Y_{n}Y_{n+2})_{p}\right)+\left(\sum_{n}n(I_{n}-Z_{n})_{p}\right)\right)$
(56)
$\displaystyle\left(\sum_{n}\frac{\sqrt{n+1}}{2}(X_{n}X_{n+1}+Y_{n}Y_{n+1})_{p+k}\right)\left(\sum_{n}\frac{\sqrt{n+1}}{2}(X_{n}X_{n+1}+Y_{n}Y_{n+1})_{p-k}\right)$
$\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|^{2}}\sum_{p,k}\frac{1}{\omega_{p}\sqrt{\omega_{p+k}\omega_{p-k}}}$
$\displaystyle\quad\times\left(\sum_{n_{1},n_{2},n_{3}}c_{n}^{(1)}(X_{n_{1}}X_{n_{1}+2}+Y_{n_{1}}Y_{n_{1}+2})_{p}(X_{n_{2}}X_{n_{2}+1}+Y_{n_{2}}Y_{n_{2}+1})_{p+k}(X_{n_{3}}X_{n_{3}+1}+Y_{n_{3}}Y_{n_{3}+1})_{p-k}\right.$
$\displaystyle\left.\qquad+\sum_{n_{1},n_{2},n_{3}}c_{n}^{(2)}(I_{n_{1}}-Z_{n_{1}})_{p}(X_{n_{2}}X_{n_{2}+1}+Y_{n_{2}}Y_{n_{2}+1})_{p+k}(X_{n_{3}}X_{n_{3}+1}+Y_{n_{3}}Y_{n_{3}+1})_{p-k}\right)$

where $c_{n}^{(1)}=\frac{\sqrt{(n_{1}+2)(n_{1}+1)(n_{2}+1)(n_{3}+1)}}{8}$,
$c_{n}^{(2)}=\frac{n_{1}\sqrt{(n_{2}+1)(n_{3}+1)}}{4}$.

######  Lemma 6.

It is possible to implement $e^{-iH_{2\varphi}t}$ using at most
$\frac{N^{3}|\Omega|^{2}}{3}$ $cR_{z}$, $\frac{8N^{3}|\Omega|^{2}}{3}$ T,
$\frac{20N^{3}|\Omega|^{2}}{3}$ CNOT and $\frac{2N^{3}|\Omega|^{2}}{3}$ H
gates per time step $t$.

The proof of the above lemma, including the detail of the circuit design has
been given in Appendix A.

#### Group III : Two distinct momentum modes : $\mathbf{p_{1}\neq p_{2},k=0}$.

We use $H_{3\varphi}$ to denote the sum of the terms with momentum mode
satisfying the given constraint. Then,

$\displaystyle:H_{3\varphi}:$ $\displaystyle=$
$\displaystyle\frac{\lambda}{24|\Omega|^{2}}\sum_{p_{1},p_{2}}\frac{1}{4\omega_{p_{1}}\omega_{p_{2}}}\left(a_{p_{1}}+a^{\dagger}_{p_{1}}\right)^{2}\left(a_{p_{2}}+a^{\dagger}_{p_{2}}\right)^{2}$
(57) $\displaystyle=$
$\displaystyle\frac{\lambda}{96}\sum_{p_{1},p_{2}}\frac{1}{\omega_{p_{1}}\omega_{p_{2}}}\left((a^{\dagger}_{p_{1}})^{2}+(a_{p_{1}})^{2}+2\hat{n}_{p_{1}}\right)\left((a^{\dagger}_{p_{2}})^{2}+(a_{p_{2}})^{2}+2\hat{n}_{p_{2}}\right).$

After applying Eq. 46, Lemma 50, Lemma 3 we obtain the following.

$\displaystyle:H_{3\varphi}:$ (58) $\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|^{2}}\sum_{p_{1},p_{2}}\frac{1}{\omega_{p_{1}}\omega_{p_{2}}}\left(\sum_{n}\frac{\sqrt{(n+1)(n+2)}}{2}(X_{n}X_{n+2}+Y_{n}Y_{n+2})_{p_{1}}+\sum_{n}n(I_{n}-Z_{n})_{p_{1}}\right)$
$\displaystyle\left(\sum_{n}\frac{\sqrt{(n+1)(n+2)}}{2}(X_{n}X_{n+2}+Y_{n}Y_{n+2})_{p_{2}}+\sum_{n}n(I_{n}-Z_{n})_{p_{2}}\right)$
$\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|^{2}}\sum_{p_{1},p_{2}}\frac{1}{\omega_{p_{1}}\omega_{p_{2}}}\left(\sum_{n_{1},n_{2}}c_{n}^{(3)}(X_{n_{1}}X_{n_{1}+2}+Y_{n_{1}}Y_{n_{1}+2})_{p_{1}}(X_{n_{2}}X_{n_{2}+2}+Y_{n_{2}}Y_{n_{2}+2})_{p_{2}}\right.$
$\displaystyle+\sum_{n_{1},n_{2}}c_{n}^{(4)}(X_{n_{1}}X_{n_{1}+2}+Y_{n_{1}}Y_{n_{1}+2})_{p_{1}}(I_{n_{2}}-Z_{n_{2}})_{p_{2}}+\sum_{n_{1},n_{2}}c_{n}^{(5)}(X_{n_{2}}X_{n_{2}+2}+Y_{n_{2}}Y_{n_{2}+2})_{p_{2}}(I_{n_{1}}-Z_{n_{1}})_{p_{1}}$
$\displaystyle\left.+\sum_{n_{1},n_{2}}n_{1}n_{2}(I_{n}-Z_{p_{1},n_{1}}-Z_{p_{2},n_{2}}+Z_{p_{1},n_{1}}Z_{p_{2},n_{2}})\right)$

where $c_{n}^{(3)}=\frac{\sqrt{(n_{1}+2)(n_{1}+1)(n_{2}+2)(n_{2}+1)}}{4}$,
$c_{n}^{(4)}=\frac{n_{2}\sqrt{(n_{1}+2)(n_{1}+1)}}{2}$,
$c_{n}^{(5)}=\frac{n_{1}\sqrt{(n_{2}+2)(n_{2}+1)}}{2}$.

######  Lemma 7.

It is possible to implement $e^{-iH_{3\varphi}t}$ using at most
$2N^{2}|\Omega|(|\Omega|-1)$ $cR_{z}$, $8N^{2}|\Omega|(|\Omega|-1)$ T,
$16N^{2}|\Omega|(|\Omega|-1)$ CNOT and $3N^{2}|\Omega|(|\Omega|-1)$ H gates
per time step $t$.

The proof of the above lemma, including the detail of the circuit design has
been given in Appendix A.

#### Group IV: All equal momentum modes : $\mathbf{p_{1}=p_{2},k=0}$.

The sum of the terms with all four equal momentum mode is denoted by
$H_{4\varphi}$. These terms must be normal ordered to be

$\displaystyle:H_{4\varphi}:$ $\displaystyle=$
$\displaystyle\frac{\lambda}{24|\Omega|}\sum_{p}\frac{1}{4(\omega_{p})^{2}}\left\\{\left(a^{\dagger}_{p}+a_{p}\right)^{4}\right\\}$
$\displaystyle=$
$\displaystyle\frac{\lambda}{24|\Omega|}\sum_{p}\frac{1}{4(\omega_{p})^{2}}\left\\{\left((a^{\dagger}_{p})^{4}+(a_{p})^{4}\right)+4\left((a_{p}^{\dagger})^{3}a_{p}+a^{\dagger}_{p}(a_{p})^{3}\right)+6\left((a^{\dagger}_{p})^{2}(a_{p})^{2}\right)\right\\}$
$\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|}\sum_{p}\frac{1}{(\omega_{p})^{2}}\left(((a^{\dagger}_{p})^{4}+(a_{p})^{4})+4\left((a_{p}^{\dagger})^{2}\hat{n}_{p}+\hat{n}_{p}(a_{p})^{2}\right)+6\left((\hat{n}_{p})^{2}-\hat{n}_{p}\right)\right);$

and after applying Eq. 46, Lemma 50, Lemma 3 we get the following.

$\displaystyle:H_{4\varphi}:$ (59) $\displaystyle=$
$\displaystyle\frac{\lambda}{96|\Omega|}\sum_{p,n}\frac{1}{(\omega_{p})^{2}}\left(\frac{\sqrt{(n+4)(n+3)(n+2)(n+1)}}{2}(X_{n}X_{n+4}+Y_{n}Y_{n+4})_{p}\right.$
$\displaystyle\left.+2n\sqrt{(n+2)(n+1)}(X_{n}X_{n+2}+Y_{n}Y_{n+2})_{p}+3(n^{2}-n)(I_{n}-Z_{n})_{p}\right)$

######  Lemma 8.

It is possible to implement $e^{-iH_{4\varphi}t}$ using at most $3N|\Omega|$
(c)-$R_{z}$, $4N|\Omega|$ CNOT and $4N|\Omega|$ H gates per time step $t$.

The proof of the above lemma, including the details of the circuit design, has
been given in Appendix A. An estimate of the total number of gates required to
implement $e^{-iHt}$ can be obtained by summing the gate costs in Lemma 4-8.
This is summarized in the following theorem, as well as in Table 2.

| $\mathbf{\\#(c)R_{z}}$ |  $\mathbf{\\#}$ T |  $\mathbf{\\#}$ CNOT |  $\mathbf{\\#}$ H |
---|---|---|---|---|---
$\mathbf{e^{-iH_{0}t}}$ | $N|\Omega|$ | - | - | - | Lemma 4
$\mathbf{e^{-iH_{1\varphi}t}}$ | $\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{48}$ | $\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{4}$ | $\frac{11N^{4}|\Omega|^{2}(|\Omega|-1)}{16}$ | $\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{24}$ | Lemma 5
$\mathbf{e^{-iH_{2\varphi}t}}$ | $\frac{N^{3}|\Omega|^{2}}{3}$ | $\frac{8N^{3}|\Omega|^{2}}{3}$ | $\frac{20N^{3}|\Omega|^{2}}{3}$ | $\frac{2N^{3}|\Omega|^{2}}{3}$ | Lemma 6
$\mathbf{e^{-iH_{3\varphi}t}}$ | $2N^{2}|\Omega|(|\Omega|-1)$ | $8N^{2}|\Omega|(|\Omega|-1)$ | $16N^{2}|\Omega|(|\Omega|-1)$ | $3N^{2}|\Omega|(|\Omega|-1)$ | Lemma 7
$\mathbf{e^{-iH_{4\varphi}t}}$ | $3N|\Omega|$ | - | $4N|\Omega|$ | $4N|\Omega|$ | Lemma 8
Table 2: Summary of the number of gates required to implement the exponentials
of the different Hamiltonian partitions in the occupation basis.

######  Theorem 9.

It is possible to implement $e^{-iHt}$ with the following number of gates per
time step.

  1. 1.

$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{48}+\frac{N^{3}|\Omega|^{2}}{3}+2N^{2}|\Omega|(|\Omega|-1)+4N|\Omega|$
(controlled)-$R_{z}$ gates.

  2. 2.

$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{4}+\frac{8N^{3}|\Omega|^{2}}{3}+8N^{2}|\Omega|(|\Omega|-1)$
additional T gates.

  3. 3.

$\frac{11N^{4}|\Omega|^{2}(|\Omega|-1)}{16}+\frac{20N^{3}|\Omega|^{2}}{3}+16N^{2}|\Omega|(|\Omega|-1)+4N|\Omega|$
CNOT gates.

  4. 4.

$\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{24}+\frac{2N^{3}|\Omega|^{2}}{3}+3N^{2}|\Omega|(|\Omega|-1)+4N|\Omega|$
H gates.

###### Proof.

The proof follows from Lemmas 5, 5,6,7, and 8, which are summarized in 2. ∎

###  IV.3 Trotter Error Analysis

In this section we discuss and bound the various errors that occur during
simulation of the occupation basis Hamiltonian. For convenience, we restate
the occupation basis Hamiltonian in the following equation.

$\displaystyle H_{occ}$ $\displaystyle=$ $\displaystyle H_{0}+H_{\lambda}$
(60) $\displaystyle:H_{0}:$ $\displaystyle=$
$\displaystyle\sum_{p}\omega_{p}a_{p}^{\dagger}a_{p}$ (61)
$\displaystyle:H_{\lambda}:$ $\displaystyle=$ $\displaystyle
H_{1\varphi}+H_{2\varphi}+H_{3\varphi}+H_{4\varphi}$ (62)

##### I. Trotter Error :

One primary source of error is the one inherent in the Trotterization
procedure, where we express the exponential of sum of operators as product of
the exponentiated operators. This is true if and only if the operators are
commuting. If there are non-commuting parts then the resultant error in the
approximation is referred to as the _Trotter error_. Quite a few bounds on the
Trotter error have been derived before [64, 65, 66], but we use the one in
[67], which shows the dependence on nested commutators.

If a Hamiltonian $H=\sum_{\gamma=1}^{\Gamma}H_{\gamma}$ is a sum of $\Gamma$
fragment Hamiltonians, then $e^{-i\tau H}$ can be approximated by product of
exponentials, using the $p^{th}$ order Trotter-Suzuki formula [68],
$\mathscr{S}_{p}(\tau)=e^{-i\tau H}+\mathscr{A}(\tau)$, where
$\|\mathscr{A}(\tau)\|\in O\left(\widetilde{\alpha}_{comm}\tau^{p+1}\right)$
if each $H_{\gamma}$ are Hermitian [67]. Here
$\widetilde{\alpha}_{comm}=\sum_{\gamma_{1},\gamma_{2},\ldots,\gamma_{p+1}=1}^{\Gamma}\|[H_{\gamma_{p+1}},\ldots[H_{\gamma_{2}},H_{\gamma_{1}}]]\|$.
In most applications, it is quite cumbersome, if not impossible, to derive a
rigorous analytic expression of the nested commutators, in order to tightly
bound higher order Trotter error. We use the following bound from [27] in this
work.

$\displaystyle\widetilde{\alpha}_{comm}\leq
2^{p-(p^{\prime}+1)}\sum_{\gamma_{i_{1}},\gamma_{i_{2}},\ldots,\gamma_{i_{p^{\prime}+1}}}\|[H_{\gamma_{p^{\prime}+1}},[\ldots[H_{\gamma_{3}},[H_{\gamma_{2}},H_{\gamma_{1}}]]\ldots]]\|\left(\sum_{\gamma=1}^{\Gamma}\|H_{\gamma}\|\right)^{p-p^{\prime}}\qquad[1\leq
p^{\prime}\leq p]$ (63)

The above formula is specifically useful in scenarios where we can compute
tight bounds till level $p^{\prime}<p$ of nesting. Then we combine the tighter
bound till level $p^{\prime}$ and a looser sum of norm bound for the remaining
levels. In this paper we consider the $2^{nd}$ order Trotter error in our
implementations. In Appendix B we derive rigorous expressions and bounds on
the first level commutators and norms. Then we use the above equation in order
to bound the higher order errors.

We primarily focus on these low level Trotter errors for two reasons. The
symmetric low-level Trotter errors provide tighter bounds [67]. They also
appear to provide more accurate results when compared to explicit calculations
[69].

######  Lemma 10.

Let $H$ be the occupation basis Hamiltonian derived in Eqs. 52-59. Let
$\mathscr{S}_{2}(\tau)$ be a $2^{\rm nd}$ order Trotter-Suzuki approximation
for $e^{-i\tau H}$. Then,

$\displaystyle\left\|e^{-iH\tau}-\mathscr{S}_{p}(\tau)\right\|\in
O\left(\widetilde{\alpha}_{comm}\tau^{3}\right),$ (64)

where

$\displaystyle\widetilde{\alpha}_{comm}$ $\displaystyle\in
O\left(\frac{\lambda^{2}N^{6}}{M^{5}}\right)$ . (65)

The details of this proof are provided in Appendix B.

###  IV.4 Total T-gate cost estimate

Here we review the various contributions to the error in phase estimation.
Based on previous simulation work, the most efficent distribution of error
seems a highly precise phase estimation calculation to offset a single Trotter
step [6]. In particular, we focus on the scheme of [70], which allows us to
obtain optimal constant factors using Fourier-based phase estimation. We use
an approximate quantum Fourier transform (AQFT) [71] in the preparation. The
number of T-gates required to approximate an $n$-qubit QFT circuit up to error
$\epsilon_{QFT}$ is

$\displaystyle
8n\log_{2}\left(\frac{n}{\epsilon_{QFT}}\right)+\log_{2}\left(\frac{n}{\epsilon_{QFT}}\right)\log_{2}\left(\frac{\log_{2}\left(\frac{n}{\epsilon_{QFT}}\right)}{\epsilon_{QFT}}\right).$
(66)

######  Lemma 11.

There exists $\epsilon_{E}^{\prime}$ such that for all sufficiently small
$\epsilon_{E}\leq\epsilon_{E}^{\prime}$ the total number of implementations of
$\mathscr{S}_{2}(\tau)$ required to estimate an eigenstate within error
$\epsilon_{E}$ is bounded above by

$2^{m}\leq\frac{\pi^{2}\sqrt{\widetilde{\alpha}_{\rm
comm}}}{\epsilon_{E}^{3/2}}.$ (67)

###### Proof.

The variance on the phase measured from QPE is given as

$\epsilon_{V}=\langle\cos((\hat{\theta}-\widetilde{\theta}))\rangle^{2}-1,$
(68)

known as the the Holevo variance, where $\widetilde{\theta}$ is the true phase
and $\hat{\theta}$ is the measured outcome. This allows simple analytic
results for a compact phase that may otherwise have artificially high variance
when peaked near the boundary. This can be computed using the phase estimation
algorithm implemented in [70]. The estimated phase can be written as

$\displaystyle\theta_{\mathrm{est}}=\theta_{\mathrm{true}}+\theta+\epsilon_{\mathrm{prep}},$
(69)

where $\theta$ is a random variable with zero mean and Holevo variance
$\epsilon_{V}$ describing the output of phase estimation and
$\epsilon_{\mathrm{prep}}$ represents the systematic errors in the phase that
arise because of approximate gate synthesis, approximate QFT and approximation
due to Trotterization. In the limit of small variance, with high probability
we have the following.

$\displaystyle\epsilon_{\theta}$
$\displaystyle=\sqrt{\mathbb{E}\left[\left(\theta_{\mathrm{est}}-\theta_{\mathrm{true}}\right)^{2}\right]}\approx\sqrt{\epsilon_{V}(\theta)+\left(\pi\epsilon_{\mathrm{QFT}}+\epsilon_{\mathrm{Trotter}}+\epsilon_{\rm
synth}\right)^{2}}$ (70)
$\displaystyle\approx\sqrt{\left(\frac{\pi}{2^{m+1}}\right)^{2}+\left(\pi\epsilon_{\mathrm{QFT}}+\epsilon_{\rm
synth}+\epsilon_{\mathrm{Trotter}}\right)^{2}},$

We select the following in order to get the above expression.

$\displaystyle\epsilon_{\rm
Trotter}\leq\frac{\sqrt{2}\epsilon_{\theta}}{4},\qquad\epsilon_{\rm
QFT}\leq\frac{\sqrt{2}\epsilon_{\theta}}{8\pi},\qquad\epsilon_{\rm
synth}\leq\frac{\sqrt{2}\epsilon_{\theta}}{8},\qquad
2^{m}\geq\frac{\pi}{\sqrt{2}\epsilon_{\theta}}$ (71)

We require that the error in the Trotter-Suzuki expansion to be at most
$\epsilon_{\rm Trotter}=\tilde{\alpha}_{\rm comm}\tau^{3}$. Further, let

$\displaystyle\epsilon_{\theta}=\epsilon_{E}\tau,$ (72)

where $\epsilon_{E}$ is the estimate in the energy that comes from rescaling
the estimate of the phase by a factor of $\tau$. Thus it suffices to choose

$\displaystyle\tau=\sqrt{\frac{\epsilon_{E}}{2^{3/2}\widetilde{\alpha}_{\rm
comm}}}$ (73)

and therefore, the error in the phase is related to the error in the energy
estimate as

$\epsilon_{\theta}=\sqrt{\frac{\epsilon_{E}^{3}}{2^{3/2}\widetilde{\alpha}_{\rm
comm}}}$ (74)

and therefore for $\pi\sqrt{\widetilde{\alpha}_{\rm
comm}}/\epsilon_{E}^{3/2}\geq\frac{1}{1-2^{-1/8}}$ at most

$\displaystyle 2^{m}=\left\lceil\frac{\pi\sqrt{\widetilde{\alpha}_{\rm
comm}}}{2^{1/8}\epsilon_{E}^{3/2}}\right\rceil\leq\frac{\pi\sqrt{\widetilde{\alpha}_{\rm
comm}}}{\epsilon_{E}^{3/2}}.$ (75)

Finally, as the analysis in [70] is tight in the limit as
$\epsilon_{E}\rightarrow 0$ let us assume that $\epsilon_{E}$ is small enough
so that

$\displaystyle\left|\sqrt{\mathbb{E}\left[\left(\theta_{\mathrm{est}}-\theta_{\mathrm{true}}\right)^{2}\right]}-\sqrt{\left(\frac{\pi}{2^{m+1}}\right)^{2}+\left(\pi\epsilon_{\mathrm{QFT}}+\epsilon_{\mathrm{Trotter}}+\epsilon_{\rm
synth}\right)^{2}}\right|\leq\epsilon_{\theta}.$ (76)

We then can ensure a sufficient value of $m$ by taking
$\epsilon_{E}\rightarrow\epsilon_{E}/2$ and using the remaining error budget
to acommodate the error due to $m$ being finite. Using the observation that
since $2^{3/2}<\pi$ it suffices to take

$\displaystyle\frac{\pi\sqrt{\widetilde{\alpha}_{\rm
comm}}}{(\epsilon_{E}/2)^{3/2}}\leq\frac{\pi^{2}\sqrt{\widetilde{\alpha}_{\rm
comm}}}{\epsilon_{E}^{3/2}}=2^{m}$ (77)

And thus this sufficient value is an upper bound on the necessary value.

∎

######  Theorem 12.

Given an eigenstate $\ket{\psi}$ of $H$ such that $H\ket{\psi}=E\ket{\psi}$,
the occupation basis Hamiltonian stated in (62), then it there exists a
quantum algorithm that outputs with probability greater than $2/3$ a value
$\hat{E}$ such that $|\hat{E}-E|\leq\epsilon_{E}$, using a number of $T$ gates
that scales as

$\mathcal{O}\left(\frac{\lambda
N^{7}|\Omega|^{3}}{M^{5/2}\epsilon_{E}^{3/2}}\log(\frac{N|\Omega|}{M\epsilon_{E}})\right)$

T gates and $\mathcal{O}\left(N|\Omega|\right)$ qubits, plus an additional
$O\left(\log{\left(\frac{\lambda N}{M\epsilon_{E}}\right)}\right)$ ancillary
qubits required for phase estimation. Here $M$ is the particle mass for the
field, the log argumenent comes from $1/\epsilon_{r}$
$\epsilon_{r}=\frac{\sqrt{2}\epsilon_{E}}{8N_{r}}\sqrt{\frac{\epsilon_{E}}{2^{3/2}\tilde{\alpha}_{comm}}}$,
$N_{r}\in\mathcal{O}\left(N^{4}|\Omega|^{3}\right)$ and
$\tilde{\alpha}_{comm}\in\mathcal{O}\left(\frac{\lambda^{2}N^{6}}{M^{5}}\right)$.

###### Proof.

Suppose we allocate $\epsilon_{r}$ as an upper bound on the permitted
synthesis error per $R_{z}$ gate. From Theorem 9 and Table 2 we find that the
total number of (controlled)-$R_{z}$ required for each Trotter step is at most

$\displaystyle
N_{r}:=\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{48}+\frac{N^{3}|\Omega|^{2}}{3}+2N^{2}|\Omega|(|\Omega|-1)+4N|\Omega|$
(78)

and so using the T-count estimate in [32] and assuming the T-count of
controlled-$R_{z}$ is at most the T-count of $R_{z}$ [36], the expected number
of T-gates from rotations is at most

$\displaystyle N_{T/R_{z}}\leq
N_{r}\left(3.067\log_{2}(2/\epsilon_{r})-4.327\right).$

Also, from Theorem 9 we require the following additional number of T-gates.

$\displaystyle
N_{+}:=\frac{N^{4}|\Omega|^{2}(|\Omega|-1)}{4}+\frac{8N^{3}|\Omega|^{2}}{3}+8N^{2}|\Omega|(|\Omega|-1).$
(79)

So, the total number of T-gates per Trotter step is

$\displaystyle\mathcal{G}_{T}$ $\displaystyle\leq$ $\displaystyle
N^{4}|\Omega|^{2}(|\Omega|-1)[0.064\log_{2}(2/\epsilon_{r})+0.16]+N^{3}|\Omega|^{2}[1.022\log_{2}(2/\epsilon_{r})+1.224]$
$\displaystyle+N^{2}|\Omega|(|\Omega|-1)[6.134\log_{2}(2/\epsilon_{r})-0.654]+4N|\Omega|[3.067\log_{2}(2/\epsilon_{r})-4.327]$
$\displaystyle\in$ $\displaystyle
O\left(N^{4}|\Omega|^{3}\log_{2}(2/\epsilon_{r})\right).$ (81)

Using Lemmas 10 and 11, the total number of T-gates required to achieve an
eigenstate within total error $\epsilon_{E}$ is at most

$\displaystyle\mathcal{G}_{T}\cdot 2^{m}$ $\displaystyle\leq$
$\displaystyle\frac{\pi^{2}\sqrt{\alpha_{\rm
comm}}}{\epsilon_{E}^{3/2}}\left(N^{4}|\Omega|^{2}(|\Omega|-1)[0.064\log_{2}(2/\epsilon_{r})+0.16]+N^{3}|\Omega|^{2}[1.022\log_{2}(2/\epsilon_{r})+1.224]\right.$
(82)
$\displaystyle\left.+N^{2}|\Omega|(|\Omega|-1)[6.134\log_{2}(2/\epsilon_{r})-0.654]+4N|\Omega|[3.067\log_{2}(2/\epsilon_{r})-4.327]\right)$
$\displaystyle\in$ $\displaystyle O\left(\frac{\lambda
N^{7}|\Omega|^{3}}{m^{5/2}\epsilon_{E}^{3/2}}\log(1/\epsilon_{r})\right).$

The total synthesis error due to approximation of the rotation gates is
$\epsilon_{synth}=\epsilon_{r}\cdot N_{r}$, where the value of
$\epsilon_{synth}$ is given in Eq. 71. Plugging in the values of the timestep
$\tau$ and $\epsilon_{\theta}$ from Eq. (73) and (72), respectively we obtain
the following bound on $\epsilon_{r}$ in order to ensure that the final error
in the estimate of the energy is at most $\epsilon_{E}$.

$\epsilon_{r}\leq\frac{\sqrt{2}\epsilon_{E}\tau}{8N_{r}}=\frac{\sqrt{2}\epsilon_{E}}{8N_{r}}\sqrt{\frac{\epsilon_{E}}{2^{3/2}\tilde{\alpha}_{comm}}}$
(83)

To obtain the total number of ancillary qubits used, we get

$\displaystyle m$ $\displaystyle\in$ $\displaystyle
O\left(\log{\left(\frac{\alpha_{comm}^{1/2}}{\epsilon_{E}^{3/2}}\right)}\right)$
(84) $\displaystyle\in$ $\displaystyle O\left(\log{\left(\frac{\lambda
N}{M\epsilon_{E}}\right)}\right)$

where we have repeatedly used $\log{(A^{m}/B^{n})}\in
O\left(\log{(A/B)}\right)$ for constant $m,n>0$. Thus the theorem is proved. ∎

###  IV.5 Field Amplitude Basis

In this section we describe algorithms to simulate the scalar $\phi^{4}$
Hamiltonian $H_{amp}$ expressed in the amplitude basis, as given in Eq. 11. We
use two types of algorithms - qubitization [72, 73, 21] and Trotterization
[68]. For the former we describe three approaches, which mainly differ in the
LCU (linear combination of unitaries) decomposition of the operators. Similar
approaches have been considered in [25]. In Section IV.5.1 we discuss a
decomposition of $\Phi$ using an equal weight LCU decomposition of the field
operator (Algorithm I). Next, in Section IV.5.2 we describe another LCU
decomposition of $\Phi$ as sum of mainly Z-operators, which not only helps us
in developing a qubitization based algorithm (Algorithm IIIa) in IV.5.3, but
it also makes the expression amenable to Trotterization (Algorithm II), as
discussed in Section IV.5.2. In Section IV.5.4 we describe another more
compact LCU decomposition of all operators using binary represenation of
integers. With this we describe another qubitization algorithm (Algorithm
IIIb). Before describing our algorithms we mention the following results which
have been used in the LCU based approaches in order to reduce the gate
complexity.

##### Recursive block encoding :

We use the following theorem to recursively block encode $H_{amp}$ using a
divide and conquer approach, as described in [27], where it has been shown
that with such an approach it is possible to block encode with a smaller
number of gates. Suppose without loss of generality, we have a Hamiltonian
$H_{i}$ expressed as a linear combination of unitaries (LCU) i.e.
$H_{i}=\sum_{j=1}^{M_{i}}h_{ij}U_{ij}$, such that
$\lambda_{i}=\sum_{j}|h_{ij}|$. In this case, we can have a
$\left(\lambda_{i},\log M_{i},0\right)$-block encoding of $H_{i}$ using an
ancilla preparation subroutine and a unitary selection subroutine, which we
denote by $\text{PREP}_{i}$ and $\text{SELECT}_{i}$ respectively.

$\displaystyle\text{PREP}_{i}\ket{0}^{\log M_{i}}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{M_{i}}\sqrt{\frac{h_{ij}}{\lambda_{i}}}\ket{j}$ (85)
$\displaystyle\text{SELECT}_{i}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{M_{i}}\ket{j}\bra{j}\otimes U_{ij}$ (86)

It can be shown that [74]

$\displaystyle\bra{0}\text{PREP}_{i}^{\dagger}\cdot\text{SELECT}_{i}\cdot\text{PREP}_{i}\ket{0}$
$\displaystyle=$ $\displaystyle\frac{H_{i}}{\lambda_{i}}.$ (87)

Suppose we have $M$ Hamiltonians denoted by $H_{1},\ldots,H_{M}$, each of
which has an LCU decomposition and for each one of them we define the
subroutines as in Eq. 85 and 86. Now we use these subroutines to define the
following,

$\displaystyle\text{PREP}\ket{0}^{\log M+\sum_{i}\log M_{i}}$ $\displaystyle=$
$\displaystyle\left(\sum_{i=1}^{M}\sqrt{\frac{w_{i}\lambda_{i}}{\mathcal{\mathcal{A}}}}\ket{i}\right)\otimes\bigotimes_{i=1}^{M}\text{PREP}_{i}$
(88) SELECT $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}\left(\ket{i}\bra{i}\otimes\bigotimes_{k=1}^{i-1}\mathbb{I}\otimes\text{SELECT}_{i}\otimes\bigotimes_{k=i+1}^{M}\mathbb{I}\right),$
(89)

where $w_{i}>0$ and $\mathcal{A}=\sum_{i=1}^{M}w_{i}\lambda_{i}$. We can use
the above two subroutines to block encode a linear combination of Hamiltonians
as follows.

######  Theorem 13 ([27]).

Let $H=\sum_{i=1}^{M}w_{i}H_{i}$ be the sum of $M$ Hamiltonians and each of
them is expressed as sum of unitaries as :
$H_{i}=\sum_{j=1}^{M_{i}}h_{ij}U_{ij}$ such that
$\lambda_{i}=\sum_{j}|h_{ij}|$, $w_{i}>0$. Each of the summand Hamiltonian is
block-encoded using the subroutines defined in Eq. 85 and 86. Then, we can
have an $(\mathcal{A},\lceil\log_{2}(M)\rceil,0)$-block encoding of $H$, where
$\mathcal{A}=\sum_{i=1}^{M}w_{i}\lambda_{i}$, using the ancilla preparation
subroutine (PREP) defined in Eq. 88 and the unitary selection subroutine
(SELECT) defined in Eq. 89.

  1. 1.

The PREP subroutine has an implementation cost of
$\mathcal{C}_{\text{PREP}}=\sum_{i=1}^{M}\mathcal{C}_{\text{PREP}_{i}}+\mathcal{C}_{w}$,
where $\mathcal{C}_{\text{PREP}_{i}}$ is the number of gates to implement
$\text{PREP}_{i}$ and $\mathcal{C}_{w}$ is the cost of preparing the state
$\sum_{i=1}^{M}\sqrt{\frac{w_{i}\lambda_{i}}{\mathcal{A}}}\ket{i}$.

  2. 2.

The SELECT subroutine can be implemented with a set of multi-controlled-X
gates -
$\\{M_{i}\text{ pairs of }C^{\log_{2}M_{i}+1}X\text{ gates }:i=1,\ldots,M\\}$,
$M$ pairs of $C^{\log M}X$ gates and $\sum_{i=1}^{M}M_{i}$ single-controlled
unitaries - $\\{cU_{ij}:j=1,\ldots,M_{i};i=1,\ldots,M\\}$.

Additionally, suppose in the above theorem, all the $H_{i}$ are the same but
they act on disjoint subspaces. In this case, each $\text{PREP}_{i}$ is the
same and so it is sufficient to keep only one copy of $\text{PREP}_{i}$ in the
PREP subroutine of Eq. 88. We can absorb $w_{i}$ in the weights of the
unitaries obtained in the LCU decomposition of $H_{i}$. Thus, in this case we
have

$\displaystyle\text{PREP}\ket{0}^{\log M+\log
M_{i}}=\left(\sqrt{\frac{1}{M}}\sum_{i=1}^{M}\ket{i}\right)\otimes\text{PREP}_{i}.$
(90)

We require only $\lceil\log M\rceil$ H gates to prepare the superposition in
the first register by padding out the number of such subspaces to be a power
of $2$. We also need to make slight modifications in the SELECT procedure.
This time, we keep an extra ancilla qubit, initialized to 0, in each subspace.
Given a particular state of the first register, we select a subspace by
flipping the qubit in the corresponding subspace. The unitaries in each
subspace are now additionally controlled on this qubit (of its own subspace).

##### Group of multi-controlled-X gate synthesis :

We can further optimize the number of gates by implementing the group of
multi-controlled-unitaries in the SELECT subroutines, using the following
theorem [27]. Here we partition the control qubits into different groups,
store intermediate information in some ancillae and then implement the
required logic using these intermediate results.

######  Theorem 14 ([27]).

Consider the unitary $U=\sum_{j=0}^{M-1}\outerproduct{j}{j}\otimes U_{j}$ for
unitary operators $U_{j}$ that can be implemented controllably. We assume $M$
is a power of 2 for simplicity. Suppose we have $\log_{2}M$ qubits and $M$
(compute-uncompute) pairs of $C^{\log_{2}M}X$ gates for selecting the $M$
basis states. Let $r_{1},\ldots,r_{n}\geq 1$ be positive fractions such that
$\sum_{i=1}^{n}\frac{1}{r_{i}}=1$ and $\frac{\log_{2}M}{r_{i}}$ are integers.
Then, $U$ can be implemented with a circuit with

$\sum_{i=1}^{n}M^{\frac{1}{r_{i}}}C^{\frac{\log_{2}M}{r_{i}}}X+MC^{n}X$

(compute-uncompute) pairs of gates, $M$ applications of controlled $U_{j}$ and
at most $\sum_{i=1}^{n}M^{\frac{1}{r_{i}}}$ ancillae.

Following the construction in [75, 76], the number of T-gates required to
implement such multiply controlled gates is

$\displaystyle\mathcal{T}_{n}=\sum_{i=1}^{n}M^{\frac{1}{r_{i}}}\left(\frac{4\log_{2}M}{r_{i}}-4\right)+M(4n-4).$
(91)

With the help of logical AND gadgets we do not require to use any T-gate for
the uncomputation part.

##### Equal weight LCU:

Note that the Hamiltonian Eq. (11) consists of 4 different families of terms,
the $\Pi^{2}$, $\Phi\Phi$, $\Phi^{2}$, and $\Phi^{4}$ terms, each of which
share the same coefficients. In principle, if we can further ensure that the
LCU decompositions of each of those terms provides the same weight to every
unitary in the LCU, then we can exploit this structure to drastically simplify
the PREP and SELECT circuits. To achieve such an equal weight LCU, we will
make use of the following lemmma

######  Lemma 15.

Consider any arbitrary diagonal matrix

$A=\sum_{j=0}^{L-1}n_{j}|j\rangle\langle j|$ (92)

where $n_{j}\in\mathbb{Z}^{+}$. Then,

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=0}^{2n_{max}-1}U^{(i)},\quad\text{where}$
$\displaystyle U_{j,k}^{(i)}$ $\displaystyle=$
$\displaystyle\delta_{jk}\left[2\Theta(n_{max}+n_{j}-i-1)-1\right]$ (93)

where the step-function is defined as

$\Theta(x)=\begin{cases}1,&\text{if $x\geq 0$}\\\ 0,&\text{if
$x<0$}\end{cases}$ (94)

###### Proof.

We simply note that
$n_{j}=\frac{1}{2}\left[(n_{max}+n_{j})-(n_{max}-n_{j})\right]$, we can
describe such a matrix as a sum consisting of unitaries with entries only in
$\pm 1$ such that the $(j,j)$-th entry in the first $(n_{max}+n_{j})$
unitaries equals $+1$, and equals $-1$ in the remaining $(n_{max}-n_{j})$
unitaries. ∎

Note that while the dimensionality of the $U_{j,k}^{i}$ is given by $L$, the
number of terms in the LCU Eq. (93) is independent of $L$.

##### LCU decomposition of integer diagonal matrices :

If a diagonal matrix consists of integers only then we can have an LCU
decomposition consisting of $O(\log m_{max})$ signature matrices, where
$m_{max}$ is the maximum absolute value of any of its entries. This can be
obtained by a binary decomposition of each integer, as stated below.

######  Lemma 16 ([27]).

Let $M_{I}$ is a $N\times N$ diagonal integer matrix, which has $N^{\prime}$
positive integers whose maximum value is $m_{max}^{\prime}$ and
$N^{\prime\prime}$ negative integers such that
$m_{max}^{\prime\prime}=\max_{i}\\{|M_{I}[i,i]|:M_{I}[i,i]<0\\}$. Then,
$M=c_{0}\mathbb{I}+\sum_{i=1}^{N^{\prime}+N^{\prime\prime}}c_{i}D_{i}$, where
$D_{i}$ are signature matrices and
$N^{\prime}\leq\lceil\log_{2}(m_{max}^{\prime}+1)\rceil=\zeta^{\prime}$ and
$N^{\prime\prime}\leq\lceil\log_{2}(m_{max}^{\prime\prime}+1)\rceil=\zeta^{\prime\prime}$.
Also, $\sum_{i=0}^{N^{\prime}+N^{\prime\prime}}|c_{i}|\leq
2^{\zeta^{\prime}}-1$.

Specifically, we can use the following lemma to have an LCU decomposition of
$\frac{\phi}{\Delta\phi}=\text{diag}\left(-k,\ldots,-1,0,1;\ldots,k-1\right)$.

######  Lemma 17 (Lemma 11 in [27]).

Let
$U=\mathbb{I}_{(\zeta)}\otimes\ldots\otimes\mathbb{I}_{(\ell+1)}\otimes\text{Z}_{(\ell)}\otimes\mathbb{I}_{(\ell-1)}\otimes\ldots\otimes\mathbb{I}_{(1)}$
is a tensor product of $\zeta$ single-qubit unitaries where Z is applied on
qubit $\ell$ and $\mathbb{I}$ on the rest. Then $U$ is a diagonal matrix of
the following form.

$\displaystyle U_{j,j}$ $\displaystyle=$ $\displaystyle 1\qquad\text{if
}\qquad j=2^{\ell}k,2^{\ell}k+1,\ldots 2^{\ell}k+2^{\ell-1}-1$
$\displaystyle=$ $\displaystyle-1\qquad\text{if }\qquad
j=2^{\ell}k+2^{\ell-1},2^{\ell}k+2^{\ell-1}+1,\ldots 2^{\ell}k+2^{\ell}-1$
$\displaystyle\qquad\text{where }\qquad k=0,1,\ldots,2^{\zeta-\ell}-1$

It follows that

$\displaystyle\frac{\phi}{\Delta\phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j=0}^{\log_{2}k}2^{j}Z_{j}-\frac{1}{2}\mathbb{I}.$
(95)

##### Error in the scalar field :

We also bound the error in the scalar field in terms of the target energy
scale as below.

######  Theorem 18.

Let $|\phi_{max}|$ be the maximum allowed value of the scalar field. Then, it
suffices to take

$|\phi_{max}|=\left(\frac{\epsilon
E_{max}}{C(M,\Lambda,d)|\Omega|}\right)^{1/4}$ (96)

where $E_{max}>0$ is the maximum energy scale we wish to allow in our
simulation of the Hamiltonian $H_{amp}$ given in Eq. (11),
$\epsilon=\Pr(H_{amp}\geq E_{max})$ is the probability that a measurement of
the Hamiltonian exceeds $E_{max}$, $|\Omega|$ is the lattice size, and
$C(M,\Lambda,d)=M^{2}+3d+\frac{\Lambda}{4!}+\frac{3}{2}$.

###### Proof.

Let $|\phi_{max}|$ be the maximum allowed value of the scalar field in our
simulation. Then, in the field amplitude basis defined by
$\hat{\Phi}\ket{\phi}=\phi\ket{\phi}$, it readily follows that
$\bra{\phi}\hat{\Phi}^{n}\ket{\phi}\leq|\phi_{max}|^{n}$ at a single site.
Similarly, since the momentum operator is related to the field operator via a
Fourier transform, $\hat{\Pi}=\mathcal{F}^{\dagger}\hat{\Phi}\mathcal{F}$,
then letting $\ket{\tilde{\phi}}=\mathcal{F}\ket{\phi}$, it also readily
follows for a single lattice site that

$\displaystyle\bra{\phi}\hat{\Pi}^{n}\ket{\phi}$ $\displaystyle=$
$\displaystyle\bra{\tilde{\phi}}\hat{\Phi}^{n}\ket{\tilde{\phi}}$ (97)
$\displaystyle\leq$ $\displaystyle|\phi_{max}|^{n}$

Now, the Hamiltonian is given by

$H_{amp}=\sum_{\vec{x}\in\Omega}\left[\frac{1}{2}\Pi^{2}(x)+(M^{2}+d+1)\Phi^{2}(\vec{x})+\frac{\Lambda}{4!}\Phi^{4}(\vec{x})-2\sum_{i=1}^{d}\Phi(\vec{x})\Phi(\vec{x}+a\hat{x}_{i})\right]$
(98)

It therefore follows that

$\displaystyle\bra{\phi}H_{amp}\ket{\phi}$ $\displaystyle=$
$\displaystyle\sum_{\vec{x}\in\Omega}\bra{\phi}\left[\frac{1}{2}\Pi^{2}(x)+(M^{2}+d+1)\Phi^{2}(\vec{x})+\frac{\Lambda}{4!}\Phi^{4}(\vec{x})-2\sum_{i=1}^{d}\Phi(\vec{x})\Phi(\vec{x}+a\hat{x}_{i})\right]\ket{\phi}$
(99) $\displaystyle\leq$
$\displaystyle\sum_{\vec{x}\in\Omega}\left[\left(\frac{1}{2}+(M^{2}+d+1)+2d\right)|\phi_{max}|^{2}+\frac{\Lambda}{4!}|\phi_{max}|^{2}\right]$
$\displaystyle\leq$
$\displaystyle\sum_{\vec{x}\in\Omega}\left[\left(\frac{1}{2}+M^{2}+3d+1+\frac{\Lambda}{4!}\right)|\phi_{max}|^{2}\right]$
$\displaystyle=$ $\displaystyle C(M,\Lambda,d)|\Omega||\phi_{max}|^{4}$

where $C(M,\Lambda,d)=M^{2}+3d+1+\frac{\Lambda}{4!}+\frac{1}{2}$. Let
$E_{max}>0$ be the maximum allowed energy scale in our simulation. By Markov’s
inequality then,

$\displaystyle\Pr(H_{amp}\geq E_{max})$ $\displaystyle\leq$
$\displaystyle\frac{\bra{\phi}H_{amp}\ket{\phi}}{E_{max}}$ (100)
$\displaystyle\leq$
$\displaystyle\frac{C(M,\Lambda,d)|\Omega||\phi_{max}|^{4}}{E_{max}}$

from which the statement of the theorem readily follows. ∎

####  IV.5.1 Algorithm I : Equal weight LCU

Here, we describe how to construct the PREP and SELECT oracles for the
$\Phi^{4}$ Hamiltonian $H_{amp}$ appearing in Eq. (11), using an equal weight
LCU for the operators appearing in the Hamiltonian. By ensuring that we
provide the same coefficient to each of the unitaries appearing in the LCU
decompositions of $\Phi$, $\Phi^{2}$, $\Phi^{4}$ and $\Pi^{2}$, we essentially
ensure that we have a total of only 4 families of terms each sharing the same
coefficient within that family of terms. This drastically simplifies the PREP
circuit, the only non-trivial part of which is to prepare a 2-qubit state,
while the rest of it is composed of transversal Hadamards. We also describe
how to efficiently carry out a SELECT circuit that does not require selecting
each of the unitaries within a family in succession, but can rather apply all
of them simultaneously using a comparator.

First, we describe the decompositions of the operators appearing in the
Hamiltonian of Eq. (11). We apply Lemma 94 to obtain these decompositions. The
$\Phi\Phi$, $\Phi^{2}$ and $\Phi^{4}$ terms are all built from the same basic
diagonal operator $\frac{\hat{\Phi}}{\Delta\Phi}=Diag(-k+1,\dots,0,\dots,k)$.

##### Decompositions

First, for the $\hat{\Phi}$ operator, we use $L=2k$, $n_{max}=k$, and
$n_{j}=j-k+1$ in Eq. (93) to get

$\displaystyle\frac{\hat{\Phi}}{\Delta\Phi}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=0}^{2k-1}U^{(i)},\quad\text{where}$
$\displaystyle U^{(i)}$ $\displaystyle=$
$\displaystyle-\sum_{j=0}^{i-1}|j\rangle\langle
j|+\sum_{j=i}^{2k-1}|i\rangle\langle
i|=\sum_{j=0}^{2k-1}\left[2\Theta(j-i)-1\right]|j\rangle\langle j|$ (101)

On any two sites $a$ and $b$, the operator $-\hat{\Phi}_{a}\hat{\Phi}_{b}$ is
then simply given by

$\displaystyle-\frac{\hat{\Phi}_{a}\hat{\Phi_{b}}}{\left(\Delta\Phi\right)^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{i_{a},j_{a}=0}^{2k-1}\sum_{i_{b},j_{b}=0}^{2k-1}\left[2\Theta(j_{a}-i_{a})-1\right]\left[1-2\Theta(j_{b}-i_{b})\right]\ket{j_{a}}\bra{j_{a}}\otimes\ket{j_{b}}\bra{j_{b}}$
(102)

For the $\hat{\Phi}^{2}$ term, we use $L=2k^{2}$, $n_{max}=k^{2}$ and
$n_{j}=(k-j-1)^{2}$ in Eq. (93) to get

$\displaystyle\left(\frac{\hat{\Phi}}{\Delta\Phi}\right)^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=0}^{2k^{2}-1}U^{(i)},\quad\text{where}$
$\displaystyle U^{(i)}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{2k-1}\left[2\Theta(k^{2}+(k-j-1)^{2}-i-1)-1\right]|j\rangle\langle
j|$ (103)

Similarly, for the $\hat{\Phi}^{4}$ term, we use $L=2k$, $n_{max}=k^{4}$ and
$n_{j}=(k-j-1)^{4}$ in Eq. (93) to get

$\displaystyle\left(\frac{\hat{\Phi}}{\Delta\Phi}\right)^{4}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=0}^{2k^{4}-1}U^{(i)},\qquad\text{where}$
$\displaystyle U^{(i)}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{2k-1}\left[2\Theta(k^{4}+(k-j-1)^{4}-i-1)-1\right]\ket{j}\bra{j}$
(104)

#####  PREP circuit

As the above decompositions illustrate, we have LCU decompositions for each of
the terms in $H_{amp}$ such that all the unitaries share the same coefficients
within each of 4 different families. In the language of Theorem 13, this means
that we have $M=4$ different families. These equal weight LCU decompositions
lead to a much simpler PREP circuit. In particular, we only need to prepare a
non-trivial $\log_{2}{M}=2$ qubit register whose amplitudes encode the shared
coefficients of the 4 groups of terms in the Hamiltonian (11). The
$\text{PREP}_{i}$ (sub-PREPARE) circuits for each of these groups, which would
control which of the $U^{(i)}$’s in the LCUs above would be selectively
applied, then simply become a layer of transversal Hadamard gates.

In other words, the PREP circuit then simplifies to

$\text{PREP}=\text{PREP}_{F}\otimes H^{\otimes\log_{2}{|\Omega|}}\otimes
H^{\otimes\log_{2}{2k^{4}}}$ (105)

where

$\text{PREP}_{F}\ket{0}^{\otimes
2}=\alpha_{\pi^{2}}\ket{00}+\alpha_{\phi^{2}}\ket{01}+\alpha_{\phi^{4}}\ket{10}+\alpha_{\phi\phi}\ket{11}$
(106)

where the amplitudes $\alpha_{j}$ encode the coefficients of the 4 families of
terms in the Hamiltonian (11)

$\displaystyle|\alpha_{\pi^{2}}|^{2}/\mathcal{N}=\frac{1}{2}$ ,
$\displaystyle|\alpha_{\phi^{2}}|^{2}=\left(M^{2}+d+1\right)/\mathcal{N}$
$\displaystyle|\alpha_{\phi^{4}}|^{2}=\frac{\Lambda}{4!}/\mathcal{N}$ ,
$\displaystyle|\alpha_{\phi\phi}|^{2}=2/\mathcal{N}$ (107)

where
$\mathcal{N}=\sum_{j\in\\{\pi^{2},\phi^{2},\phi^{4},\phi\phi\\}}|\alpha_{j}|^{2}$.
Note that this is not the same as the coefficient 1-norm of the Hamiltonian,
which instead is given by

$|\alpha|_{1}=|\Omega|\left[\frac{k^{4}\Delta_{\phi}^{4}\Lambda}{4!}+k^{2}\Delta_{\phi}^{2}\left(M^{2}+3d+\frac{3}{2}\right)\right]$
(108)

where we take $\Delta_{\phi}=\sqrt{2\pi/2^{n_{Q}}}=\sqrt{\pi/k}$ where
$n_{Q}=\log_{2}{2k}$ is the number of qubits required to implement the field
register at a single lattice site. This choice ensures that the free part of
the Hamiltonian corresponds to the usual eigenspectrum of a harmonic
oscillator.

Thus, the PREP circuit will contribute at most $O(1)$ T gates, by preparing a
superposition of the four non-trivial coefficients in the $\phi^{4}$
Hamiltonian, since the LCUs for all the operators appearing in the Hamiltonian
are all equal weight. We need only prepare a non-trivial 2-qubit state, with
real amplitudes. It was shown by Vatan and Williams [77] that such an $SO(4)$
operator requires at most 12 RZ gates and 2 CNOTs.

We now use the result of [28] that gives us a T gate cost of
$3.067\log_{2}(2/\epsilon_{synth})-4.327$ to compile an RZ gate to target
precision $\epsilon_{synth}$. Other such cost estimates include the result of
[78], which estimates that on average, the T gate cost of a
synthesized/compiled unitary $V$ $\epsilon_{synth}$-approximating a single RZ
gate $U$ (i.e. $|U-V|\leq\epsilon_{synth}$) is
$1.15\log_{2}(1/\epsilon_{synth})$, while [79] estimates a worst case upper
bound of $10+4\log_{2}(1/\epsilon_{synth})$ T gates per RZ gate. Since the
errors are additive, the total T gate cost increases by the number of RZ gates
times the compilation cost of a single RZ gate. Thus, the total T gate cost of
the PREP circuit is given by

$\displaystyle Count(T)_{\text{PREP}}$ $\displaystyle=$ $\displaystyle
36.804\log_{2}(2/\epsilon_{synth})-51.924$ (109)

##### Simultaneous SELECT

In principle, with the simplification of the PREP circuit described above, we
could selectively apply each of the unitaries in a particular LCU by
controlling on the state of some register prepared in equal superposition of
all bitstrings. However, this decomposition also allows us to simultaneously
apply all the $\text{SELECT}_{i}$ circuits corresponding to these unitaries.
In particular, by making use of the comparator CMP operation, defined as

$\texttt{CMP}\ket{i}\ket{j}\ket{0}=\ket{i}\ket{j}\ket{j<i}$ (110)

we can apply the following lemma

######  Lemma 19.

For integers $i$ and $j$, we have

$\texttt{CMP}^{\dagger}\left(\mathbb{I}\otimes\mathbb{I}\otimes
Z\right)\texttt{CMP}\ket{i}\ket{j}\ket{0}=\left(2\Theta(j-i)-1\right)\ket{i}\ket{j}\ket{0}$
(111)

where the step function $\Theta(x)$ is defined in Eq. (94).

from which we readily obtain the following corollary

######  Corollary 20.

Let $\ket{\psi}=\sum_{j=0}^{D-1}\alpha_{j}\ket{j}$ denote the state of the
scalar field at some lattice point, and
$\ket{\overline{+}}=\sum_{i=0}^{D-1}\frac{1}{\sqrt{D}}\ket{i}$ an equal
superposition register. Then,

$\ket{\overline{+}}\otimes\frac{\hat{\Phi}}{\Phi_{max}}\ket{\psi}\otimes\ket{0}=\texttt{CMP}^{\dagger}\left(\mathbb{I}\otimes\mathbb{I}\otimes
Z\right)\texttt{CMP}\ket{\overline{+}}\ket{\psi}\ket{0}$ (112)

where $\Phi_{max}=k\Delta\Phi$.

In this way, we obviate the need to selectively apply each of the unitaries
$U^{(i)}$ appearing in Eq. (101) individually. Instead, we can apply all of
them in superposition by a single application of the comparator and its
inverse. This simplifies the $\text{SELECT}_{i}$ circuits necessary to
implement each of the terms in the Hamiltonian. Fig 3 demonstrates this
circuit equivalence, and the advantage it confers. Applying each of the
equally weighted unitaries individually would require $O(|\Omega|k)$ many
operations controlled on the state of $O(\log{|\Omega|k})$ many qubits,
leading to a $T$ gate count of $O(|\Omega|k\log{|\Omega|k})$. Instead, we have
$O(|\Omega|)$ many single-qubit operations controlled on the state of
$O(\log{|\Omega|})$ many qubits, leading to a significantly reduced $T$ gate
count of $O(|\Omega|\log{|\Omega|})$ for the SELECT operation.

Note that instead of the comparator, we are also free to apply any other
operator $\texttt{CMP}^{\prime}$ of the form

$\displaystyle\texttt{CMP}^{\prime}\ket{i}\ket{j}\ket{scratch}\ket{0}$
$\displaystyle=$ $\displaystyle\ket{\Psi_{i,j,scratch}}\ket{j<i}$ (113)

which may significantly reduce the gate cost. An example of such an operation
would be to compute the high carry bit of $\overline{\overline{j}+i}=j-i$,
which is 1 iff $j<i$, but leaves the state of the $\ket{i}$, $\ket{j}$ and
other ancillary qubits $\ket{scratch}$ entangled. Once the appropriate phase
has been extracted from the $\ket{j<i}$ qubit, we simply run the inverse of
this operation. Noting also that

$\displaystyle\ket{\mathbf{x}}\bra{\mathbf{x}}\otimes\left(U^{\dagger}VU\right)+\ket{\mathbf{x}^{\perp}}\bra{\mathbf{x}^{\perp}}\otimes\mathbb{I}$
$\displaystyle=$ $\displaystyle\left(\mathbb{I}\otimes
U^{\dagger}\right)\left(\ket{\mathbf{x}}\bra{\mathbf{x}}\otimes
V+\ket{\mathbf{x}^{\perp}}\bra{\mathbf{x}^{\perp}}\otimes\mathbb{I}\right)\left(\mathbb{I}\otimes
U\right)$ (114)

this means that the controlled version of the entire operation $H^{\otimes
m}\texttt{CMP}^{\prime\dagger}Z_{anc}\texttt{CMP}^{\prime}H^{\otimes m}$
simply needs to control the $Z_{anc}$ operation. An example circuit involving
the $\texttt{CMP}^{\prime}$ operation is depicted in Fig. 4.

We can now split up the entire SELECT circuit as

$\text{SELECT}=\sum_{i\in\\{\pi^{2},\phi^{2},\phi^{4},\phi\phi\\}}\ket{\alpha_{i}}\bra{\alpha_{i}}\otimes\sum_{a=0}^{|\Omega|-1}\ket{a}\bra{a}\otimes\text{SELECT}_{i,a}$
(115)

where the left-most register controls the family of terms in the Hamiltonian,
the middle register controls which site we apply an operator to, and finally
the $\text{SELECT}_{i}$ circuits applies the relevant transformations for each
of the family of terms as described below, and further detailed in Appendix
C.1.2.

Figure 3: Circuit equivalence demonstrating how the comparator CMP can
implement the entire sum of unitaries in Eq. (101). The “family” register
controls which of the 4 families of terms in the Hamiltonian we want to apply,
the “site” register controls which of the lattice sites we want to apply the
operator to, and the “LCU term” register controls which of the unitaries of
Eq. (101) we wish to apply. Figure 4: An example circuit for $k=4$
implementing the operation
$\texttt{CMP}^{\prime\dagger}Z_{anc}\texttt{CMP}^{\prime}$ using the logical
AND construction of [76]. The field register
$\ket{\phi}=\ket{a_{2}}\ket{a_{1}}\ket{a_{0}}$ is compared against an equal
superposition of bitstrings using $\texttt{CMP}^{\prime}$, the comparison
recorded onto an ancilla using $Z_{anc}$, followed by the uncomputation
$\texttt{CMP}^{\prime\dagger}$ which involves measurements and post-
measurement Clifford operations.

#####  $\text{SELECT}_{\phi\phi}$ term

To implement the $-\Phi\Phi$ term, we can simply make use of Corollary (20)
and the fact that $XZX=-Z$ to obtain

$\ket{\overline{+}}\otimes\left(-\frac{\hat{\Phi}}{\Phi_{max}}\right)\ket{\psi}\otimes\ket{0}=\texttt{CMP}^{\dagger}\left(\mathbb{I}\otimes\mathbb{I}\otimes
XZX\right)\texttt{CMP}\ket{\overline{+}}\ket{\psi}\ket{0}$ (116)

Then, for any two lattice sites $a$ and $b$, the operator
$\left(-\hat{\Phi}_{a}\hat{\Phi}_{b}\right)/(\Delta\Phi)^{2}$ is then
implemented by the product of operations given in Eqs. (20) and (116). In $D$
spatial dimensions, a lattice site has $2D$ neighbors. Since the operation
$\left(\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}\right)$ squares to the
identity, this means we only need to apply the operation
$\left(\texttt{CMP}^{\dagger}X_{anc}Z_{anc}X_{anc}\texttt{CMP}\right)$ on each
of the neighbors of a given lattice site, without applying a similar operation
on the given lattice site itself. Thus, the entire sub-select unitary for this
family of terms is then given by

$\displaystyle\text{SELECT}_{\phi\phi,a}$ $\displaystyle=$
$\displaystyle\sum_{b=0}^{2D-1}\ket{b}\bra{b}\otimes\left(\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}\right)_{a}\left(\texttt{CMP}^{\dagger}X_{anc}Z_{anc}X_{anc}\texttt{CMP}\right)_{b}$
(117)

where the $\ket{b}$ register, consisting of $\log_{2}{2D}$ qubits, controls
which of the neighbors we apply the operation to. For notational simplicity,
in the expression above, the subscripts $a$ and $b$ denote not just the
respective lattice sites, but also the associated ancilla required to
implement the CMP, as well as the $X_{anc}$ and $Z_{anc}$ operations. As
described earlier, the only operation inside the parantheses above that needs
to be controlled is the $Z_{anc}$ operation, since all other operations on
either side of it are self-conjugate. Any given lattice site will receive one
application of $\left(\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}\right)$ and
$2D$ applications of
$\left(\texttt{CMP}^{\dagger}X_{anc}Z_{anc}X_{anc}\texttt{CMP}\right)$ coming
from its neighbors. Between these successive applications, CMP and its inverse
multiply to the identity, as do several applications of the $X$ operator.
Therefore, per lattice site, we require one application each of CMP and its
conjugate, two applications of the $X_{anc}$ operator, and $2D+1$ applications
of a controlled $Z_{anc}$ operator, each of which is controlled on the state
of $\log_{2}{(8D|\Omega|)}$ ($=2+\log_{2}{|\Omega|}+\log_{2}{(2D)}$) many
qubits, where the factor of 2 comes from the PREP circuit, the factor of
$\log_{2}{|\Omega|}$ from the lattice size, and $\log_{2}{(2D)}$ from the
number of neighboring lattice sites in $D$ spatial dimensions.

A $C^{n}$-$Z$ gate can be implemented [80, 76] with $n-2$ ancillas initialized
in the $\ket{0}$ state, using $n-1$ logical AND computations (and just as many
uncomputations) along with 2 Hadamard gates. Each logical AND computation
requires 4 T gates. Therefore, the controlled $Z_{anc}$ operations naively
require $O(D|\Omega|\log_{2}{(D|\Omega|)})$ T gates across the entire lattice,
and in each of the $D$ spatial dimensions.

A more efficient approach using unary iteration is also described in [70],
which brings the cost for implementing such a SELECT operation consisting of
$L$ operations down to $4L-4$ T gates, instead of $O(L\log_{2}{L})$. Briefly,
this method exploits the fact that successive control patterns differ in at
most a constant number of bits, so that the collective circuit can be
simplified. In our case, we have $L=|\Omega|(2D+1)$ many controlled operations
to apply across the entire lattice, which incurs a T gate cost of
$4L-4=4|\Omega|(2D+1)-4$.

As described before, the CMP operations may be replaced by
$\texttt{CMP}^{\prime}$ operations described by Eq. (113). For two $n$-bit
numbers $i$ and $j$, the operation $\texttt{CMP}^{\prime}$ records the Boolean
value of $j<i$ onto an ancillary qubit, and requires $n$ many logical AND
computations. Each of these in turn can be implemented using 4 T gates. The
$\texttt{CMP}^{\prime\dagger}$ computations are simply the logical AND
uncomputations that require measurements on the ancilla and classically
conditioned post-measurement Clifford operations. In our case, we require the
comparison of two $\log_{2}{(2k)}$-bit numbers. Thus, the
$\texttt{CMP}^{\prime}$ operations contribute a total T gate count of
$4|\Omega|\log_{2}{(2k)}$ across the entire lattice.

The total T gate count for the $\phi\phi$ part of the Hamiltonian is therefore

$\displaystyle\text{Count(T)}_{\phi\phi}$ $\displaystyle=$ $\displaystyle
4|\Omega|\left(2D+\log_{2}{k}+2\right)-4$ (118)

#####  $\text{SELECT}_{\Phi^{2}}$ term

The sub-select circuit for $\Phi^{2}$ term can be constructed quite similarly.
The comparator is still used to extract the relevant phase. The only
difference is that now we must massage an ancillary register with some initial
set of operations that we summarize as $U_{initial}^{\phi^{2}}$.

Specifically, given $\ket{\psi}=\sum_{j}c_{j}\ket{j}$, denoting the field
register at some lattice point, and another ancillary register $\ket{k}$, we
require

$\displaystyle
U_{initial}^{\phi^{2}}\left(\sum_{j}c_{j}\ket{j}\right)\ket{k}\ket{scratch}\ket{0}$
$\displaystyle=$ $\displaystyle\ket{\Psi^{\prime}}\ket{k^{2}+(j+1-k)^{2}}.$
(119)

Once the value $k^{2}+(j+1-k)^{2}$ has been written onto a register, we can
repeat the procedure described earlier of using a comparator, extracting out
the relevant phase, then undoing the comparison. In sum, for each lattice site
$a$

$\displaystyle\text{SELECT}_{\Phi^{2},a}$ $\displaystyle=$
$\displaystyle\left(U_{initial}^{\phi^{2},\dagger}\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}U_{initial}^{\phi^{2}}\right)_{a}.$
(120)

As shown in the Appendix , the cost of the $U_{initial}^{\phi^{2}}$ operation
and its inverse incurs a total T gate cost of
$8|\Omega|\left(4\left(\log_{2}{2k}\right)^{2}+2\log_{2}{2k}-6\right)$, where
the leading term comes from multiplying two $\log_{2}{2k}$-bit numbers, while
other terms come from other arithmetic operations, such as addition.

The $\texttt{CMP}^{\prime}$ operation (and its inverse), replacing the CMP
operation as before, compares two $\log_{2}{2k^{2}}$-qubit numbers and costs
$\log_{2}{2k^{2}}$ many logical ANDs. Thus, the comparisons require
$|\Omega|\log_{2}{2k^{2}}$ many logical ANDs across the entire lattice,
contributing to a T gate count of $4|\Omega|\log_{2}{2k^{2}}$. In addition, we
have $|\Omega|$ many multi-controlled $Z$ gates to apply across the entire
lattice, which incurs a cost of $4|\Omega|-4$ T gates using the unary
iteration method of [70].

In all, the total T gate count for this family of terms is given by

$\displaystyle Count(T)_{\phi^{2}}$ $\displaystyle=$ $\displaystyle
8|\Omega|\left(4\left(\log_{2}{2k}\right)^{2}+2\log_{2}{2k}-6\right)+4|\Omega|\left(\log_{2}{2k^{2}}+1\right)-4$
(121) $\displaystyle=$ $\displaystyle
8|\Omega|\left(4\log_{2}^{2}{k}+11\log_{2}{k}+1\right)-4$

#####  $\text{SELECT}_{\pi^{2}}$ term

Using the relation
$\pi=\mathcal{F}^{\dagger}\phi\mathcal{F}\Rightarrow\pi^{2}=\mathcal{F}^{\dagger}\phi^{2}\mathcal{F}$,
where $\mathcal{F}$ denotes the discrete Fourier transform, we see that the
$\text{SELECT}_{\pi^{2}}$ circuit is essentially the same as the
$\text{SELECT}_{\phi^{2}}$ except that at each lattice site, we have
additional (uncontrolled) $\log_{2}{2k}$-qubit Fourier transforms and their
inverses at every lattice site. For each lattice site $a$, we have

$\displaystyle\text{SELECT}_{\pi^{2},a}$ $\displaystyle=$
$\displaystyle\left(\mathcal{F}^{\dagger}U_{initial}^{\phi^{2},\dagger}\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}U_{initial}^{\phi^{2}}\mathcal{F}\right)_{a}$
(122)

where $\mathcal{F}$ now represents the quantum Fourier transform (QFT)
circuit. Here, we use the approximate QFT (AQFT) circuit of [71], which incurs
a T gate cost of approximately

$Count(T)_{AQFT}\approx
8n\log_{2}\left(\frac{n}{\epsilon_{AQFT}}\right)+\log_{2}{\left(\frac{n}{\epsilon_{AQFT}}\right)}\log_{2}{\left(\frac{\log_{2}{\left(\frac{n}{\epsilon_{AQFT}}\right)}}{\epsilon_{AQFT}}\right)}$
(123)

for a target precision of $\epsilon_{AQFT}$. One application each of the AQFT
and its inverse at each lattice site incurs a total cost of $2|\Omega|$ times
the expression in Eq. (123) with $n=\log_{2}{2k}$. The total T gate cost for
this family of terms is then

$\displaystyle Count(T)_{\pi^{2}}$ $\displaystyle=$ $\displaystyle
Count(T)_{\phi^{2}}+Count(T)_{AQFT}$ $\displaystyle=$ $\displaystyle
8|\Omega|\left(4\log_{2}^{2}{k}+11\log_{2}{k}+1\right)-4+16|\Omega|\log_{2}{2k}\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}+\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}\log_{2}{\left(\frac{\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}}{\epsilon_{AQFT}}\right)}$

#####  $\text{SELECT}_{\phi^{4}}$ term

The sub-select circuit for the $\Phi^{4}$ term can be constructed similarly to
that of the $\Phi^{2}$ term. As before, we apply an initial transformation to
record an appropriate value to an ancillary register, to which the comparator
is applied, and subsequently the appropriate phase extracted to match the
entries of the signature matrices of the equal weight LCU Eq. (104).

Assuming $\ket{\psi}=\sum_{j}c_{j}\ket{j}$ denotes the field register at some
lattice point, and we have another ancillary register initialized to
$\ket{k}$, along with an entire $\ket{scratch}$ ancillary register to assist
with various sub-routines, as well as another ancillary register initialized
to $\ket{0}$, we require

$\displaystyle
U_{initial}^{\phi^{4}}\left(\sum_{j}c_{j}\ket{j}\right)\ket{k}\ket{scratch}\ket{0}$
$\displaystyle=$
$\displaystyle\ket{\Psi^{\prime}}\ket{k^{4}+\left(j+1-k\right)^{4}}$ (125)

Once the value $k^{4}+\left(j+1-k\right)^{4}$ has been recorded in
superposition, we can apply the usual comparator technique as described above
to extract the relevant phase and produce the desired LCU. In all, we have for
each lattice site $a$

$\displaystyle\text{SELECT}_{\Phi^{4},a}$ $\displaystyle=$
$\displaystyle\left(U_{initial}^{\phi^{4},\dagger}\texttt{CMP}^{\dagger}Z_{anc}\texttt{CMP}U_{initial}^{\phi^{4}}\right)_{a}$
(126)

As before, the $U_{initial}^{\phi^{4}}$ circuit requires the implementation of
a few basic arithmetic sub-routines, including binary multiplication of two
$O(\log{k})$-bit numbers, which contributes the dominant cost of
$O(|\Omega|(\log{k})^{2})$ T gates. As shown in the appendix, this part of the
circuit contributes a total T gate cost of
$8|\Omega|\left(20\log_{2}^{2}{(2k)}-10\log_{2}{(2k)}+1\right)$.

The cost of the $\texttt{CMP}^{\prime}$ (simplifying the CMP operation), which
now compares two $\log_{2}{2k^{4}}$-bit numbers, is
$4|\Omega|\log_{2}{2k^{4}}$. As in the case of the $\phi^{2}$ term, we have
$|\Omega|$ many multi-controlled $Z$ gates to apply across the entire lattice,
which contributes a cost of $4|\Omega|-4$ T gates using the unary iteration
method of [70].

In all, the total T gate count for this family of terms is

$\displaystyle Count(T)_{\phi^{4}}$ $\displaystyle=$ $\displaystyle
8|\Omega|\left(\log_{2}^{2}{k}-6\log_{2}{k}-7\right)-4$ (127)

The explicit construction of the $U_{initial}^{\phi^{4}}$ circuit, as well as
a detailed counting of the $T$ gates, is provided in the Appendix.

##### Total cost for Hamiltonian simulation

The constructions above create a block encoding of the Hamiltonian via
$\left(\bra{G}\otimes\mathbb{I}\right)U\left(\ket{G}\otimes\mathbb{I}\right)=H$
where $\ket{G}=U_{PREP}\ket{\overline{0}}$ and $U=U_{SELECT}$. Since
$U_{SELECT}^{2}=\mathbb{I}$, by Corollary 9 of [5], we have that the walk
operator

$W=\left(\left(2U_{PREP}\ket{0}\bra{0}U_{PREP}^{\dagger}-I\right)\otimes
I\right)U_{SELECT}$ (128)

also provides a block encoding of the Hamiltonian, i.e.
$\left(\bra{G}\otimes\mathbb{I}\right)W\left(\ket{G}\otimes\mathbb{I}\right)=H$
but in an $SU(2)$ invariant subspace containing $\ket{G}$, i.e. in qubitized
form. It is this qubitized block encoding that we shall employ to perform
phase estimation in order to subsequently estimate scattering matrix elements.

The cost of this block encoding is then given by the sum of twice the cost of
the PREP oracle, the cost of the SELECT oracle, and the cost of the reflection
operator $2\ket{0}\bra{0}-\mathbb{I}$. The cost of PREP circuit is given by
Eq. (109), while the cost of the SELECT circuit is given by the sum of Eqs.
(118), (121), (LABEL:eqn:tgate_cost_pisqr) and (127). The cost of the
reflection operator is at most logarithmic in all parameters, and is
exponentially sub-dominant to the other costs, so we neglect it in our
analysis. In all, we have

$\displaystyle Cost(T)_{W}$ $\displaystyle=$ $\displaystyle
72|\Omega|\log_{2}^{2}{k}+168|\Omega|\log_{2}{k}+8|\Omega|D-32|\Omega|+73.608\log_{2}{\left(2/\epsilon_{synth}\right)}$
(129)
$\displaystyle+16|\Omega|\log_{2}{2k}\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}+\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}\log_{2}{\left(\frac{\log_{2}{\left(\frac{\log_{2}{2k}}{\epsilon_{AQFT}}\right)}}{\epsilon_{AQFT}}\right)}-119.848$

In order to employ this block-encoding for Hamiltonian simulation using QSVT,
we need

$\Theta\left(|\alpha|_{1}t+\frac{\log{1/\epsilon}}{\log{\left(e+\frac{\log{(1/\epsilon)}}{|\alpha|_{1}t}\right)}}\right)$
(130)

many queries to the walk operator $W$, where the coefficient 1-norm is given
by Eq. (108).

##### Ancilla count

To implement the operations described above, we need several ancillary qubits.
First, we need 2 ancillae for the 4 families of terms in the Hamiltonian,
$\log_{2}{|\Omega|}$ ancillae for the $|\Omega|$ sites on the lattice, and
$\log_{2}{(2D)}$ ancillae for the neighbors of a given lattice site. These are
common across all the $\text{SELECT}_{i,a}$ operations, even though the last
register is only used for the $\text{SELECT}_{\phi\phi,a}$ operation.

Second, the comparison operations can only be carried out in succession, and
may thus be recycled for use at each lattice site (or lattice site-neighbor
pair in the case of the $\text{SELECT}_{\phi\phi,a}$ operation). The largest
comparisons we require are for the $\text{SELECT}_{\phi^{4},a}$ operation.
Each such comparison requires $\log_{2}{(2k^{4})}$ qubits to compare against,
another $\log_{2}{(2k^{4})}$ qubits to hold temporary carry values, and 1 more
qubit for the $Z_{anc}$ or $(XZX)_{anc}$ operation. The controlled-$Z$
operation further requires at most $\log_{2}{(|\Omega|D)}$ many ancillary
qubits.

Similarly, the ancillae involved in the various $U_{initial}$ operations can
be recycled for use across each site and across each family of terms. The
dominant cost is for the $\phi^{4}$ term, and provided in the appendix. In
addition to the ancilla qubits, we also have $|\Omega|\log_{2}{(2k)}$ many
qubits to hold the values of the scalar field across the entire lattice.
Adding all these counts together, we find

$\text{Count(Qubits)}=|\Omega|\log_{2}{(2k)}+18\log^{2}_{2}{k}+60\log_{2}{k}+\log_{2}{(|\Omega|D)}+29\in
O\left(|\Omega|\log{k}+\log^{2}{k}\right)$ (131)

####  IV.5.2 Algorithm II : Trotterization with Z operators

In this section we describe a decomposition of the operators
$\Phi,\Phi^{2},\Phi^{4}$ as a function of Z operators, thus making it amenable
to Trotterization. In the next section we will use these LCU decompositions to
estimate the simulation cost using QSVT or QSP. The sum of the $\ell_{1}$ norm
of the coefficients in these decompositions is referred to as the $\ell_{1}$
norm of the decomposition or operator. This is also an upper bound on the
spectral norm of the operators and is useful to bound the complexity of LCU
based simulation algorithms.

##### Decomposition of $\Phi$, $\Phi^{2}$ and $\Phi^{4}$ :

We decompose $\frac{\Phi}{\Delta\Phi}$ as follows, as done in Eq. 95.

$\displaystyle\frac{\Phi}{\Delta\Phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j=0}^{\log_{2}k}2^{j}\text{Z}_{j}-\frac{1}{2}\mathbb{I}$

Thus,

$\displaystyle\left(\frac{\Phi}{\Delta\Phi}\right)^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\Phi}{\Delta\Phi}+\frac{\mathbb{I}}{2}\right)^{2}-\left(\frac{\Phi}{\Delta\Phi}+\frac{\mathbb{I}}{2}\right)+\frac{\mathbb{I}}{4}$
(132) $\displaystyle=$
$\displaystyle\frac{1}{6}\left(2^{2\log_{2}k+1}+1\right)\mathbb{I}+\sum_{j=0}^{\log_{2}k-1}\sum_{k>j}^{\log_{2}k}2^{j+k-1}\text{Z}_{j}\text{Z}_{k}+\frac{1}{2}\sum_{j=0}^{\log_{2}k}2^{j}\text{Z}_{j}$

and so $\left(\frac{\Phi}{\Delta\Phi}\right)^{4}$ can be obtained by squaring
the above equation and hence can be expressed as sum of Z operators as well.
The $\ell_{1}$ norm of the decomposition of
$\left(\frac{\Phi}{\Delta\Phi}\right)$,
$\left(\frac{\Phi}{\Delta\Phi}\right)^{2}$ and
$\left(\frac{\Phi}{\Delta\Phi}\right)^{4}$ is $k$, $k^{2}$ and $k^{4}$,
respectively. The detail of these calculations have been provided in Appendix
C.2.

Plugging in the above LCU decomposition of the operators into Eq. 11, we
obtain the following decomposition of $H_{amp}$. We keep in mind that we have
a lattice with $|\Omega|$ vertices, where each vertex has $D$ neighbors. Let
us denote the set of edges by $E_{D}$.

$\displaystyle H_{amp}$ $\displaystyle=$ $\displaystyle H_{\pi}+H_{\phi}$
(133) $\displaystyle H_{\pi}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}H_{\pi\mathbf{x}}:=\sum_{\mathbf{x}\in\Omega}\frac{1}{2}\Pi^{2}(\mathbf{x})=\sum_{\mathbf{x}\in\Omega}\mathcal{F}\left(\gamma_{0}\mathbb{I}+\sum_{j}\gamma_{j}Z_{j}+\sum_{j,k}\gamma_{jk}Z_{j}Z_{k}\right)\mathcal{F}^{\dagger}$
(134) $\displaystyle H_{\phi^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}\left(\mathbb{I}+\sum_{j}\alpha_{j}Z_{j}+\sum_{j,k}\alpha_{jk}Z_{j}Z_{k}+\sum_{j,k,l}\alpha_{jkl}Z_{j}Z_{k}Z_{l}+\sum_{j,k,l,m}\alpha_{jklm}Z_{j}Z_{k}Z_{l}Z_{m}\right)$
(135) $\displaystyle H_{\phi^{\prime\prime}}$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{x},\mathbf{x}^{\prime}\\\
(\mathbf{x},\mathbf{x}^{\prime})\in
E_{D}\end{subarray}}\sum_{j,j^{\prime}}\beta_{jj^{\prime}}\left(Z_{j}\right)_{\mathbf{x}}\left(Z_{j^{\prime}}\right)_{\mathbf{x}^{\prime}}$
(136) $\displaystyle H_{\phi}$ $\displaystyle=$ $\displaystyle
H_{\phi^{\prime}}+H_{\phi^{\prime\prime}}$ (137)

In the above set of equations the coefficients
$\gamma_{0},\gamma_{j},\gamma_{jk}$, $\alpha_{j},\alpha_{jk},\alpha_{jkl}$,
$\alpha_{jklm},\beta_{jj^{\prime}}$ are determined by the LCU decomposition of
the operators as well as Eq. 11. We ignore the identity term, since if
$H_{amp}=H_{amp^{\prime}}+\beta_{0}\mathbb{I}$, where $\beta_{0}$ is a
constant, then $e^{-iH_{amp}^{\prime}\tau}\propto e^{-iH_{amp}\tau}$ upto some
global phase. So the gate complexity to implement the exponentials is the
same.

We have two non-commuting components, $H_{\phi}$ and $H_{\pi}$. We first
analyze the gate complexity to implement $e^{-iH_{\phi}\tau}$ and
$e^{-iH_{\pi}\tau}$, where $\tau$ is one Trotter step. As mentioned,
Trotterization becomes quite straight-forward with this decomposition. We
require 1 $R_{z}$ gate and few CNOT gates (to compute the parity) for each
exponentiated Z operator. Thus we need at most

$\displaystyle|\Omega|\left(\binom{\log_{2}k+1}{1}+\binom{\log_{2}k+1}{2}+\binom{\log_{2}k+1}{3}+\binom{\log_{2}k+1}{4}\right)+E_{D}(\log_{2}k+1)^{2}$
(138)

number of $R_{z}$ gates and at most

$\displaystyle
2|\Omega|\left(\binom{\log_{2}k+1}{2}+2\binom{\log_{2}k+1}{3}+3\binom{\log_{2}k+1}{4}\right)+2E_{D}(\log_{2}k+1)^{2}$
(139)

number of CNOT gates. We observe that in order to implement the required
exponential of the operators we require a number of CNOT that will realize the
necessary paritites i.e. the XOR of all combinations of 2, 3 and 4 qubits.
Such a circuit is usually referred to as a ’parity network’. The $R_{z}$ gates
are placed at that point of the circuit where the corresponding parity is
realized. The number of CNOT gates can be optimized using algorithms like [81,
82]. Additionally we require $2|\Omega|$ number of $\log_{2}2k$-qubit QFT. If
we use the approximate QFT (AQFT) circuit of [71] then we require

$\displaystyle
2|\Omega|\left(8(\log_{2}2k)\log_{2}\left(\frac{\log_{2}2k}{\epsilon_{QFT}}\right)+\log_{2}\left(\frac{\log_{2}2k}{\epsilon_{QFT}}\right)\log_{2}\left(\frac{\log_{2}\left(\frac{\log_{2}2k}{\epsilon_{QFT}}\right)}{\epsilon_{QFT}}\right)\right)$
(140)

number of additional T gates.

Now we bound the error using the $2^{nd}$ order Trotter formula [68, 67], in a
similar manner as discussed in Section IV.3. We need to calculate the second
level nested commutators. Here, for convenience we fragment $H_{\phi}$ as
follows.

$\displaystyle H_{\phi}$ $\displaystyle:=$ $\displaystyle
H_{1\phi}+H_{2\phi}+H_{3\phi}$ (141) $\displaystyle H_{1\phi}$
$\displaystyle:=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}H_{1\phi\mathbf{x}}=\sum_{\mathbf{x}\in\Omega}\frac{1}{2}(M^{2}+d+1)\Phi^{2}(\mathbf{x})$
(142) $\displaystyle H_{2\phi}$ $\displaystyle:=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}H_{2\phi\mathbf{x}}=\sum_{\mathbf{x}\in\Omega}\frac{\Lambda}{24}\Phi^{4}(\mathbf{x})$
(143) $\displaystyle H_{3\phi}$ $\displaystyle:=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}H_{3\phi\mathbf{x}}=\sum_{\mathbf{x}\in\Omega}\left(-2\sum_{\begin{subarray}{c}\mathbf{x}^{\prime}:\\\
(\mathbf{x},\mathbf{x}^{\prime})\in
E_{D}\end{subarray}}\Phi(\mathbf{x})\Phi(\mathbf{x}^{\prime})\right):=\sum_{\mathbf{x}\in\Omega}\sum_{\begin{subarray}{c}\mathbf{x}^{\prime}:\\\
(\mathbf{x},\mathbf{x}^{\prime})\in
E_{D}\end{subarray}}H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}$ (144)

We observe the following commutation relations between $H_{\pi\mathbf{x}}$
(Eq. 134) and the Hamiltonian fragments in the above equations.

$\displaystyle\left[H_{1\phi\mathbf{x}},H_{2\phi\mathbf{x}^{\prime}}\right],\left[H_{1\phi\mathbf{x}},H_{3\phi\mathbf{x}^{\prime}}\right],\left[H_{3\phi\mathbf{x}},H_{2\phi\mathbf{x}^{\prime}}\right]$
$\displaystyle=$ $\displaystyle 0\qquad\text{for
all}\quad\mathbf{x},\mathbf{x}^{\prime},$
$\displaystyle\left[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}^{\prime}}\right],\left[H_{\pi\mathbf{x}},H_{2\phi\mathbf{x}^{\prime}}\right]$
$\displaystyle\neq$ $\displaystyle 0\qquad\text{if and only if
}\quad\mathbf{x}=\mathbf{x}^{\prime},$
$\displaystyle\text{and}\qquad\left[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}^{\prime}}\right]$
$\displaystyle\neq$ $\displaystyle 0\qquad\text{if and only if
}\quad(\mathbf{x},\mathbf{x}^{\prime})\in E_{D}.$ (145)

Now assuming
$H_{1\phi}^{\prime}=\sum_{\mathbf{x}\in\Omega}H_{1\phi\mathbf{x}}^{\prime}=\sum_{\mathbf{x}\in\Omega}H_{1\phi\mathbf{x}}+H_{2\phi\mathbf{x}}=H_{1\phi}+H_{2\phi}$,
we have

$\displaystyle\widetilde{\alpha}_{comm}$ $\displaystyle\leq$
$\displaystyle\left(\|[[H_{\pi},H_{1\phi}^{\prime}],H_{\pi}]+[[H_{\pi},H_{1\phi}^{\prime}],H_{1\phi}^{\prime}]\|\right)$
(146)
$\displaystyle+\left(\|[[H_{\pi},H_{1\phi}^{\prime}],H_{3\phi}]+[[H_{\pi},H_{3\phi}],H_{\pi}]+[[H_{\pi},H_{3\phi}],H_{1\phi}^{\prime}]+[[H_{\pi},H_{3\phi}],H_{3\phi}]\|\right)$
$\displaystyle:=$ $\displaystyle S_{1}+S_{2}.$

For the first sum, we use Eq. (145) and observe that Eq. (63) can be applied
along with triangle inequality in order to have the following bound.

$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{x}\in\Omega}\|[[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}],H_{\pi\mathbf{x}}]\|+\|[[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}],H_{1\phi\mathbf{x}}^{\prime}]\|$
$\displaystyle\leq$
$\displaystyle\sum_{\mathbf{x}\in\Omega}\|[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}]\|\left(\|H_{1\phi\mathbf{x}}^{\prime}\|+\|H_{\pi\mathbf{x}}\|\right)\leq|\Omega|\|[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}]\|\left(\|H_{1\phi\mathbf{x}}^{\prime}\|+\|H_{\pi\mathbf{x}}\|\right)$
$\displaystyle\leq$
$\displaystyle\sum_{\mathbf{x}\in\Omega}\left(\frac{M^{2}+d+1}{4}\|[\Pi^{2}(\mathbf{x}),\Phi^{2}(\mathbf{x})]\|+\frac{\Lambda}{48}\|[\Pi^{2}(\mathbf{x}),\Phi^{4}(\mathbf{x})]\|\right)\left(\frac{1}{2}(M^{2}+d+1)\|\Phi^{2}\|+\frac{\Lambda}{24}\|\Phi^{4}\|+\frac{1}{2}\|\Pi^{2}\|\right)$
$\displaystyle\leq$
$\displaystyle|\Omega|\left(\frac{M^{2}+d+1}{2}\|\Pi^{2}(\mathbf{x})\|\|\Phi^{2}(\mathbf{x})\|+\frac{\Lambda}{24}\|\Pi^{2}(\mathbf{x})\|\|\Phi^{4}(\mathbf{x})\|\right)\left(\frac{1}{2}(M^{2}+d+1)\|\Phi^{2}\|+\frac{\Lambda}{24}\|\Phi^{4}\|+\frac{1}{2}\|\Pi^{2}\|\right)$

For the second sum $S_{2}$ we again use Eq. (145) and triangle inequality and
obtain the following bound.

$\displaystyle S_{2}$ $\displaystyle\leq$
$\displaystyle\sum_{\mathbf{x}\in\Omega}\sum_{\begin{subarray}{c}\mathbf{x}^{\prime}:\\\
(\mathbf{x},\mathbf{x}^{\prime})\in
E_{D}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{x}^{\prime\prime}:(\mathbf{x},\mathbf{x}^{\prime\prime})\in
E_{D}\\\ \text{or }(\mathbf{x}^{\prime},\mathbf{x}^{\prime\prime})\in
E_{D}\end{subarray}}\left(\|[[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}],H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}]\|+\|[[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}],H_{\pi\mathbf{x}^{\prime\prime}}]\|+\|[[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}],H_{1\phi\mathbf{x}^{\prime\prime}}^{\prime}]\|\right.$
(148)
$\displaystyle+\left.\|[[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}],H_{3\phi\mathbf{x}\mathbf{x}^{\prime\prime}}]\|+\|[[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}],H_{3\phi\mathbf{x}^{\prime}\mathbf{x}^{\prime\prime}}]\|\right)$
$\displaystyle\leq$
$\displaystyle|\Omega|d^{2}\left(\|[H_{\pi\mathbf{x}},H_{1\phi\mathbf{x}}^{\prime}]\|\|H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}\|+\|[H_{\pi\mathbf{x}},H_{3\phi\mathbf{x}\mathbf{x}^{\prime}}]\|\left(\|H_{\pi\mathbf{x}^{\prime\prime}}\|+\|H_{1\phi\mathbf{x}^{\prime\prime}}^{\prime}\|+\|H_{3\phi\mathbf{x}\mathbf{x}^{\prime\prime}}\|+\|H_{3\phi\mathbf{x}^{\prime}\mathbf{x}^{\prime\prime}}\|\right)\right)$

We already mentioned that the $\ell_{1}$ norm for the decomposition of