# Simulating collider physics on quantum computers using effective field
theories

Christian W. Bauer<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Marat Freytsis
<EMAIL_ADDRESS>NHETC, Department of Physics and Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  Physics Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Benjamin Nachman
<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA

###### Abstract

Simulating the full dynamics of a quantum field theory over a wide range of
energies requires exceptionally large quantum computing resources. Yet for
many observables in particle physics, perturbative techniques are sufficient
to accurately model all but a constrained range of energies within the
validity of the theory. We demonstrate that effective field theories (EFTs)
provide an efficient mechanism to separate the high energy dynamics that is
easily calculated by traditional perturbation theory from the dynamics at low
energy and show how quantum algorithms can be used to simulate the dynamics of
the low energy EFT from first principles. As an explicit example we calculate
the expectation values of vacuum-to-vacuum and vacuum-to-one-particle
transitions in the presence of a time-ordered product of two Wilson lines in
scalar field theory, an object closely related to those arising in EFTs of the
Standard Model of particle physics. Calculations are performed using
simulations of a quantum computer as well as measurements using the IBMQ
Manhattan machine.

It is well known that quantum computers can in principle simulate the time
evolution of quantum field theories Jordan:2011ne . The main technique
involves disretizing the spatial degrees of freedom by introducing a lattice
Kogut:1974ag ; Kogut:1979wt ; Bronzan:1984xb , and digitizing the field values
at a given lattice point Jordan:2011ci ; Macridin:2018gdw ; Macridin:2018oli ;
Klco:2018zqz ; Hackett:2018cel ; Yeter-Aydeniz:2018mix ; Kreshchuk:2020dla ;
Kreshchuk:2020kcz ; Haase:2020kaj . This turns the uncountably infinite
dimensional Hilbert space of standard quantum field theories into a finite
dimensional Hilbert space of dimension

$\displaystyle n_{H}=n_{\phi}^{\pqty{N^{d}}}\,,$ (1)

where $n_{\phi}$ denotes the dimensionality of the Hilbert space for a given
lattice point, $N$ is the total number of lattice points in each spatial
direction, and $d$ represents the number of spatial dimensions. The physical
volume of the lattice is determined by the distance between adjacent lattice
points $\variation{x}$ in each direction and is given by

$\displaystyle V=(N\variation{x})^{d}\equiv L^{d}\,.$ (2)

The total number of qubits required for such a simulation is given by

$\displaystyle
n_{Q}=\mathcal{O}(\log_{2}n_{H})=\mathcal{O}(N^{d}\log_{2}n_{\phi})\,.$ (3)

The discretization and finite volume of space introduce upper and lower
cutoffs to the energies $E$ over which the resulting lattice field theory is a
good approximation for the continuum. In particular, one finds

$\displaystyle\frac{1}{N\variation{x}}\lesssim
E\lesssim\frac{1}{\variation{x}}\,,$ (4)

which implies that the range of energies that can be described is directly
proportional to the number of lattice sites per dimension. In principle, to
have access to the full dynamics of the Large Hadron Collinder (LHC), one
would need to describe the energy range between
$\mathcal{O}($10\text{\,}\mathrm{MeV}$)$ (the smallest resolvable transverse
momentum between hadrons) and $7\text{\,}\mathrm{TeV}$ (the beam energy of the
LHC). To fully capture this energy range would require a lattice with
$\mathcal{O}(${10}^{6}$)$ lattice points in each dimension, and more than
${10}^{18}$ qubits to reproduce the resulting physical system. Even if one
does not require the full energy range up to the LHC center-of-mass energy,
the number of qubits required for a full simulation will clearly remain beyond
the realm of feasibility for a long time to come.

For many observables of interest at the LHC, physics at short distances is
reliably computed in fixed order perturbation theory, and high precision can
be reached with existing techniques. Physics at lower energies introduces new
significant challenges. Asymptotic freedom Gross:1973id ; Politzer:1973fx
implies that the strong coupling constant becomes large at low energies,
increasing the production of additional particles. An energetic particle can
thus easily radiate additional soft and collinear particles, resulting in a
collimated jet of particles. This means that fixed order perturbative
calculations are of little use to describe lower energy effects such as the
precise makeup of jets, and other techniques such as resummation
Collins:1984kg ; Catani:1992ua ; Bonciani:2003nt  and parton showers
Marchesini:1983bm ; Bengtsson:1986et  have to be used to make predictions. At
even lower energies, the fully nonperturbative dynamics of hadronization
dominate. While classical lattice methods provide spectral information about
the strongly-coupled limit of the theory, the dynamics of hadronization is
currently only understood through phenomenological models and form factors
extracted from data. While these techniques have been quite successful and it
has even been shown that quantum algorithms can be used to include quantum
interference effects in some models of parton showers Bauer:2019qxa , it would
be a significant breakthrough were it possible to simulate these low energy
dynamics from first principles.

In this paper we show that this is indeed possible by using effective field
theories (EFTs) which have been designed to reproduce the desired long
distance physics. (For work on simulating EFTs not related to collider
physics, see Dumitrescu:2018njn ; Lu:2018pjk ; Cervia:2019res ;
Roggero:2019myu ; Cervia:2020fkk ; Stetina:2020abi ). The dynamics of these
EFTs can then be simulated directly and from first principles on a quantum
computer. Since the energy range that needs to be simulated is much smaller
than that of the full theory, the resource requirements are orders of
magnitudes smaller than those required to simulate the full theory.
Additionally, the EFT simplifies the description of the initial and final
state and significantly reduces the resource requirements of their
implementation. For example, consider an observable measured on two jets, each
with an energy of $\mathcal{O}($1\text{\,}\mathrm{TeV}$)$ and an invariant
mass of $\mathcal{O}($50\text{\,}\mathrm{GeV}$)$. Even restricting ourselves
to observables insensitive to hadronization, a generic observable of the
momenta of final state particles in the full theory requires a simulation of
the range $$1\text{\,}\mathrm{GeV}$\lesssim
E\lesssim$2\text{\,}\mathrm{TeV}$$. As we discuss below, the EFT only needs to
capture $$1\text{\,}\mathrm{GeV}$\lesssim E\lesssim$50\text{\,}\mathrm{GeV}$$.
The number of qubits required to simulate the EFT is smaller by a factor of
$(2000/50)^{3}\approx${10}^{5}$$, and this factor increases rapidly as the
range of the full theory is increased.

The EFT relevant for hadronic jet physics is Soft-Collinear Effective Theory
(SCET) Bauer:2000ew ; Bauer:2000yr ; Bauer:2001ct ; Bauer:2001yt . Using this
EFT, a typical cross section describing a physical observable at the LHC can
be factorized into separate pieces Bauer:2002nz ; Bauer:2002ie ; Bauer:2003di
; Manohar:2003vb

$\displaystyle\sigma=H\otimes J_{1}\otimes\dotsb\otimes J_{n}\otimes S\,.$ (5)

Here $H$ denotes a coefficient describing the short distance physics which can
be computed reliably in standard perturbation theory, while $J_{i}$ and $S$
denote squares of matrix elements of operators in SCET. The jet function
$J_{i}$ denotes physics arising from collinear degrees of freedom that all
have small momentum relative to each other (while moving collectively with a
large energy), with the different jet functions completely decoupled from one
another. The soft function $S$ denotes physics arising from soft degrees of
freedom, which all have small absolute momentum (in the frame of the
collider). These matrix elements describe the transition of a simple initial
state produced in the short distance interaction $H$ into a final state, with
the dynamics described by the EFT Hamiltonian. The final states that arise can
contain a large number of particles. State of the art techniques use operator
renormalization to compute the overall scale dependence of the jet and soft
functions, and then compute the relevant matrix elements perturbatively, such
that the effects arising from the high multiplicity final states or
hadronization can not be properly included. Classical lattice techniques are
also not suitable for the computation of matrix elements of SCET operators,
since the long-distance dynamics is governed by massless modes which are
inherently Minkowskian in nature. Having a simulation of the dynamics of SCET
(or a different EFT for other problems) will allow for a full non-perturbative
calculation of any matrix elements.

The EFT reproduces the full theory result up to power corrections, whose size
depends on the kinematics of the process studied. For many observables of
interest, these power corrections are considerably smaller than the
perturbative corrections present in either the full theory or EFT matrix
elements.

Since individual jet functions do not interact with one another, their
dynamics can be simulated in the reference frame where all particles have
small absolute momentum Becher:2006qw ; Goerke:2017ioi . This requires
simulating the dynamics of the original field theory (with any degrees of
freedom that only contribute to short distance effects removed), but with a
much smaller energy range required to be simulated. To achieve this on a
Quantum Computer, one needs a setup very similar to that for the full theory,
but reliable calculations can be achieved with much coarser lattices than
those one would need for a simulation of the full theory.

In the soft function, on the other hand, the energetic particles no longer
contribute to the dynamics of the theory. Instead, their effect is captured by
a so-called Wilson line, which describes the interaction of a charge moving
along a fixed world-line with the bath of soft degrees of freedom. The physics
underlying this observation is that soft particles cannot change the direction
of energetic particles in a meaningful way. This implies that energetic
particles can be included in the EFT as a static object (Wilson line)
described by the relevant quantum numbers (charge, color, etc.) moving along a
fixed world line along a light-like direction. Matrix elements in the soft
theory therefore need to compute the dynamics of a Hamiltonian describing the
soft bath of particles in the presence of such Wilson lines. For gauge
theories, such as those describing the three fundamental forces contained in
the Standard Model, Wilson lines are given by path ordered exponentials of the
relevant gauge fields Bauer:2001yt

$\displaystyle
Y_{n}=\mathrm{P}\exp\bqty{ig\int_{0}^{\infty}\\!\\!\differential{s}\,n\\!\cdot\\!A(x^{\mu}=n^{\mu}s)}\,.$
(6)

Here $A^{\mu}$ is the soft gauge field and $n^{\mu}$ describes the direction
of the light-like direction $n^{\mu}=(1,\mathbf{n})$ with $\mathbf{n}^{2}=1$
such that $n_{\mu}n^{\mu}=0$. Thus, the gauge field is evaluated on a path
going from the origin to spatial infinity along the world-line characterized
by the direction $n^{\mu}$. The dynamics of the soft gauge field is described
by its full theory Hamiltonian.

The soft function for a process containing two energetic particles of zero
total charge moving back-to-back requires two Wilson lines, $Y_{n}$ for the
particle moving in the $n$ direction and $Y_{\bar{n}}^{\dagger}$ for the
antiparticle moving in the ${\bar{n}}^{\mu}=(1,-\mathbf{n})$ direction. The
matrix elements required in the soft sector are given by

$\displaystyle\matrixelement{X}{\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]}{{\Omega}}\,,$
(7)

