Simulating collider physics on quantum computers using effective field theories

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 [1]. The main technique involves disretizing the spatial degrees of freedom by introducing a lattice [2][3][4], and digitizing the field values at a given lattice point [5][6][7][8][9][10][11][12][13]. This turns the uncountably infinite dimensional Hilbert space of standard quantum field theories into a finite dimensional Hilbert space of dimension where n φ 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 δx in each direction and is given by The total number of qubits required for such a simulation is given by 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 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 O(10 MeV) (the smallest resolvable transverse momentum between hadrons) and 7 TeV (the beam energy of the LHC). To fully capture this energy range would require a lattice with 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 [14,15] 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 [16][17][18] and parton showers [19,20] 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 [21], it would be a significant breakthrough were it possible to simulate these low arXiv:2102.05044v1 [hep-ph] 9 Feb 2021 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 [22][23][24][25][26][27]). 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 O(1 TeV) and an invariant mass of O(50 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 GeV E 2 TeV. As we discuss below, the EFT only needs to capture 1 GeV E 50 GeV. The number of qubits required to simulate the EFT is smaller by a factor of (2000/50) 3 ≈ 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) [28][29][30][31]. Using this EFT, a typical cross section describing a physical observable at the LHC can be factorized into separate pieces [32][33][34][35] 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 [36,37]. 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 [31] Y n = P exp ig ∞ 0 ds n·A(x µ = n µ s) .
Here A µ is the soft gauge field and n µ describes the direction of the light-like direction n µ = (1, n) with n 2 = 1 such that n µ n µ = 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 µ . 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 † n for the antiparticle moving in thē n µ = (1, −n) direction. The matrix elements required in the soft sector are given by where T denotes the time-ordering operator, |Ω denotes the ground state of a Hilbert space containing only the gauge degree of freedom, and |X is the final state. If reproducing multi-particle weakly coupled final states, |X will contain a fixed number of gauge momentum modes with given momenta k µ 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 [38][39][40][41][42][43][44] can only compute these matrix elements for final states |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 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 , . . . , x N −1 . To eliminate the zero-momentum mode of the theory, we impose twisted boundary conditions [45][46][47][48]. The result is a theory defined at discrete positions x and momenta p given by x i = x min + i δx and p i = p min + i δp with x min = −(N −1) δx /2, p min = −π/ δx and δp = 2π/N δx, Writing φ i ≡ φ(x i ), the twisted boundary conditions correspond to the condition φ i+N = −φ i . The Hamiltonian becomes [5] where the lattice operator ∇ 2 is defined through its action on a field as [ The Wilson line operators can be written as 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 [8]. Choosing n Q qubits per lattice site allows for n φ ≡ 2 n Q different field values. For each lattice point, the possible field values are chosen to be by φ The value of φ max has to be chosen to optimize the digitized description, which for free fields is accomplished by Forω = 1, as is the case for a single lattice site with ω = 1, corresponding to a single harmonic oscillator, Eq. (11) reproduces the empirical numerical values obtained in [8].
To implement the Wilson line operator we first rewrite the time-ordered product of the two Wilson lines as 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 max = n 0 δ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 → ∞, δx → 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 [5]. 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 nonlocal 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 [8] and uses the fact that the digitized field φ (k) i can be written in terms of sums of σ z operators. This implies that the exponential of products of fields φ i can be implemented through combinations of CNOT gates and R Z rotations [6][7][8]. The exponential of the conjugate operatorφ 2 can be implemented by taking a quantum Fourier transformation of the exponential of the operator φ 2 . These can then be combined via the Suzuki-Trotter formula [49][50][51]. 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) [52]. To identify states of definite multiplicity and momentum in |X one can follow the general ideas laid out in [1,5].
Our quantum circuit has been implemented in Qiskit [53] and is available from the authors upon request. In this first exploratory paper we compute the foundational quantities, namely for |X = |Ω and |X = |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 where |l n denotes the register of qubits for the n th lattice site. This creates the multivariate Gaussian vacuum state from the initial state with all qubits zero using U Ω , 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 |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 δx = 1. With only 3 lattice sites, the Wilson line operator simplifies to 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 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 [5] will be necessary; the effect of such approximations is presented in the Supplementary Material. 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 [54] as well as CNOT gate noise mitigation [55]. For more details, including References , 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 %. This would be reduced to at most 1.5 % 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 % 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 % of the analytical result, and by using noise-mitigation techniques, uncertainties due to running on present-day hardware can be reduced to about 5 %.

Simulating effective field theories on quantum computers Supplementary Material
Christian W. Bauer, Marat Freytsis, Benjamin Nachman

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 [5], while the calculations involving the Wilson lines have not appeared anywhere else to our knowledge.