where $\mathrm{T}$ denotes the time-ordering operator, $\ket{{\Omega}}$
denotes the ground state of a Hilbert space containing only the gauge degree
of freedom, and $\ket{X}$ is the final state. If reproducing multi-particle
weakly coupled final states, $\ket{X}$ will contain a fixed number of gauge
momentum modes with given momenta $k^{\mu}_{i}$. To compute the soft function
for a given observable one needs to sum over all possible final states that
can contribute to this observable. Traditional perturbative calculations
Becher:2005pd ; Hoang:2008fs ; Jouttenus:2011wh ; Kelley:2011ng ; Monni:2011gb
; Boughezal:2015eha ; Moult:2018jzp  can only compute these matrix elements
for final states $\ket{X}$ containing a small number of particles and at low
order in perturbation theory, since the complexity of the calculation
increases factorially with the power of the coupling constant $g$. They
therefore do not compute the full matrix element for a given observable, but
only a perturbative approximation. For large coupling constants, which is the
relevant case for Quantum Chromodynamics (QCD) at low energies, such
perturbative calculations can give rise to large uncertainties.

Simulating the full dynamics of the field theory would provide the full non-
perturbative result for the soft matrix element. To evaluate the matrix
element on a quantum computer one first needs to define circuits that perform
time evolution of the system, as well as circuits that can create and measure
the ground and exited states of the theory. In addition, one needs to create
circuits that can correctly interleave the implementation of nonlocal Wilson
line operators with the evolution of the system to reproduce their time-
ordered product.

In the following we discuss how to compute matrix elements of Wilson line
operators $Y_{n}$ analogous to Eq. 7, but for a massless scalar, rather than
gauge, theory. The EFT is insensitive to the precise origin of the Wilson
lines, but a particualy straightforward realization would result from a pair
of highly-energetic fermions coupled to massless scalars through a Yukawa
coupling. When constructing the explicit circuit we also limit ourselves to
(1+1) dimensions, mainly to restrict the quantum resources required such that
it can be implemented on currently existing hardware. This allows us to omit
some technical complications that arise when dealing with gauge theories
(gauge transformations, the existence of unphysical polarizations, _etc_.),
while capturing all the resulting simplification of working within an EFT.

To be precise, we consider a massless field theory in (1+1) dimensions with
Hamiltonian and Wilson lines defined by

$\displaystyle H$
$\displaystyle=\int\differential{x}\frac{1}{2}\pqty{\dot{\phi}^{2}-\phi\,\partial^{2}\phi}\,,\quad$
$\displaystyle Y_{n}$
$\displaystyle=\mathrm{P}\exp\bqty{ig\int_{0}^{\infty}\\!\\!\differential{s}\,\phi(x^{\mu}=n^{\mu}s)}\,.$
(8)

In the Supplemental Material, we discuss how working in (1+1) dimensions gives
rise to several effects not present in higher dimensions, and how these
generalize to higher dimensions.

We discretize the position $x$ into an odd number of lattice points, labeling
the positions by $x_{0},\dotsc,x_{N-1}$. To eliminate the zero-momentum mode
of the theory, we impose twisted boundary conditions Lin:2001zz ;
Sachrajda:2004mi ; Bedaque:2004kc ; Briceno:2013hya . The result is a theory
defined at discrete positions $x$ and momenta $p$ given by
$x_{i}=x_{\mathrm{min}}+i\variation{x}$ and
$p_{i}=p_{\mathrm{min}}+i\variation{p}$ with
$x_{\mathrm{min}}=-(N-1)\variation{x}/2$,
$p_{\mathrm{min}}=-\pi/\variation{x}$ and $\variation{p}=2\pi/N\variation{x}$,
Writing $\phi_{i}\equiv\phi(x_{i})$, the twisted boundary conditions
correspond to the condition $\phi_{i+N}=-\phi_{i}$. The Hamiltonian becomes
Jordan:2011ci

$\displaystyle
H=\frac{\variation{x}}{2}\sum_{i=0}^{N-1}\bqty{\dot{\phi}_{i}^{2}-\phi_{i}\,[\nabla^{2}\phi]_{j}}\,,$
(9)

where the lattice operator $\nabla^{2}$ is defined through its action on a
field as
$[\nabla^{2}\phi]_{i}=(2\phi_{i}-\phi_{i-1}-\phi_{i+1})/\variation{x}^{2}$.
Due to the twisted boundary conditions
$[\nabla^{2}\phi]_{0}=(2\phi_{0}+\phi_{N-1}-\phi_{1})/\variation{x}^{2}$ and
$[\nabla^{2}\phi]_{N-1}=(2\phi_{N-1}-\phi_{N-2}+\phi_{0})/\variation{x}^{2}$.
The Wilson line operators can be written as

$\displaystyle Y_{n}$
$\displaystyle=\mathrm{P}\exp\bqty{ig\variation{x}\sum_{i=n_{0}}^{2n_{0}}\phi_{x_{i}}(t=x_{i}-n_{0})}\,,\quad$
$\displaystyle Y_{{\bar{n}}}^{\dagger}$
$\displaystyle=\mathrm{P}\exp\bqty{-ig\variation{x}\sum_{i=0}^{n_{0}}\phi_{x_{i}}(t=n_{0}-x_{i})}\,,$
(10)

where $n_{0}=(N-1)/2$ denotes the point at the center of the lattice.

We represent the field theory through the field values at each lattice
position, and in order to describe the theory on a digital quantum computer
one needs to digitize the continuous field value at each lattice point
Klco:2018zqz . Choosing $n_{Q}$ qubits per lattice site allows for
$n_{\phi}\equiv 2^{n_{Q}}$ different field values. For each lattice point, the
possible field values are chosen to be by
$\phi_{i}^{(k)}=-\phi_{\mathrm{max}}+k\,\delta\phi$, with
$\delta\phi=2\phi_{\mathrm{max}}/(n_{\phi}-1)$. The value of
$\phi_{\mathrm{max}}$ has to be chosen to optimize the digitized description,
which for free fields is accomplished by

$\displaystyle\phi_{\mathrm{max}}=\frac{1}{\sqrt{\delta
x\,\bar{\omega}}}\sqrt{\frac{\pi}{2}\frac{(n_{\phi}-1)^{2}}{n_{\phi}}}\,,$
(11)

where

$\displaystyle\bar{\omega}=\frac{1}{N}\sum_{i}\omega_{i}\,,\quad\omega_{i}=\frac{2}{\variation{x}}\absolutevalue{\sin\frac{p_{i}\variation{x}}{2}}\,.$
(12)

For $\bar{\omega}=1$, as is the case for a single lattice site with
$\omega=1$, corresponding to a single harmonic oscillator, Eq. 11 reproduces
the empirical numerical values obtained in Klco:2018zqz .

To implement the Wilson line operator we first rewrite the time-ordered
product of the two Wilson lines as

$\displaystyle\mathrm{T}[Y_{n}\,Y_{\bar{n}}^{\dagger}]$
$\displaystyle=e^{-iH\,n_{0}\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{x_{2n_{0}}}-\phi_{x_{0}}}}$
(13) $\displaystyle\qquad\times
e^{iH\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{x_{2n_{0}-1}}-\phi_{x_{1}}}}$
$\displaystyle\qquad\times\dotsm\times
e^{iH\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{x_{n_{0}}}-\phi_{x_{n_{0}}}}}\,,$

where we have used the time translation operator to make the time dependence
on the field operators explicit. Thus, the Wilson line operator consists of a
sequence of time-evolution operators for a time interval corresponding to the
lattice spacing and exponentials of the field operator. The last time
evolution evolves the state back from the largest time to which the Wilson
lines can be sensibly evolved, namely $t_{\mathrm{max}}=n_{0}\variation{x}$,
to $t=0$ at which all states are defined.

Ultimately, to make contact with the continuum field theory any such
simulation will have to be performed on a series of increasing lattices, and
the result extrapolated to the $N\to\infty$, $\variation{x}\to 0$ limit. Any
parameters of the theory present in the continuum must be suitably matched for
this procedure to yield meaningful results. For local terms in the
Hamiltonian, this procedure is discussed in detail in Jordan:2011ci . Dealing
with a massless theory simplifies this procedure since only local interactions
(of which in the present case there are none) need to be matched. However, the
EFT will also require the matching of Wilson line operators, which is
complicated by their non-local nature and sensitivity to total lattice size,
as discussed in the Supplemental Materials. In this letter, we work at fixed
lattice size and we leave the detailed investigation of these issues to future
work.

The implementation of the exponential of the field operator, as well as the
time evolution operator, follows the discussion in Klco:2018zqz  and uses the
fact that the digitized field $\phi_{i}^{(k)}$ can be written in terms of sums
of $\sigma_{z}$ operators. This implies that the exponential of products of
fields $\phi_{i}$ can be implemented through combinations of CNOT gates and
$R_{Z}$ rotations Macridin:2018gdw ; Macridin:2018oli ; Klco:2018zqz . The
exponential of the conjugate operator $\dot{\phi}^{2}$ can be implemented by
taking a quantum Fourier transformation of the exponential of the operator
$\phi^{2}$. These can then be combined via the Suzuki–Trotter formula
10.2307/2033649 ; Suzuki1976 ; 1976PThPh..56.1454S . For details, see the
Supplemental Materials. The initial ground state of the scalar field theory is
a multi-variate Gaussian distribution, which can be created using the approach
of Kitaev and Webb (KW) kitaev2009wavefunction . To identify states of
definite multiplicity and momentum in $\ket{X}$ one can follow the general
ideas laid out in Jordan:2011ne ; Jordan:2011ci .

Our quantum circuit has been implemented in Qiskit Qiskit  and is available
from the authors upon request. In this first exploratory paper we compute the
foundational quantities, namely

$\displaystyle\mathcal{Y}_{X}=\absolutevalue{\\!\matrixelement{X}{\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]}{{\Omega}}}^{2}$
(14)

for $\ket{X}=\ket{{\Omega}}$ and $\ket{X}=\ket{p_{i}}$, the one-particle
momentum eigenstates of the theory. It should be noted that these quantities
are not infrared (IR) safe, and will therefore depend on the IR scale in the
problem, the lattice size $L$. However, as discussed in more detail in the
Supplemental Materials, there is no non-trivial IR safe observable that can be
defined in (1+1) dimensions, and these transition rates are therefore
representative quantities of what can be computed in this theory.

The quantum circuit for this measurement can be represented as

$\displaystyle\Qcircuit@C=0.5em@R=0.8em@!R$ /
$\displaystyle\lstick{\ket{\dotso}}$ $\displaystyle{/}\qw$
$\displaystyle\ghost{U_{\Omega}}$ $\displaystyle\qw$
$\displaystyle\ghost{U_{Y}}$ $\displaystyle\qw$
$\displaystyle\ghost{U_{X}^{\dagger}}$ $\displaystyle\qw$
$\displaystyle\lstick{\ket{l_{N-1}}}$ $\displaystyle{/}\qw$
$\displaystyle\ghost{U_{\Omega}}$ $\displaystyle\qw$
$\displaystyle\ghost{U_{Y}}$ $\displaystyle\qw$
$\displaystyle\ghost{U_{X}^{\dagger}}$ $\displaystyle\qw$

where $\ket{l_{n}}$ denotes the register of qubits for the $n^{\text{th}}$
lattice site. This creates the multivariate Gaussian vacuum state from the
initial state with all qubits zero using $U_{\Omega}$, acts on this vacuum
with the time ordered product of the two Wilson lines using $U_{Y}$, and
finally applies the inverse of the state preparation of state $\ket{X}$. The
details of these various circuits can be found in the Supplementary Material.

For our numerical results, we work with an $N=3$ site lattice with spacing
$\variation{x}=1$. With only 3 lattice sites, the Wilson line operator
simplifies to

$\displaystyle\mathcal{Y}_{X}=\absolutevalue{\\!\matrixelement{X}{\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]}{{\Omega}}}^{2}=\absolutevalue{\\!\matrixelement{X}{e^{ig\variation{x}\pqty{\phi_{x_{2}}-\phi_{x_{0}}}}}{{\Omega}}}^{2}\,,$
(15)

since all time evolution operators act on the initial or final eigenstate of
the Hamiltonian only and can therefore be neglected as contributing an overall
phase. In Fig. 1 we show the dependence of the expectation values
$\mathcal{Y}_{X}$ on the coupling $g$ for $n_{Q}=2$ qubits per lattice site
for different final states, and compare them against the analytical results,
shown by black lines. Results are given for both a quantum simulation and from
the 65-qubit IMBQ Manhattan quantum computer. The operators for implementing
all states are exact, as the resources for doing so on a small lattice are
modest. On a larger lattice approximate methods, such as KW ground state
approximation and the excited state preparation techniques of Jordan:2011ci
will be necessary; the effect of such approximations is presented in the
Supplementary Material.

Figure 1: Result of transition rates from the valuum of the Wilson line for 3
lattice sites and $n_{Q}=2$ qubits per site to the vacuum and the lowest-lying
single excited state. The solid lines shows the analytical result with no
field digitization while the dashed lines represents the result from a quantum
simulator of our circuit. The black data points show the result from the
65-qubit IBMQ Manhattan quantum computer, corrected both for readout errors
and CNOT gate errors. We only show result from the Manhattan computer for
$X={\Omega}$, since the circuit to measure the excited state was too deep to
give reliable results.

Errors in the quantum circuits, especially readout errors and CNOT gate errors
are quite large on existing hardware. As discussed in the Supplemental
Materials, the exponential of the field operator at a given position requires
only $n_{Q}$ single qubit gates, such that the operator in Eq. 15 requires no
CNOT gates. For $n_{Q}=2$, he state preparation requires 6 CNOT gates for
gates. Note that for more than 3 lattice sites the time evolution operator is
required, which requires a much larger amount of gates, although the resulting
circuits are known. For example, even for 3 lattice sites the standard
implementation of a single Trotter step of our Hamiltonian requires 60 CNOT
gates. We have applied both readout error mitigation as described in
ReadoutCorrection  as well as CNOT gate noise mitigation He:2020udd . For more
details, including References futureKW ; DAgostini:1994fjx ; 1974AJ…..79..745L
; Richardson:72 ; 2010.07496 ; geller_efficient_2020 ; song_10-qubit_2017 ;
gong_genuine_2019 ; hamilton2020scalable ; wei_verifying_2020 ; pyquil ; cirq
; arute2020quantum ; xacc ; alex2019xacc ; mccaskey_quantum_2019 ; ignis ;
1904.11935 ; 1907.08518 ; Dumitrescu:2018 ; PhysRevX.8.031027 ;
PhysRevLett.119.180509 , see the Supplemental Materials. One can see that the
digitized result with 2 qubits per lattice site differs from the analytic
calculation by up to $10\text{\,}\mathrm{\char 37\relax}$. This would be
reduced to at most $1.5\text{\,}\mathrm{\char 37\relax}$ by adding just a
single qubit per lattice site, since the resulting digitization errors fall
exponentially with the number of qubits. The quantum computer reproduces the
simulated result to about $5\text{\,}\mathrm{\char 37\relax}$ accuracy.

In conclusion, EFT are well known to be able to describe the low energy
dynamics of field theories and, through short distance, perturbatively
computable matching coefficients, can be used to describe the dynamics of a
full underlying quantum field theory. We have argued that the dynamics of a
low energy EFT can be simulated with significantly smaller quantum resources
than the dynamics of the corresponding full theory. In SCET the interactions
of highly energetic particles with soft particles of low energy are described
through operators containing Wilson lines, and we have shown in detail how the
dynamics of an analogous scalar soft theory can be described using quantum
algorithms. Using Wilson lines of free scalar fields in (1+1) dimension, we
have computed the simplest matrix elements in this soft theory, namely the
transition matrix elements from the vacuum to itself and the lowest-lying
excited states of two Wilson lines in opposite directions, using 3 lattice
sites. We have compared the computations on a quantum computer to analytical
results that can be obtained for this simple theory. Using only 2 qubits per
lattice we obtain results within $10\text{\,}\mathrm{\char 37\relax}$ of the
analytical result, and by using noise-mitigation techniques, uncertainties due
to running on present-day hardware can be reduced to about
$5\text{\,}\mathrm{\char 37\relax}$.

## Acknowledgements

We would like to thank Dorota Grabowska, Bert de Jong, Michael Kreshchuk, Pier
Monni, John Preskill, Martin Savage and Miro Urbanek for useful discussions.
CWB and BN are supported by the U.S. Department of Energy (DOE), Office of
Science under contract DE-AC02-05CH11231. In particular, support comes from
Quantum Information Science Enabled Discovery (QuantISED) for High Energy
Physics (KA2401032) and the Office of Advanced Scientific Computing Research
(ASCR) through the Accelerated Research for Quantum Computing Program. MF is
supported by the DOE under grant DE-SC0010008. This research used resources of
the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science
User Facility supported under Contract DE-AC05-00OR22725.