Definition of the lattice
We consider a lattice of d dimensions, N lattice sites for each dimension and lattice spacing δx. Points on the latice are indexed by d integers ranging from 0 to N − 1. We denote this set of integers by with the the positions of the lattice in each dimension given by The corresponding momentum values in each dimension are then given by (s denotes the same set of integers as r) 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 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 x r = r δx and p s = (s + ∆) δ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 Note that we have introduced a short-hand form for the sum over lattice sites as where r and s denote the positions on the regular and conjugate lattice, respectively. One has the usual completeness relation were we have again used a short hand notation and we have chosen the factor of 2π 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 ∆. 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 (S10) Thus, non-integer values of ∆ i give an extra phase in the boundary condition, with half-integer values adding a relative minus sign, giving twisted boundary conditions.
The free field theory The free scalar field theory on the lattice is given by the Hamiltonian 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 where frequencies are given by (S13) Introducing the creation and annihilation operators aṡ which satisfy the commutation relations one finds for the Hamiltonian (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 × n occupation numbers l (0) · · · l (n−1) , where the occupation number at each lattice site is determined by d integers ≥ 0. Thus, a general state is given by with the ground state given by The energy of the ground state |Ω is given by while the energy of excited states are easily determined in terms of their occupation numbers as Wilson Line expectation values Using this notation, one can compute the expectation values of the Wilson line operator T [Y n Y † n ]. 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 The Wilson line operator is therefore given by where we have used thatn = −n and x n0+m = −x n0−m . In the last line we defined as well as t ≡ n 0 δx, δt ≡ δx. Note that since x n0 = 0 one has U 0 (t) = e −iHt . Using Eq. (S14) we can write 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 which satisfies the relation and displacement operators acting on different momentum modes commute The action of the time evolution on the displacement operator is Combining all this information, one finds To compute expectation values of this operator, we use In particular, the magnitude square of the vacuum expectation value is given by where the sum over x and y runs over the lattice position and we have defined One can verify that this expression reproduces the perturbative result to order g 2 . To that order, one can write Taking the square of this result one finds which is the continuum limit of Eq. (S33).

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 with Hamiltonian The momenta of the 3-site lattice are given by The corresponding frequencies are obtained from Eq. (12) 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/ δx. The excited states are simply given by the ground state energy increased by up to n φ additions of each frequency. Thus, the lowest lying states are 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 %. This is in agreement with the findings of [8].
In Fig. 2 we show the dependence of the vacuum expectation value 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 [52] and used in the rest of the text with a modified form discussed in Section . The blue solid lines are the analytical result for transitions to the vacuum (Ω) and the lowest-lying one one-particle state (p1). Digitization limits the number of states the field can take on for nQ = 1 (red) to 2,nQ = 2 (orange) to 4 and nQ = 3 green to 8. In each case, the exact ground state is digitized and evolved.

CONSEQUENCES OF THE NON-LOCAL NATURE OF WILSON LINES SPECIAL CONSIDERATIONS IN (1+1) DIMENSIONS
The two Wilson lines Y n and Y † n 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 ∞ along the directions n µ = (1, n) andn µ = (1, −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 δx. operators. The presence of the finite lattice produces terms L n · p and Ln · p, regulating divergences as n · p andn · p tend to zero. On the other hand, for finite lattice spacing one has ω − |p| ∼ |p| 2 δ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 · p → 0 with p → ∞). 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 orn, which means that every momentum satisfies n · p = 0 or n · 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 µ andn µ . 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.

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 φ = 2 n Q distinct values. One can write the field as Representing the integer |k i through its binary representation of the n Q qubits at lattice site i, we can definê The value of φ 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 |k i is represented as usual through the bitstring of the n Q qubits at the given lattice site The full Hilbert space is then represented through the tensor product of the states |k on each lattice site Our explicit circuit constructions in this paper will only use a single spatial direction, such that we use The free Hamiltonian The construction of the Hamiltonian of free massless scalar field theory follows previous work [6][7][8]. One can easily convince oneself that the operatorφ i can be written through its action on the n Q qubits at each lattice sitê where the operatorσ z,i is a single σ z Pauli matrix applied to the jth qubit of the ith 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 π i ≡φ i , and one can write As discussed before, the φ i [∇ 2 φ] 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) e −iHt n = e iHπt/n e iH φ t/n n .
To construct the exponential of the H π operator, we use that the φ and π are related by a Fourier transform on the φ register. Thus, we can write where QFT denotes the (symmetrized) Quantum Fourier Transform, which was discussed in [8]. We do not repeat its circuit here. Given this, on needs to find a circuit representation exp[iθφ i φ j ] for general i and j, from which one can construct both the circuits for exp[iH π t] and exp[iH φ t].
Using Eq. (S47), one can writê This allows us to write The action exp i 2 (k+l) θ σ i and q (k) i are equal and exp −i2 (k+l) θ if they are opposite. Thus, it can be implemented by the circuit and the full operator exp[iθφ iφ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.

Ground state preparation
The ground state of a massless scalar field theory is given by a multivariate Gaussian 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 [52]. 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) 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 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 [52]. 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 ≥ 3 cases [56]. 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 .

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 δx φ 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 which can easily be written in circuit form

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 where 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 |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 |Ω 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. 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.
which is sensitive to the energy E Ω . 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.

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 qubit possible states |i and measure the frequency of observing the state |f . These probabilities are tabulated in a 2 n qubit × 2 n 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 [54,[57][58][59]: where t 0 i = 1/2 n 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. S59 avoids pathologies from other forms of regularized matrix inversion and can be further improved with variations like readout rebalancing [60] and subexponential approximations [61][62][63][64][65]. Matrix inversion approaches from the quantum simulators pyQuil [66] (by Rigetti), Cirq [67,68] (by Google), and XACC [69][70][71] and the least squares method from Qiskit by IBM [53,72] (see also Ref. [73,74]) produce similar results [54]. In current IBMQ machines, the dominant gate noise is from multiqubit