## References

  * (1) S. P. Jordan, K. S. M. Lee and J. Preskill, Science 336, 1130-1133 (2012) [arXiv:1111.3633 [quant-ph]].
  * (2) J. B. Kogut and L. Susskind, Phys. Rev. D 11, 395-408 (1975)
  * (3) J. B. Kogut, Rev. Mod. Phys. 51, 659 (1979)
  * (4) J. B. Bronzan, Phys. Rev. D 31, 2020-2028 (1985)
  * (5) S. P. Jordan, K. S. M. Lee and J. Preskill, Quant. Inf. Comput. 14, 1014-1080 (2014) [arXiv:1112.4833 [hep-th]].
  * (6) A. Macridin, P. Spentzouris, J. Amundson and R. Harnik, Phys. Rev. Lett. 121, no.11, 110504 (2018) [arXiv:1802.07347 [quant-ph]].
  * (7) A. Macridin, P. Spentzouris, J. Amundson and R. Harnik, Phys. Rev. A 98, no.4, 042312 (2018) [arXiv:1805.09928 [quant-ph]].
  * (8) N. Klco and M. J. Savage, Phys. Rev. A 99, no.5, 052335 (2019) [arXiv:1808.10378 [quant-ph]].
  * (9) D. C. Hackett, K. Howe, C. Hughes, W. Jay, E. T. Neil and J. N. Simone, Phys. Rev. A 99, no.6, 062341 (2019) [arXiv:1811.03629 [quant-ph]].
  * (10) K. Yeter-Aydeniz, E. F. Dumitrescu, A. J. McCaskey, R. S. Bennink, R. C. Pooser and G. Siopsis, Phys. Rev. A 99, no.3, 032306 (2019) [arXiv:1811.12332 [quant-ph]].
  * (11) M. Kreshchuk, W. M. Kirby, G. Goldstein, H. Beauchemin and P. J. Love, [arXiv:2002.04016 [quant-ph]].
  * (12) M. Kreshchuk, S. Jia, W. M. Kirby, G. Goldstein, J. P. Vary and P. J. Love, [arXiv:2009.07885 [quant-ph]].
  * (13) J. F. Haase, L. Dellantonio, A. Celi, D. Paulson, A. Kan, K. Jansen and C. A. Muschik, [arXiv:2006.14160 [quant-ph]].
  * (14) D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343-1346 (1973)
  * (15) H. D. Politzer, Phys. Rev. Lett. 30, 1346-1349 (1973)
  * (16) J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B 250, 199-224 (1985)
  * (17) S. Catani, L. Trentadue, G. Turnock and B. R. Webber, Nucl. Phys. B 407, 3-42 (1993)
  * (18) R. Bonciani, S. Catani, M. L. Mangano and P. Nason, Phys. Lett. B 575, 268-278 (2003) [arXiv:hep-ph/0307035 [hep-ph]].
  * (19) G. Marchesini and B. R. Webber, Nucl. Phys. B 238, 1-29 (1984)
  * (20) M. Bengtsson and T. Sjostrand, Nucl. Phys. B 289, 810-846 (1987)
  * (21) C. W. Bauer, W. A. De Jong, B. Nachman and D. Provasoli, [arXiv:1904.03196 [hep-ph]].
  * (22) E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean and P. Lougovski, Phys. Rev. Lett. 120, no.21, 210501 (2018) [arXiv:1801.03897 [quant-ph]].
  * (23) H. H. Lu, N. Klco, J. M. Lukens, T. D. Morris, A. Bansal, A. Ekström, G. Hagen, T. Papenbrock, A. M. Weiner and M. J. Savage, et al. Phys. Rev. A 100, no.1, 012320 (2019) doi:10.1103/PhysRevA.100.012320 [arXiv:1810.03959 [quant-ph]].
  * (24) M. J. Cervia, A. V. Patwardhan, A. B. Balantekin, S. N. Coppersmith and C. W. Johnson, Phys. Rev. D 100, no.8, 083001 (2019) [arXiv:1908.03511 [hep-ph]].
  * (25) A. Roggero, A. C. Y. Li, J. Carlson, R. Gupta and G. N. Perdue, Phys. Rev. D 101, no.7, 074038 (2020) [arXiv:1911.06368 [quant-ph]].
  * (26) M. J. Cervia, A. B. Balantekin, S. N. Coppersmith, C. W. Johnson, P. J. Love, C. Poole, K. Robbins and M. Saffman, [arXiv:2011.04097 [hep-th]].
  * (27) T. F. Stetina, A. Ciavarella, X. Li and N. Wiebe, [arXiv:2101.00111 [quant-ph]].
  * (28) C. W. Bauer, S. Fleming and M. E. Luke, Phys. Rev. D 63, 014006 (2000) [arXiv:hep-ph/0005275 [hep-ph]].
  * (29) C. W. Bauer, S. Fleming, D. Pirjol and I. W. Stewart, Phys. Rev. D 63, 114020 (2001) [arXiv:hep-ph/0011336 [hep-ph]].
  * (30) C. W. Bauer and I. W. Stewart, Phys. Lett. B 516, 134-142 (2001) [arXiv:hep-ph/0107001 [hep-ph]].
  * (31) C. W. Bauer, D. Pirjol and I. W. Stewart, Phys. Rev. D 65, 054022 (2002) [arXiv:hep-ph/0109045 [hep-ph]].
  * (32) C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothstein and I. W. Stewart, Phys. Rev. D 66, 014017 (2002) [arXiv:hep-ph/0202088 [hep-ph]].
  * (33) C. W. Bauer, A. V. Manohar and M. B. Wise, Phys. Rev. Lett. 91, 122001 (2003) [arXiv:hep-ph/0212255 [hep-ph]].
  * (34) C. W. Bauer, C. Lee, A. V. Manohar and M. B. Wise, Phys. Rev. D 70, 034014 (2004) [arXiv:hep-ph/0309278 [hep-ph]].
  * (35) A. V. Manohar, Phys. Rev. D 68, 114019 (2003) [arXiv:hep-ph/0309176 [hep-ph]].
  * (36) T. Becher and M. Neubert, Phys. Lett. B 637, 251-259 (2006) [arXiv:hep-ph/0603140 [hep-ph]].
  * (37) R. Goerke and M. Luke, JHEP 02, 147 (2018) [arXiv:1711.09136 [hep-ph]].
  * (38) T. Becher and M. Neubert, Phys. Lett. B 633, 739-747 (2006) [arXiv:hep-ph/0512208 [hep-ph]].
  * (39) A. H. Hoang and S. Kluth, [arXiv:0806.3852 [hep-ph]].
  * (40) T. T. Jouttenus, I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, Phys. Rev. D 83, 114030 (2011) [arXiv:1102.4344 [hep-ph]].
  * (41) R. Kelley, M. D. Schwartz, R. M. Schabinger and H. X. Zhu, Phys. Rev. D 84, 045022 (2011) [arXiv:1105.3676 [hep-ph]].
  * (42) P. F. Monni, T. Gehrmann and G. Luisoni, JHEP 08, 010 (2011) [arXiv:1105.4560 [hep-ph]].
  * (43) R. Boughezal, X. Liu and F. Petriello, Phys. Rev. D 91, no.9, 094035 (2015) [arXiv:1504.02540 [hep-ph]].
  * (44) I. Moult and H. X. Zhu, JHEP 08, 160 (2018) [arXiv:1801.02627 [hep-ph]].
  * (45) C. Lin, F. H. Zong and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001) [arXiv:cond-mat/0101339 [cond-mat.stat-mech]].
  * (46) C. T. Sachrajda and G. Villadoro, Phys. Lett. B 609, 73-85 (2005) [arXiv:hep-lat/0411033 [hep-lat]].
  * (47) P. F. Bedaque, Phys. Lett. B 593, 82-88 (2004) [arXiv:nucl-th/0402051 [nucl-th]].
  * (48) R. A. Briceno, Z. Davoudi, T. C. Luu and M. J. Savage, Phys. Rev. D 89, no.7, 074509 (2014) [arXiv:1311.7686 [hep-lat]].
  * (49) H. F. Trotter, Proceedings of the American Mathematical Society 10 (1959) 545.
  * (50) Masuo Suzuki, Communications in Mathematical Physics 51 (1976) 183.
  * (51) Masuo Suzuki, Progress of Theoretical Physics 56 (1976) 1454.
  * (52) A. Kitaev and W. A. Webb [arXiv:0801.0342 [quant-ph]].
  * (53) H. Abraham et. al. 10.5281/zenodo.2562110
  * (54) B. Nachman, M. Urbanek, W. A. de Jong, C. W. Bauer, npj Quantum Inf 6, 84 (2020) [arXiv:1910.01969 [quant-ph]].
  * (55) A. He, B. Nachman, W. A. de Jong and C. W. Bauer, Phys. Rev. A 102, no.1, 012426 (2020) [arXiv:2003.04941 [quant-ph]].
  * (56) C. W. Bauer, P. Deliyannis, M. Freytsis, B. Nachman, _forthcoming_.
  * (57) G. D’Agostini, Nucl. Instrum. Methods Phys. Res. A 362 (1995) 487.
  * (58) L. B. Lucy, Astron. J. 79 (1974) 745.
  * (59) W. H. Richardson, J. Opt. Soc. Am. 62 (1972) 55.
  * (60) R. Hicks, C. Bauer, and B. Nachman, Phys. Rev. A 103, 022407 (2021) [arXiv:2010.07496 [quant-ph]].
  * (61) M. Geller and M. Sun, [arXiv:2001.09980 [quant-ph]].
  * (62) C. Song et al., Phys. Rev. Lett. 119 (2017) 180511.
  * (63) M. Gong et al., Phys. Rev. Lett 122 (2019) 110501.
  * (64) K. E. Hamilton et al., [arXiv:2006.01805 [quant-ph]].
  * (65) K. X. Wei et al., Phys. Rev. A 101 (2020) 032343, [arXiv:1905.05720 [quant-ph]].
  * (66) Rigetti Forest, http://docs.rigetti.com, 2020.
  * (67) The Cirq Contributors, https://github.com/quantumlib/Cirq, 2020.
  * (68) F. Arute et al., [arXiv:2004.04197 [quant-ph]].
  * (69) The XACC Contributors, https://xacc.readthedocs.io, 2020.
  * (70) A. J. McCaskey et al., [arXiv:1911.02452 [quant-ph]].
  * (71) A. J. McCaskey et al., npj Quantum Information 5 (2019) 99.
  * (72) IBM Research, https://qiskit.org/ignis, 2019.
  * (73) Y. Chen, M. Farahzad, S. Yoo, and T. Wei, [arXiv:1904.11935 [quant-ph]].
  * (74) F. B. Maciejewski, Z. Zimboras, and M. Oszmaniec, [arXiv:1907.08518 [quant-ph]].
  * (75) E. F. Dumitrescu et al., Phys. Rev. Lett. 120 (2018) 210501.
  * (76) S. Endo, S. C. Benjamin, and Y. Li, Phys. Rev. X 8 (2018) 031027.
  * (77) K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119 (2017) 180509.

Simulating effective field theories on quantum computers

Supplementary Material Christian W. Bauer, Marat Freytsis, Benjamin Nachman

##  I Analytical calculations in the lattice field theory

In this appendix we provide analytical calculations for the field theory under
considerations that our quantum algorithms can be compared against. The
lattice field theory and its quantization agree with the analogous treatment
in Jordan:2011ci , while the calculations involving the Wilson lines have not
appeared anywhere else to our knowledge.

###  I.1 Definition of the lattice

We consider a lattice of $d$ dimensions, $N$ lattice sites for each dimension
and lattice spacing $\variation{x}$. Points on the latice are indexed by $d$
integers ranging from $0$ to $N-1$. We denote this set of integers by

$\displaystyle\mathbf{r}=(r_{1},\dotsc,r_{d})\,,\quad
r_{i}\in\\{0,\dotsc,N-1\\}$ (S1)

with the the positions of the lattice in each dimension given by

$\displaystyle\mathbf{x}_{\mathbf{r}}=-\mathbf{x}_{\textrm{max}}+\mathbf{r}\variation{x}\,,\quad\mathbf{x}_{\textrm{max}}=(n_{0}\variation{x},\dotsc,n_{0}\variation{x})\,\mbox{\quad
using\quad}n_{0}=\left\lfloor\frac{N-1}{2}\right\rfloor\,.$ (S2)

The corresponding momentum values in each dimension are then given by
($\mathbf{s}$ denotes the same set of integers as $\mathbf{r}$)

$\displaystyle\mathbf{p}_{\mathbf{s}}=-\mathbf{p}_{\mathrm{max}}+(\mathbf{s}-{\bf\it\Delta})\variation{p}\,,\quad\mathbf{p}_{\mathrm{max}}=(n_{0}\variation{p},\dotsc,n_{0}\variation{p})\,\mbox{\quad
using\quad}\variation{p}=\frac{2\pi}{N\variation{x}}\equiv\frac{2\pi}{L}\,.$
(S3)

Periodic boundary conditions correspond to ${\bf\it\Delta}=(0,\dotsc,0)$,
while twisted boundary conditions (which are equal to periodic boundary
conditions with an extra factor of $(-1)$ for each wrap around the lattice),
correspond to ${\bf\it\Delta}=(1/2,\dotsc,1/2)$.

To denote the sum over all lattice positions or conjugate lattice momenta we
define

$\displaystyle\sum_{\mathbf{r}}\equiv\sum_{u=1}^{d}\sum_{r_{u}=0}^{N-1}\,.$
(S4)

Note that the twisted boundary conditions ensure that the mode with vanishing
momentum is not part of the spectrum. For each dimension every momentum has a
partner with negative momentum

$\displaystyle p_{m}=-p_{2n_{0}-m+1}\,\mbox{\quad for\quad}m>0\,,$ (S5)

except for an odd number of lattice sites, in which case the momentum $p_{0}$
does not have corresponding partner. (Note that one can of course translate
the lattice and conjugate lattice such that
$\mathbf{x}_{\mathbf{r}}=\mathbf{r}\variation{x}$ and
$\mathbf{p}_{\mathbf{s}}=(\mathbf{s}+{\bf\it\Delta})\variation{p}$.)

To go between fields defined on the position spaced lattice to those defined
on the momentum spaced lattice one uses the discrete Fourier transform which
is given by

$\displaystyle\phi_{\mathbf{x}}=\frac{1}{(2\pi)^{d}}\sum_{\mathbf{p}}e^{-i\mathbf{p}\cdot\mathbf{x}}\phi_{\mathbf{p}}\,,\quad\phi_{\mathbf{p}}=\sum_{\mathbf{x}}e^{i\mathbf{p}\cdot\mathbf{x}}\phi_{\mathbf{x}}\,.$
(S6)

Note that we have introduced a short-hand form for the sum over lattice sites
as

$\displaystyle\sum_{\mathbf{x}}\equiv(\variation{x})^{d}\sum_{\mathbf{r}}\,,\quad\sum_{\mathbf{p}}\equiv(\variation{p})^{d}\sum_{\mathbf{s}}\,,$
(S7)

where $\mathbf{r}$ and $\mathbf{s}$ denote the positions on the regular and
conjugate lattice, respectively.

One has the usual completeness relation

$\displaystyle\frac{1}{(2\pi)^{d}}\sum_{p}e^{i\mathbf{p}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}=\delta_{\mathbf{x}\mathbf{x}^{\prime}}\,,\quad\sum_{x}e^{i\mathbf{x}\cdot(\mathbf{p}-\mathbf{p}^{\prime})}=(2\pi)^{d}\delta_{\mathbf{p}\mathbf{p}^{\prime}}\,,$
(S8)

were we have again used a short hand notation

$\displaystyle\delta_{\mathbf{x}\mathbf{x}^{\prime}}\equiv\frac{1}{(\delta
x)^{d}}\delta_{\mathbf{r}\mathbf{r}^{\prime}}\,,\quad\delta_{\mathbf{p}\mathbf{p}^{\prime}}\equiv\frac{1}{(\delta
p)^{d}}\delta_{\mathbf{s}\mathbf{s}^{\prime}}\,,$ (S9)

and we have chosen the factor of $2\pi$ to reflect the standard choice made in
continuum field theory.

Before we conclude this section, we briefly discuss the relationship between
the boundary conditions and the value of ${\bf\it\Delta}$. From the definition
of the Fourier transformation in Eq. S6, one can immediately see that
transforming a field by the lattice size in some dimension yields

$\displaystyle\phi_{\mathbf{x}+L\mathbf{n}_{i}}=e^{i(2\pi{\bf\it\Delta}_{i})}\phi_{\mathbf{x}}\,.$
(S10)

Thus, non-integer values of ${\bf\it\Delta}_{i}$ give an extra phase in the
boundary condition, with half-integer values adding a relative minus sign,
giving twisted boundary conditions.

###  I.2 The free field theory

The free scalar field theory on the lattice is given by the Hamiltonian

$\displaystyle
H=\frac{1}{2}\sum_{\mathbf{x}}\bqty{\dot{\phi}_{\mathbf{x}}^{2}-\phi_{\mathbf{x}}\,[\gradient^{2}\,\phi]_{\mathbf{x}}}\,,$
(S11)

where the sums run over the $n$ possible integers (lattice sites) for each
dimension. The d’Alembertian operator is defined through its action on a
field, analogously to the second derivative operator in the main text.
Performing a Fourier transform, one can write this expression as

$\displaystyle
H=\frac{1}{2}\frac{1}{(2\pi)^{d}}\sum_{\mathbf{p}}\bqty{\dot{\phi}_{\mathbf{p}}\dot{\phi}_{-\mathbf{p}}+\omega_{\mathbf{p}}^{2}\,\phi_{\mathbf{p}}\phi_{-\mathbf{p}}}\,,$
(S12)

where frequencies are given by

$\displaystyle\omega_{\mathbf{p}}=\frac{2}{\variation{x}}\sqrt{\sum_{i=1}^{d}\sin^{2}\pqty{\frac{\mathbf{p}_{i}\variation{x}}{2}}}\,.$
(S13)

Introducing the creation and annihilation operators as

$\displaystyle\dot{\phi}_{\mathbf{p}}=-i\sqrt{\frac{\omega_{\mathbf{p}}}{2}}\pqty{a_{\mathbf{p}}-a^{\dagger}_{-\mathbf{p}}}\,,\quad\phi_{\mathbf{p}}=\frac{1}{\sqrt{2\omega_{\mathbf{p}}}}\pqty{a_{\mathbf{p}}+a^{\dagger}_{-\mathbf{p}}}\,,$
(S14)

which satisfy the commutation relations

$\displaystyle\commutator{a_{\mathbf{p}}}{a^{\dagger}_{\mathbf{q}}}=\pqty{\frac{2\pi}{\variation{p}}}^{d}\delta_{\mathbf{p}\mathbf{q}}\,,$
(S15)

one finds for the Hamiltonian

$\displaystyle
H=\frac{1}{(2\pi)^{d}}\sum_{\mathbf{p}}\omega_{\mathbf{p}}\bqty{a_{\mathbf{p}}a^{\dagger}_{\mathbf{p}}+\frac{1}{2}\pqty{\frac{2\pi}{\variation{p}}}^{d}}\,.$
(S16)

Thus, in momentum space the free field theory on a lattice is simply a
collection of uncoupled harmonic oscillators, and one can therefore compute
the spectrum of the theory rather easily. A general state is therefore labeled
by its $d\times n$ occupation numbers
$\mathbf{l}^{(0)}\dotsb\mathbf{l}^{(n-1)}$, where the occupation number at
each lattice site is determined by $d$ integers $\geq 0$. Thus, a general
state is given by

$\displaystyle\ket{\Psi_{l_{\mathbf{k}}}}=\bigotimes_{l_{\mathbf{k}}}\ket{l_{\mathbf{k}}}\,,$
(S17)

with the ground state given by

$\displaystyle\ket{{\Omega}}=\bigotimes_{l_{\mathbf{k}}}\ket{0}\,.$ (S18)

The energy of the ground state $\ket{{\Omega}}$ is given by

$\displaystyle H\ket{{\Omega}}=E_{0}\ket{{\Omega}}\,,\quad
E_{0}=\frac{1}{2}\sum_{\mathbf{s}}\omega_{\mathbf{s}}\,,$ (S19)

while the energy of excited states are easily determined in terms of their
occupation numbers as

$\displaystyle
H\ket{\Psi_{l_{\mathbf{k}}}}=E_{l_{\mathbf{k}}}\ket{\Psi_{l_{\mathbf{k}}}}\,,\quad
E_{l_{\mathbf{k}}}=E_{0}+\sum_{\mathbf{s}}l_{\mathbf{s}}\,\omega_{\mathbf{s}}$
(S20)

###  I.3 Wilson Line expectation values

Using this notation, one can compute the expectation values of the Wilson line
operator $T[Y_{n}Y_{\bar{n}}^{\dagger}]$. The Wilson lines originate from the
point $x=0$, which is labeled by the lattice position $n_{0}$, defined in Eq.
S2. Note that with this notation $x_{k}$ with $k<n_{0}$ are negative, while
$x_{k}$ with $k>n_{0}$ are positive. Note that in this section we assume that
we have an odd number of lattice sites, such that we have as many lattice
points with $x>0$ as with $x<0$. This means that $2n_{0}=N-1$.

As already discussed in the main paper, the Wilson lines on the lattice can be
written as

$\displaystyle\begin{split}Y_{n}&=\mathrm{P}\exp\bqty{ig\variation{x}\sum_{i=n_{0}}^{N-1}\phi_{\mathbf{n}x_{i}}(t=x_{i}-n_{0})}\,,\\\
Y_{{\bar{n}}}^{\dagger}&=\mathrm{P}\exp\bqty{-ig\variation{x}\sum_{i=0}^{n_{0}}\phi_{\mathbf{n}x_{i}}(t=n_{0}-x_{i})}\,.\end{split}$
(S21)

The Wilson line operator is therefore given by

$\displaystyle\begin{split}\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]&=e^{-iH\variation{x}\,n_{0}}\exp\bqty{ig\variation{x}\pqty{\phi_{\mathbf{n}x_{2n_{0}}}-\phi_{\mathbf{n}x_{0}}}}\\\
&\qquad\times
e^{iH\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{\mathbf{n}x_{2n_{0}-1}}-\phi_{\mathbf{n}x_{1}}}}\\\
&\qquad\times\dotsm\times
e^{iH\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{\variation{n}x_{n_{0}+1}}-\phi_{\variation{n}x_{n_{0}-1}}}}\\\
&\qquad\times
e^{iH\variation{x}}\exp\bqty{ig\variation{x}\pqty{\phi_{\variation{n}x_{n_{0}}}-\phi_{\variation{n}x_{n_{0}}}}}\\\
&\equiv U_{n_{0}}(-t)\,U_{n_{0}-1}(\variation{t})\dotsb
U_{1}(\variation{t})U_{0}(\variation{t})\,.\end{split}$ (S22)

where we have used that $\mathbf{{\bar{n}}}=-\mathbf{n}$ and
$x_{n_{0}+m}=-x_{n_{0}-m}$. In the last line we defined

$\displaystyle U_{m}(t)\equiv
e^{iHt}\exp\bqty{ig\variation{x}\pqty{\phi_{\mathbf{n}x_{n_{0}+m}}-\phi_{\mathbf{n}x_{n_{0}-m}}}}\,,$
(S23)

as well as $t\equiv n_{0}\variation{x}$, $\variation{t}\equiv\variation{x}$.
Note that since $x_{n_{0}}=0$ one has $U_{0}(t)=e^{-iHt}$.

Using Eq. S14 we can write

$\displaystyle\begin{split}U_{m}(t)&=e^{iHt}\prod_{\mathbf{s}}\exp\\!\bqty{ig\variation{x}\pqty{\frac{\variation{p}}{2\pi}}^{d}\pqty{e^{i\mathbf{n}\cdot\mathbf{p}_{\mathbf{s}}m\variation{x}}-e^{-i\mathbf{n}\cdot\variation{p}_{\mathbf{s}}m\variation{x}}}\frac{1}{\sqrt{2\omega_{\mathbf{p}_{\mathbf{s}}}}}\pqty{a^{\dagger}_{\mathbf{p}_{\mathbf{s}}}-a_{\mathbf{p}_{\mathbf{s}}}}}\\\
&=e^{iHt}\prod_{\mathbf{s}}\exp\\!\bqty{-2g\variation{x}\pqty{\frac{\variation{p}}{2\pi}}^{d}\sqrt{\frac{1}{2\omega_{\mathbf{p}_{\mathbf{s}}}}}\sin(\mathbf{n}\cdot\mathbf{p}_{\mathbf{s}}\,m\variation{x})\pqty{a^{\dagger}_{\mathbf{p}_{\mathbf{s}}}-a_{\mathbf{p}_{\mathbf{s}}}}}\\\
&=e^{iHt}\prod_{\mathbf{s}}D_{\mathbf{p}_{\mathbf{s}}}\pqty{\alpha_{\mathbf{p}_{\mathbf{s}}}(m)}\,,\end{split}$
(S24)

with

$\displaystyle\alpha_{\mathbf{p}_{\mathbf{s}}}(m)=2g\variation{x}\pqty{\frac{\variation{p}}{2\pi}}^{d}\sqrt{\frac{1}{2\omega_{\mathbf{p}_{\mathbf{s}}}}}\sin(\mathbf{n}\cdot\mathbf{p}_{\mathbf{s}}\,m\variation{x})\,.$
(S25)

In the last line we have written the result in terms of the displacement
operator that is well know from the theory of coherent states

$\displaystyle
D_{\mathbf{p}}(\alpha_{\mathbf{p}})=e^{\alpha_{\mathbf{p}}a_{\mathbf{p}}^{\dagger}-\alpha_{\mathbf{p}}^{*}a_{\mathbf{p}}}\,,$
(S26)

which satisfies the relation

$\displaystyle
D_{\mathbf{p}}(\alpha_{\mathbf{p}})D_{\mathbf{p}}(\beta_{\mathbf{p}})=(2\pi)^{d}\,D_{\mathbf{p}}(\alpha_{\mathbf{p}}+\beta_{\mathbf{p}})\,e^{i\imaginary(\alpha_{\mathbf{p}}\beta^{*}_{\mathbf{p}})\pqty{\frac{2\pi}{\variation{p}}}^{d}}\,,$
(S27)

and displacement operators acting on different momentum modes commute

$\displaystyle[D_{\mathbf{p}}(\alpha_{\mathbf{p}}),D_{\mathbf{q}}(\beta_{\mathbf{q}})]=0\,.$
(S28)

The action of the time evolution on the displacement operator is

$\displaystyle
e^{iHt}D_{\mathbf{p}}[\alpha_{\mathbf{p}}]=D[\alpha_{\mathbf{p}}(t)]\,\mbox{\quad
where\quad}\alpha_{\mathbf{p}}(t)=\alpha_{\mathbf{p}}\,e^{i\omega_{\mathbf{p}}t}\,.$
(S29)

Combining all this information, one finds

$\displaystyle\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]$
$\displaystyle=\prod_{\mathbf{s}}D_{\mathbf{p}_{\mathbf{s}}}\bqty{\sum_{m=0}^{n_{0}}\alpha_{\mathbf{p}_{\mathbf{s}}}(m)e^{i(n_{0}-m)\omega_{s}\variation{t}}}e^{i\phi_{\mathbf{p}_{\mathbf{s}}}}\,,$
(S30)

with

$\displaystyle\phi_{\mathbf{p}_{\mathbf{s}}}=\pqty{\frac{2\pi}{\variation{p}}}^{d}\sum_{n>m}\imaginary\bqty{\alpha_{\mathbf{p}_{\mathbf{s}}}(n)\,\alpha^{*}_{\mathbf{p}_{\mathbf{s}}}(m)}\,.$
(S31)

To compute expectation values of this operator, we use

$\displaystyle\matrixelement{n_{\mathbf{p}}}{D_{\mathbf{p}}[\alpha_{\mathbf{p}}]}{{\Omega}}=\exp\\!\bqty{-\pqty{\frac{2\pi}{\variation{p}}}^{d}\frac{\absolutevalue{\alpha_{\mathbf{p}}}^{2}}{2}}\frac{\alpha_{\mathbf{p}}^{n_{p}}}{\sqrt{n_{\mathbf{p}}!}}\,.$
(S32)

In particular, the magnitude square of the vacuum expectation value is given
by

$\displaystyle\begin{split}\absolutevalue{\\!\matrixelement{{\Omega}}{\mathrm{T}[Y_{n}\,Y_{\bar{n}}^{\dagger}}{{\Omega}}}^{2}&=\prod_{\mathbf{s}}\exp\\!\bqty{-\pqty{\frac{2\pi}{\variation{p}}}^{d}\absolutevalue{\sum_{m=0}^{n_{0}}\alpha_{\mathbf{p_{s}}}(m)e^{i(n_{0}-m)\omega_{\mathbf{p}_{\mathbf{s}}}\variation{t}}}^{2}}\\\
&=\prod_{\mathbf{s}}\exp\\!\bqty{-4g^{2}(\variation{x})^{2}\pqty{\frac{2\pi}{\variation{p}}}^{d}\frac{1}{2\omega_{\mathbf{p}_{\mathbf{s}}}}\left|\sum_{m=0}^{n_{0}}e^{i(n_{0}-m)\omega_{\mathbf{p}_{\mathbf{s}}}\delta
t}\sin\left(\mathbf{n\cdot p_{s}}\,m\,\delta x\right)\right|^{2}}\\\
&=\exp\\!\bqty{-4\frac{g^{2}}{(2\pi)^{d}}\sum_{\mathbf{p}}\frac{1}{2\omega_{\mathbf{p}}}\absolutevalue{\sum_{x}e^{-i\omega_{\mathbf{p}}x}\sin(\mathbf{n}\cdot\mathbf{p}\,x)}^{2}}\\\
&=\exp\\!\bqty{-8\frac{g^{2}}{(2\pi)^{d}}\sum_{\mathbf{p}}\frac{1}{2\omega_{\mathbf{p}}}\sum_{x\geq
y}\cos(\omega_{\mathbf{p}}(x-y))\sin(\mathbf{n}\cdot\mathbf{p}\,x)\sin(\mathbf{n}\cdot\mathbf{p}\,y)}\,,\end{split}$
(S33)

where the sum over $x$ and $y$ runs over the lattice position and we have
defined

$\displaystyle\sum_{x\geq
y}f(x,y)\equiv\frac{1}{2}\sum_{x}f(x,x)+\sum_{x>y}f(x,y)\,.$ (S34)

One can verify that this expression reproduces the perturbative result to
order $g^{2}$. To that order, one can write

$\displaystyle\matrixelement{{\Omega}}{\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]}{{\Omega}}$
$\displaystyle=1-\frac{g^{2}}{2}\int_{0}^{L}\differential{s}\int_{0}^{L}\differential{t}\bqty{D_{F}(\mathbf{n}s-\mathbf{n}t)+D_{F}(\mathbf{{\bar{n}}}s-\mathbf{{\bar{n}}}t)-D_{F}(\mathbf{n}s-\mathbf{{\bar{n}}}t)-D_{F}(\mathbf{{\bar{n}}}s-\mathbf{n}t)}$
$\displaystyle=1-g^{2}\int_{0}^{L}\differential{s}\int_{0}^{s}\differential{t}\bqty{D(\mathbf{n}s-\mathbf{n}t)+D(\mathbf{{\bar{n}}}s-\mathbf{{\bar{n}}}t)-D(\mathbf{n}s-\mathbf{{\bar{n}}}t)-D(\mathbf{{\bar{n}}}s-\mathbf{n}t)}$
$\displaystyle=1-g^{2}\int_{0}^{L}\differential{s}\int_{0}^{s}\differential{t}\int\frac{\differential[d]{p}}{(2\pi)^{d}}\frac{1}{2\omega_{\mathbf{p}}}e^{-i\omega_{\mathbf{p}}(s-t)}\bqty{e^{i\mathbf{n}\cdot\mathbf{p}(s-t)}+e^{-i\mathbf{n}\cdot\mathbf{p}(s-t)}-e^{i\mathbf{n}\cdot\mathbf{p}(s+t)}-e^{-i\mathbf{n}\cdot\mathbf{p}(s+t)}}$
$\displaystyle=1-2g^{2}\int_{0}^{L}\differential{s}\int_{0}^{s}\differential{t}\int\frac{\differential[d]{p}}{(2\pi)^{d}}\frac{1}{2\omega_{\mathbf{p}}}e^{-i\omega_{\mathbf{p}}(s-t)}\bqty{\cos(\mathbf{n}\cdot\mathbf{p}(s-t))-\cos(\mathbf{n}\cdot\mathbf{p}(s+t))}$
$\displaystyle=1-4g^{2}\int_{0}^{L}\differential{s}\int_{0}^{s}\differential{t}\int\frac{\differential[d]{p}}{(2\pi)^{d}}\frac{1}{2\omega_{\mathbf{p}}}e^{-i\omega_{\mathbf{p}}(s-t)}\sin(\mathbf{n}\cdot\mathbf{p}\,s)\sin(\mathbf{n}\cdot\mathbf{p}\,t)$
(S35)

Taking the square of this result one finds

$\displaystyle\absolutevalue{\\!\matrixelement{{\Omega}}{\mathrm{T}[Y_{n}\,Y_{{\bar{n}}}^{\dagger}]}{{\Omega}}}^{2}=1-8g^{2}\int_{0}^{L}\differential{s}\int_{0}^{s}\differential{t}\int\frac{\differential[d]{p}}{(2\pi)^{d}}\frac{1}{2\omega_{\mathbf{p}}}\cos(\omega_{\mathbf{p}}(s-t))\sin(\mathbf{n}\cdot\mathbf{p}\,s)\sin(\mathbf{n}\cdot\mathbf{p}\,t)$
(S36)

which is the continuum limit of Eq. S33.

##  II The effect of digitization

In this section we provide some details of the effect of digitization on the
Hamiltonian of a free massless field theory. For this, we consider a (1+1)
dimensional field theory with 3 lattice sites at points

$\displaystyle x_{0}=-\variation{x}\,,\quad x_{1}=0\,,\quad
x_{2}=\variation{x}\,,$ (S37)

with Hamiltonian

$\displaystyle
H=\frac{\variation{x}}{2}\sum_{i=0}^{2}\bqty{\dot{\phi}_{i}^{2}-\phi_{i}\,[\nabla^{2}\phi]_{j}}\,,$
(S38)

The momenta of the 3-site lattice are given by

$\displaystyle p_{0}=-\frac{\pi}{3\variation{x}}\,,\quad
p_{1}=-\frac{\pi}{\variation{x}}\,,\quad p_{2}=\frac{\pi}{\variation{x}}\,.$
(S39)

The corresponding frequencies are obtained from Eq. 12

$\displaystyle\omega_{0}=\frac{2}{\variation{x}}\,,\quad\omega_{1}=\omega_{2}=\frac{1}{\variation{x}}\,.$
(S40)

The exact eigenvalues of this Hamiltonian are now easily obtained. The energy
of the ground state is given by half of the sum of frequencies, and therefore
$E_{0}=2/\variation{x}$. The excited states are simply given by the ground
state energy increased by up to $n_{\phi}$ additions of each frequency. Thus,
the lowest lying states are

$\displaystyle 1\times\frac{2}{\variation{x}}\,,\quad
2\times\frac{3}{\variation{x}}\,,\quad 4\times\frac{4}{\variation{x}}\,,\quad
6\times\frac{5}{\variation{x}}\,,\quad\dotsc$ (S41)

We compare the results obtained for the 15 lowest lying states with the
eigenvalues of the above Hamiltonian using $n_{q}=2$, $n_{q}=3$ and $n_{q}=4$
qubits per lattice site in Fig. 1. One can see that for $n_{q}=4$ these lowest
states are reproduced with an accuracy of ${10}^{-7}$, for $n_{q}=3$ with an
accuracy of ${10}^{-3}$ and for $n_{q}=2$ with an accuracy of about
$10\text{\,}\mathrm{\char 37\relax}$. This is in agreement with the findings
of Klco:2018zqz .

Figure 1: The fractional error in the 15 lowest-lying eigenstates on an $N=3$
lattice. The dependence on the number of qubits per lattice site is shown for
$n_{Q}=2$ (orange), $n_{Q}=3$ (green), and $n_{Q}=4$ (blue). The upper (lower)
part of the plot displays eigenvalues greater (smaller) than their undigitized
values. For $n_{Q}=3$, all 2 particles states are already reproduced at better
than $1\text{\,}\mathrm{\char 37\relax}$ accuracy. (If modes with
$\gtrsim$10\text{\,}\mathrm{\char 37\relax}$$ deviations from their continuum
values are excluded, this increases to all 5 particle states.) The exponential
improvement with number of qubits per site is clearly visible.  Figure 2: The
effect of field digitization on transition rates in the presence of the Wilson
line for $N=3$. The blue solid lines are the analytical result for transitions
to the vacuum (${\Omega}$) and the lowest-lying one one-particle state
($p_{1}$). Digitization limits the number of states the field can take on for
$n_{Q}=1$ (red) to 2,$n_{Q}=2$ (orange) to 4 and $n_{Q}=3$ green to 8. In each
case, the exact ground state is digitized and evolved.

Figure 3: The effect of approximated state preparations on transition rates in
the presence of the Wilson line for $N=3$ lattice sites. The blue solid lines
are the analytical result for transitions to the vacuum (${\Omega}$) and the
lowest-lying one one-particle state ($p_{1}$). We compare the digitization on
the exact state with the original and modified KW approximation (see text for
details) which can be prepared with polynomial resources in number of qubits.
The plots show results for $n_{Q}=2$ (left) and $n_{Q}=3$ (right).

In Fig. 2 we show the dependence of the vacuum expectation value
$\mathcal{Y}_{X}$ on the coupling $g$ on the number of qubits per lattice
site. One can see that already for $n_{Q}=2$ qubits per lattice site one gets
a very good approximation to the exact lattice result, and for $n_{Q}=3$ the
results including digitization are indistinguishable from the exact results if
using the exact correlations between sites in the ground state. Figure 3 shows
the dependence on the ground state, with the exact ground state shown by the
solid squares and the result for the KW approximation shown by stars for
$n_{Q}=2$ and 3. One can see that using the KW approximation introduces an
additional error, with the resulting emission probability reduced compared to
the exact result, but that the resulting difference is sub-percent level by
the time $n_{Q}=3$. For these plots only, we compare the KW approximation
originally presented in kitaev2009wavefunction  and used in the rest of the
text with a modified form discussed in Section IV.3.

##  III Consequences of the non-local nature of Wilson lines special
considerations in (1+1) dimensions

The two Wilson lines $Y_{n}$ and $Y_{\bar{n}}^{\dagger}$ describe how the
dynamical soft degrees of freedom in the effective theory (bosons in our case)
interact with the fermions that move at the speed of light through the system.
These Wilson lines contain an integral over the scalar fields along the world
lines of the fermions, going from point 0 to $\infty$ along the directions
$n^{\mu}=(1,\mathbf{n})$ and ${\bar{n}}^{\mu}=(1,-\mathbf{n})$, and are
therefore non-local in nature, with the non-locality extending all the way to
infinity. This non-locality along the light-like directions gives rise to the
well known collinear singularities in the matrix elements. An important
consequence of the non-local nature of the operators and their infinite extent
is that their lattice implementation necessarily makes them directly sensitive
to the lattice volume $L$ and separation $\variation{x}$. The fact that the
Wilson lines can only extend up to the edge of the lattice provides a
regulator for the angle between the momentum of bosons and the direction
$\mathbf{n}$, in addition to the IR regulator on the energy of each boson.

As we will now show, this fact gives rise to mixed IR–UV divergences in the
lattice definition of the Wilson line operators. The presence of the finite
lattice produces terms $L\,n\cdot p$ and $L\,{\bar{n}}\cdot p$, regulating
divergences as $n\cdot p$ and ${\bar{n}}\cdot p$ tend to zero. On the other
hand, for finite lattice spacing one has
$\omega-\absolutevalue{\mathbf{p}}\sim\absolutevalue{\mathbf{p}}^{2}\variation{x}^{2}$,
such that the UV regulator is also prohibiting the momenta of the bosons to go
on-shell, which in turn also regulates the collinear divergence. This implies
that collinear divergences in matrix elements give rise to mixed UV-IR
divergences, which would be absent in the continuum theory (even though the
continuum theory also contains mixed UV-IR divergences coming from the the
limit $n\cdot p\to 0$ with $p\to\infty$). Once infrared and collinear (IRC)
safe physical observables are being calculated, in which all collinear
divergences cancel, these lattice artifacts will also cancel, leaving only
pure UV divergences related to the normalization of operators.

There is an additional subtlety in a (1+1) dimensional field theory as is used
in the numerical results of this paper. With only a single spatial dimension
one can of course only define a single direction. This means that in such a
theory all momenta have to be aligned with the direction $n$ or ${\bar{n}}$,
which means that every momentum satisfies $n\cdot p=0$ or ${\bar{n}}\cdot p=0$
and therefore gives rise to a collinear divergence. In other words, each state
of the quantum field theory is actually collinear to the directions $n^{\mu}$
and ${\bar{n}}^{\mu}$. It should be recalled that the underlying short
distance process we are describing is the production of two energetic back-to-
back particles. The differential distribution in the total energy of this
process is described by short distance physics captured by the hard function,
while the jet function and the soft function discussed in this paper describe
the distribution of collinear and soft radiation originating from these
energetic particles. An IRC safe physical observable needs to include all such
collinear radiation, which means that in (1+1) dimensions any IRC safe
observable has to include all states of the theory, and therefore be
completely inclusive of any additional radiation. Therefore, the only IRC safe
observables are those defined by the hard kinematics, which are completely
inclusive over the soft radiation calculated here. Since the Wilson line
operator is unitary, the square of the completely inclusive operator is unity
and therefore trivial. In other words, any non-trivial observable in this
theory is by construction IRC unsafe, and we therefore choose the simplest
such observables, namely the expectation value of the Wilson line operator
between the vacuum states, and the transition from the vacuum to a state with
a single momentum mode.

##  IV Constructing the quantum circuits

###  IV.1 Representation of the Hilbert space and relationship of qubit
representation to field values

To represent the Hilbert space of the massless scalar field theory we
discretize the spatial coordinates on a lattice of size $L$ and $N$ sites per
dimension. The value of the scalar field on each lattice site is digitized and
represented by $n_{\phi}=2^{n_{Q}}$ distinct values. One can write the field
as

$\displaystyle\phi_{i}^{(k)}=-\phi_{\mathrm{max}}+k\,\delta\phi\mbox{\quad
for\quad}k=0,\dotsc,n_{\phi}-1\,,\quad\delta\phi=\frac{2\phi_{\mathrm{max}}}{n_{\phi}-1}\,.$
(S42)

Representing the integer $\ket{k}_{i}$ through its binary representation of
the $n_{Q}$ qubits at lattice site $i$, we can define

$\displaystyle\hat{\phi}_{i}\ket{k}_{i}=\phi_{i}^{(k)}\ket{k}_{i}\,.$ (S43)

The value of $\phi_{\mathrm{max}}$ should be chosen to minimize the minimize
the error due to the digitization and depends on the Hamiltonian implemented
on the lattice. The integer state $\ket{k}_{i}$ is represented as usual
through the bitstring of the $n_{Q}$ qubits at the given lattice site

$\displaystyle\ket{k}_{i}=\ket{q_{0}}_{i}\dotsb\ket{q_{n_{Q}-1}}_{i}\,.$ (S44)

The full Hilbert space is then represented through the tensor product of the
states $\ket{k}$ on each lattice site

$\displaystyle\ket{\Psi}=\ket{k}_{0}\dotsb\ket{k}_{N^{d}}\,.$ (S45)

Our explicit circuit constructions in this paper will only use a single
spatial direction, such that we use

$\displaystyle\ket{\Psi}=\ket{k}_{0}\dotsb\ket{k}_{N}\,.$ (S46)

###  IV.2 The free Hamiltonian

The construction of the Hamiltonian of free massless scalar field theory
follows previous work Macridin:2018gdw ; Macridin:2018oli ; Klco:2018zqz . One
can easily convince oneself that the operator $\hat{\phi}_{i}$ can be written
through its action on the $n_{Q}$ qubits at each lattice site

$\displaystyle\hat{\phi}_{i}=\sum_{j=0}^{n_{Q}-1}2^{j}\hat{\sigma}_{z,i}^{(j)}\,,$
(S47)

where the operator $\hat{\sigma}_{z,i}^{(j)}$ is a single $\sigma_{z}$ Pauli
matrix applied to the $j$th qubit of the $i$th lattice register.

The Hamiltonian is a sum over two pieces that do not commute with one another.
The time derivative of the field is the conjugate field
$\pi_{i}\equiv\dot{\phi}_{i}$, and one can write

$\displaystyle H=H_{\pi}+H_{\phi}$ (S48)

with

$\displaystyle H_{\pi}=\variation{x}\sum_{i}\pi_{i}^{2}\,,\quad
H_{\phi}=\variation{x}\sum_{i=0}^{N-1}\phi_{i}[\nabla^{2}\phi]_{i}\,.$ (S49)

As discussed before, the $\phi_{i}[\nabla^{2}\phi]_{i}$ operator only requires
nearest neighbor interactions on the lattice. The time evolution is then
written in terms of the Suzuki–Trotter formula (we give the first order
expression here)

$\displaystyle\bqty{e^{-iHt}}_{n}=\bqty{e^{iH_{\pi}t/n}\,e^{iH_{\phi}t/n}}^{n}\,.$
(S50)

To construct the exponential of the $H_{\pi}$ operator, we use that the $\phi$
and $\pi$ are related by a Fourier transform on the $\phi$ register. Thus, we
can write

$\displaystyle
e^{iH_{\pi}t}=\mathrm{QFT}^{-1}\,e^{i\variation{x}\,t\phi_{i}^{2}}\,\mathrm{QFT}\,,$
(S51)

where $\mathrm{QFT}$ denotes the (symmetrized) Quantum Fourier Transform,
which was discussed in Klco:2018zqz . We do not repeat its circuit here.

Given this, on needs to find a circuit representation
$\exp[i\theta\phi_{i}\phi_{j}]$ for general $i$ and $j$, from which one can
construct both the circuits for $\exp[iH_{\pi}t]$ and $\exp[iH_{\phi}t]$.

Using Eq. S47, one can write

$\displaystyle\hat{\phi}_{i}\hat{\phi}_{j}$
$\displaystyle=\bqty{\sum_{l=0}^{n_{Q}-1}2^{l}\hat{\sigma}_{z,i}^{(l)}}\bqty{\sum_{k=0}^{n_{Q}-1}2^{l}\hat{\sigma}_{z,j}^{(k)}}=\sum_{l=0}^{n_{Q}-1}\sum_{k=0}^{n_{Q}-1}2^{(l+k)}\sigma_{z,i}^{(l)}\sigma_{z,j}^{(k)}\,.$
(S52)

This allows us to write

$\displaystyle\exp[i\theta\hat{\phi}_{i}\hat{\phi}_{j}]$
$\displaystyle=\prod_{l=0}^{n_{Q}-1}\prod_{k=0}^{n_{Q}-1}\exp\\!\bqty{i\,2^{(l+k)}\theta\,\sigma_{z,i}^{(l)}\sigma_{z,j}^{(k)}}\,.$
(S53)

The action
$\exp\\!\bqty{i\,2^{(k+l)}\theta\,\sigma_{z,i}^{(l)}\sigma_{z,j}^{(k)}}$ is
equal to $\exp[i2^{(k+l)}\theta]$ if qubits $q_{i}^{(l)}$ and $q_{i}^{(k)}$
are equal and $\exp[-i2^{(k+l)}\theta]$ if they are opposite. Thus, it can be
implemented by the circuit and the full operator
$\exp\\!\bqty*{i\theta\hat{\phi}_{i}\hat{\phi}_{j}}$ for $n_{q}$ qubits per
lattice site is therefore implemented by stringing together the
$n_{Q}(n_{Q}-1)/2$ different possible 2-qubit circuits shown above.

###  IV.3 Ground state preparation

The ground state of a massless scalar field theory is given by a multivariate
Gaussian

$\displaystyle\ket{\Psi}=\exp\bqty{-\frac{1}{2}\hat{\phi}_{i}G_{ij}\hat{\phi}_{j}}\ket{k_{0}}\dotsb\ket{k}_{n}$
(S54)

While Qiskit provides a function to generate an aribtrary Gaussian
multivariate distribution, the number of gates in the resulting circuit scale
exponentially with the number of qubits in the system. However, an algorithm
with polynomial scaling was derived by Kitaev and Webb kitaev2009wavefunction
. While this algorithm does not produce the exact multivariate distribution in
its digitized form, it approaches the correct limit as the number of qubits
per lattice site becomes large.

The KW algorithm relies on the LDL or square-root-free Cholesky decomposition,
which rewrites the correlation matrix in terms of a diagonal matrix $D$ and a
lower unit-triangular matrix $L$ (a lower triangular matrix with 1’s on the
diagonal)

$\displaystyle G=LDL^{\dagger}$ (S55)

An arbitrary multi-dimensional Gaussian distribution can then be created by
generating a series of uncorrelated Gaussian distributions according to the
diagonal matrix $D$, and then applying shear matrices $L$ through a remapping
of the basis states of the Hilbert space. For details of this algorithm, we
refer the reader to the original paper. In Section II we mentioned a modified
version of the KW procedure, which only rounds the shear matrices applied to
the uncorrelated states to the nearest digitized field value after applying
the the full shear matrix rather than after every individual shearing
operation, as was originally proposed in kitaev2009wavefunction . While
requiring more ancilla qubits for memory, this does not affect the polynomial
scaling of the approximation and results in exponentially greater fidelity
with the exact ground state for the $n_{Q}\geq 3$ cases futureKW .

For 2 qubits per lattice site, the second step of the KW algorithm does not
actually change the basis states, such that KW state preparation is equal to
the production of a uncorrelated Gaussians at each lattice site with width
given by the diagonal entries of the matrix $D$. In this work, we use the
Cholesky decomposition of the correlation matrix, but rather than using the
full KW algorithm to produce the uncorrelated Gaussian, we use the built in
functionality of Qiskit, even though it has exponential scaling with $n_{Q}$.

###  IV.4 Wilson line operator

The Wilson line operator on a lattice was given in Eq. S22. From this
expression on can see that it is determined through the successive application
of the operator $\exp[ig\variation{x}\phi_{n}]$ and Hamiltonian evolution.
Given the circuit for Hamiltonian evolution derived above, one therefore only
needs the circuit for the exponential of a single field. This is easily
derived using a similar derivation to the Hamiltonian case, and one writes

$\displaystyle\exp\\!\bqty*{i\theta\hat{\phi}_{i}}=\prod_{j=0}^{n_{Q}-1}\exp\\!\bqty{i2^{j}\theta\sigma_{z,i}^{(j)}}\,,$
(S56)

which can easily be written in circuit form

##  V Validation of the circuits for the Hamiltonian

As a validation of the quantum circuits we first check the implementation of
the free field theory (ground state preparation and time evolution). In
particular, we check to what degree the ground state is an eigenstate of the
Hamiltonian, and how well the energy of the ground state agrees with the known
analytical value. We begin by computing the overlap

$\displaystyle
f(t)=\absolutevalue{\\!\matrixelement{{\Omega}}{\bqty{e^{-iHt}}_{n}}{{\Omega}}}^{2}\,,$
(S57)

where $\bqty{e^{-iHt}}_{n}$ is the Trotterized Hamiltonian with $n$ steps
given in Eq. S50. This is implemented with the circuit

and the function $f(t)$ is obtained by the fraction of measurements where all
qubits are back the initial $\ket{0}$ state. If the ground state is indeed an
eigenstate of the Hamiltonian, and the Trotterized Hamiltonian is equal to the
full Hamiltonian, this function should be identically 1, and any deviation
from this value should be due to Trotterization errors or because the state
$\ket{{\Omega}}$ not being equal to the true ground state. In Fig. 5 we show
this overlap for the exact ground state on the left and for the KW ground
state on the right. The result on the left confirms that with more Trotter
steps we approach unity, while the result on the right shows that even for a
large amount of Trotter steps the KW approximation leads to deviations from
unity.

Figure 4: The value of
$\absolutevalue{\\!\matrixelement{{\Omega}}{\bqty{e^{-iHt}}_{n}}{{\Omega}}}^{2}$
for two different values of $n$. The blue line shows the result for $n=1$,
while the red line shows the result for $n=\lceil t/0.2\rceil$. Measurements
from a noiseless simulation of 512 shots/point are overlayed. On the left, we
show the result for the exact ground state, while on the right we show the
result for the KW approximation. For the exact ground state, the deviation
from unity is only due to the Trotterization of the Hamiltonian, and gets
smaller with more Trotter steps (and more quantum gates). Since the KW
approximation is not an eigenstate of the full Hamiltonian, the deviation from
unity is due to both this approximation and the Trotterization.

While this measurements tests to what degree the ground state is an eigenstate
of the Trotterized Hamiltonian, it can not check for the energy of the ground
state since that manifests itself only as a pure phase in the above circuit. A
slight variation of this circuit

can measure the energy of the ground state. The fraction of measurements with
all qubits in the $\ket{0}$ state is given by

$\displaystyle
f_{\mathrm{ctr}}(t)=\frac{1}{4}\absolutevalue{1+\matrixelement{{\Omega}}{\bqty{e^{-iHt}}_{n}}{{\Omega}}}^{2}\,,$
(S58)

which is sensitive to the energy $E_{\Omega}$. One can see that using the
digitization of exact ground state and sufficient Trotter steps one produces
the analytically expected dependence on the ground state energy up to the
small difference in period due to the shift in ground state energy due to
digitization. Conversely, with the KW states one also sees a small reduction
in probability due to leakage out of the approximate ground state, as
expected.

Figure 5: The value of
$\absolutevalue{1+\matrixelement{{\Omega}}{\bqty{e^{-iHt}}_{n}}{{\Omega}}}^{2}/4$
for two different values of $n$. The blue line shows the result for $n=1$,
while the red line shows the result for $n=\lceil t/0.2\rceil$. Measurements
from a noiseless simulation of 512 shots/point are overlayed. On the left, we
show the result for the exact ground state, while on the right we show the
result for the Kitaev–Webb approximation. In dashed black we show the exact
result.

##  VI Error Mitigation on IBMQ

We mitigate both readout errors and gate errors. Readout error mitigation
proceeds with a classical post-processing step. We prepare all
$2^{n_{\text{qubit}}}$ possible states $\ket{i}$ and measure the frequency of
observing the state $\ket{f}$. These probabilities are tabulated in a
$2^{n_{\text{qubit}}}\times 2^{n_{\text{qubit}}}$ matrix called the response
matrix $R$. The observations from the main experiment are represented as a
vector $m_{i}$ and then readout error mitigation proceeds iteratively
ReadoutCorrection ; DAgostini:1994fjx ; 1974AJ…..79..745L ; Richardson:72 :

$\displaystyle t_{i}^{n+1}$ $\displaystyle=\sum_{j}\Pr(\text{truth is
$i$}|\text{measure $j$})\times m_{j}$
$\displaystyle=\sum_{j}\frac{R_{ji}t_{i}^{n}}{\sum_{k}R_{jk}t_{k}^{n}}\times
m_{j}\,,$ (S59)

where $t_{i}^{0}=1/2^{n_{\text{qubit}}}$ is the uniform prior and the number
of iterations $n$ is chosen to be 20. The result is not sensitive to small
variations in these choices. The iterative procedure in Eq. VI avoids
pathologies from other forms of regularized matrix inversion and can be
further improved with variations like readout rebalancing 2010.07496  and
subexponential approximations geller_efficient_2020 ; song_10-qubit_2017 ;
gong_genuine_2019 ; wei_verifying_2020 ; hamilton2020scalable . Matrix
inversion approaches from the quantum simulators pyQuil pyquil  (by Rigetti),
Cirq cirq ; arute2020quantum  (by Google), and XACC alex2019xacc ; xacc ;
mccaskey_quantum_2019  and the least squares method from Qiskit by IBM Qiskit
; ignis  (see also Ref. 1904.11935 ; 1907.08518 ) produce similar results
ReadoutCorrection . In current IBMQ machines, the dominant gate noise is from
multiqubit operations, which are only the CNOT gates in our experiments. A
common model for CNOT errors is the quantum depolarizing noise channel. This
error is mitigated by systematically increasing the noise and then
extrapolating to zero noise. In particular, we remove the leading order
depolarizing noise by repeating the nominal measurement with a circuit where
every CNOT gate is replaced with 3 CNOT gates. Following Ref. He:2020udd , we
then set the mitigated result as $3/2$ times the nominal circuit minus $1/2$
the result from the circuit result with triple the CNOT gates. This fixed
identity insertion method (FIIM) was first proposed in Ref. Dumitrescu:2018
and multiple variations exist that improve the fidelity of the extrapolation
or are more resource efficient PhysRevX.8.031027 ; PhysRevLett.119.180509 ;
He:2020udd .