# Simulating collider physics on quantum computers using effective field theories

SSimulating collider physics on quantum computers using eﬀective ﬁeld theories

Christian W. Bauer ∗ and Benjamin Nachman † Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

Marat Freytsis ‡ NHETC, Department of Physics and Astronomy,Rutgers University, Piscataway, NJ 08854, USA andPhysics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

Simulating the full dynamics of a quantum ﬁeld theory over a wide range of energies requiresexceptionally large quantum computing resources. Yet for many observables in particle physics,perturbative techniques are suﬃcient to accurately model all but a constrained range of energieswithin the validity of the theory. We demonstrate that eﬀective ﬁeld theories (EFTs) provide aneﬃcient mechanism to separate the high energy dynamics that is easily calculated by traditionalperturbation theory from the dynamics at low energy and show how quantum algorithms can beused to simulate the dynamics of the low energy EFT from ﬁrst principles. As an explicit example wecalculate the expectation values of vacuum-to-vacuum and vacuum-to-one-particle transitions in thepresence of a time-ordered product of two Wilson lines in scalar ﬁeld theory, an object closely relatedto those arising in EFTs of the Standard Model of particle physics. Calculations are performed usingsimulations of a quantum computer as well as measurements using the IBMQ Manhattan machine.

It is well known that quantum computers can in prin-ciple simulate the time evolution of quantum ﬁeld theo-ries [1]. The main technique involves disretizing the spa-tial degrees of freedom by introducing a lattice [2–4], anddigitizing the ﬁeld values at a given lattice point [5–13].This turns the uncountably inﬁnite dimensional Hilbertspace of standard quantum ﬁeld theories into a ﬁnite di-mensional Hilbert space of dimension n H = n ( N d ) φ , (1)where n φ denotes the dimensionality of the Hilbert spacefor a given lattice point, N is the total number of lat-tice points in each spatial direction, and d represents thenumber of spatial dimensions. The physical volume ofthe lattice is determined by the distance between adja-cent lattice points δx in each direction and is given by V = ( N δx ) d ≡ L d . (2)The total number of qubits required for such a simulationis given by n Q = O (log n H ) = O ( N d log n φ ) . (3)The discretization and ﬁnite volume of space introduceupper and lower cutoﬀs to the energies E over which theresulting lattice ﬁeld theory is a good approximation forthe continuum. In particular, one ﬁnds1 N δx (cid:46) E (cid:46) δx , (4)which implies that the range of energies that can bedescribed is directly proportional to the number of lat-tice sites per dimension. In principle, to have access tothe full dynamics of the Large Hadron Collinder (LHC), one would need to describe the energy range between O (10 MeV) (the smallest resolvable transverse momen-tum between hadrons) and 7 TeV (the beam energy of theLHC). To fully capture this energy range would require alattice with O (10 ) lattice points in each dimension, andmore than 10 qubits to reproduce the resulting physi-cal system. Even if one does not require the full energyrange up to the LHC center-of-mass energy, the numberof qubits required for a full simulation will clearly remainbeyond the realm of feasibility for a long time to come.For many observables of interest at the LHC, physicsat short distances is reliably computed in ﬁxed order per-turbation theory, and high precision can be reached withexisting techniques. Physics at lower energies introducesnew signiﬁcant challenges. Asymptotic freedom [14, 15]implies that the strong coupling constant becomes largeat low energies, increasing the production of additionalparticles. An energetic particle can thus easily radiateadditional soft and collinear particles, resulting in a col-limated jet of particles. This means that ﬁxed order per-turbative calculations are of little use to describe lowerenergy eﬀects such as the precise makeup of jets, andother techniques such as resummation [16–18] and par-ton showers [19, 20] have to be used to make predic-tions. At even lower energies, the fully nonperturba-tive dynamics of hadronization dominate. While clas-sical lattice methods provide spectral information aboutthe strongly-coupled limit of the theory, the dynamics ofhadronization is currently only understood through phe-nomenological models and form factors extracted fromdata. While these techniques have been quite success-ful and it has even been shown that quantum algorithmscan be used to include quantum interference eﬀects insome models of parton showers [21], it would be a signif-icant breakthrough were it possible to simulate these low a r X i v : . [ h e p - ph ] F e b energy dynamics from ﬁrst principles.In this paper we show that this is indeed possible byusing eﬀective ﬁeld theories (EFTs) which have been de-signed to reproduce the desired long distance physics.(For work on simulating EFTs not related to colliderphysics, see [22–27]). The dynamics of these EFTs canthen be simulated directly and from ﬁrst principles on aquantum computer. Since the energy range that needsto be simulated is much smaller than that of the full the-ory, the resource requirements are orders of magnitudessmaller than those required to simulate the full theory.Additionally, the EFT simpliﬁes the description of theinitial and ﬁnal state and signiﬁcantly reduces the re-source requirements of their implementation. For exam-ple, consider an observable measured on two jets, eachwith an energy of O (1 TeV) and an invariant mass of O (50 GeV). Even restricting ourselves to observables in-sensitive to hadronization, a generic observable of themomenta of ﬁnal state particles in the full theory re-quires a simulation of the range 1 GeV (cid:46) E (cid:46) (cid:46) E (cid:46)

50 GeV. The number of qubits re-quired to simulate the EFT is smaller by a factor of(2000 / ≈ , and this factor increases rapidly asthe range of the full theory is increased.The EFT relevant for hadronic jet physics is Soft-Collinear Eﬀective Theory (SCET) [28–31]. Using thisEFT, a typical cross section describing a physical ob-servable at the LHC can be factorized into separatepieces [32–35] σ = H ⊗ J ⊗ · · · ⊗ J n ⊗ S . (5)Here H denotes a coeﬃcient describing the short dis-tance physics which can be computed reliably in standardperturbation theory, while J i and S denote squares ofmatrix elements of operators in SCET. The jet function J i denotes physics arising from collinear degrees of free-dom that all have small momentum relative to each other(while moving collectively with a large energy), with thediﬀerent jet functions completely decoupled from one an-other. The soft function S denotes physics arising fromsoft degrees of freedom, which all have small absolutemomentum (in the frame of the collider). These ma-trix elements describe the transition of a simple initialstate produced in the short distance interaction H intoa ﬁnal state, with the dynamics described by the EFTHamiltonian. The ﬁnal states that arise can contain alarge number of particles. State of the art techniques useoperator renormalization to compute the overall scale de-pendence of the jet and soft functions, and then computethe relevant matrix elements perturbatively, such thatthe eﬀects arising from the high multiplicity ﬁnal statesor hadronization can not be properly included. Classicallattice techniques are also not suitable for the compu-tation of matrix elements of SCET operators, since thelong-distance dynamics is governed by massless modes which are inherently Minkowskian in nature. Having asimulation of the dynamics of SCET (or a diﬀerent EFTfor other problems) will allow for a full non-perturbativecalculation of any matrix elements.The EFT reproduces the full theory result up to powercorrections, whose size depends on the kinematics of theprocess studied. For many observables of interest, thesepower corrections are considerably smaller than the per-turbative corrections present in either the full theory orEFT matrix elements.Since individual jet functions do not interact with oneanother, their dynamics can be simulated in the refer-ence frame where all particles have small absolute mo-mentum [36, 37]. This requires simulating the dynamicsof the original ﬁeld theory (with any degrees of freedomthat only contribute to short distance eﬀects removed),but with a much smaller energy range required to be sim-ulated. To achieve this on a Quantum Computer, oneneeds a setup very similar to that for the full theory, butreliable calculations can be achieved with much coarserlattices than those one would need for a simulation of thefull theory.In the soft function, on the other hand, the energeticparticles no longer contribute to the dynamics of the the-ory. Instead, their eﬀect is captured by a so-called Wilsonline, which describes the interaction of a charge movingalong a ﬁxed world-line with the bath of soft degrees offreedom. The physics underlying this observation is thatsoft particles cannot change the direction of energeticparticles in a meaningful way. This implies that ener-getic particles can be included in the EFT as a staticobject (Wilson line) described by the relevant quantumnumbers (charge, color, etc.) moving along a ﬁxed worldline along a light-like direction. Matrix elements in thesoft theory therefore need to compute the dynamics of aHamiltonian describing the soft bath of particles in thepresence of such Wilson lines. For gauge theories, such asthose describing the three fundamental forces containedin the Standard Model, Wilson lines are given by pathordered exponentials of the relevant gauge ﬁelds [31] Y n = P exp (cid:20) ig (cid:90) ∞ d s n · A ( x µ = n µ s ) (cid:21) . (6)Here A µ is the soft gauge ﬁeld and n µ describes the di-rection of the light-like direction n µ = (1 , n ) with n = 1such that n µ n µ = 0. Thus, the gauge ﬁeld is evaluatedon a path going from the origin to spatial inﬁnity alongthe world-line characterized by the direction n µ . Thedynamics of the soft gauge ﬁeld is described by its fulltheory Hamiltonian.The soft function for a process containing two ener-getic particles of zero total charge moving back-to-backrequires two Wilson lines, Y n for the particle moving inthe n direction and Y † ¯ n for the antiparticle moving in the¯ n µ = (1 , − n ) direction. The matrix elements required inthe soft sector are given by (cid:104) X | T[ Y n Y † ¯ n ] | Ω (cid:105) , (7)where T denotes the time-ordering operator, | Ω (cid:105) denotesthe ground state of a Hilbert space containing only thegauge degree of freedom, and | X (cid:105) is the ﬁnal state. If re-producing multi-particle weakly coupled ﬁnal states, | X (cid:105) will contain a ﬁxed number of gauge momentum modeswith given momenta k µi . To compute the soft function fora given observable one needs to sum over all possible ﬁnalstates that can contribute to this observable. Traditionalperturbative calculations [38–44] can only compute thesematrix elements for ﬁnal states | X (cid:105) containing a smallnumber of particles and at low order in perturbation the-ory, since the complexity of the calculation increases fac-torially with the power of the coupling constant g . Theytherefore do not compute the full matrix element for agiven observable, but only a perturbative approximation.For large coupling constants, which is the relevant casefor Quantum Chromodynamics (QCD) at low energies,such perturbative calculations can give rise to large un-certainties.Simulating the full dynamics of the ﬁeld theory wouldprovide the full non-perturbative result for the soft ma-trix element. To evaluate the matrix element on a quan-tum computer one ﬁrst needs to deﬁne circuits that per-form time evolution of the system, as well as circuits thatcan create and measure the ground and exited states ofthe theory. In addition, one needs to create circuits thatcan correctly interleave the implementation of nonlocalWilson line operators with the evolution of the system toreproduce their time-ordered product.In the following we discuss how to compute matrix el-ements of Wilson line operators Y n analogous to Eq. (7),but for a massless scalar, rather than gauge, theory. TheEFT is insensitive to the precise origin of the Wilsonlines, but a particualy straightforward realization wouldresult from a pair of highly-energetic fermions coupledto massless scalars through a Yukawa coupling. Whenconstructing the explicit circuit we also limit ourselvesto (1+1) dimensions, mainly to restrict the quantum re-sources required such that it can be implemented on cur-rently existing hardware. This allows us to omit sometechnical complications that arise when dealing withgauge theories (gauge transformations, the existence ofunphysical polarizations, etc .), while capturing all theresulting simpliﬁcation of working within an EFT.To be precise, we consider a massless ﬁeld theory in(1+1) dimensions with Hamiltonian and Wilson lines de-ﬁned by H = (cid:90) d x (cid:16) ˙ φ − φ ∂ φ (cid:17) ,Y n = P exp (cid:20) ig (cid:90) ∞ d s φ ( x µ = n µ s ) (cid:21) . (8) In the Supplemental Material, we discuss how working in(1+1) dimensions gives rise to several eﬀects not presentin higher dimensions, and how these generalize to higherdimensions.We discretize the position x into an odd number oflattice points, labeling the positions by x , . . . , x N − . Toeliminate the zero-momentum mode of the theory, weimpose twisted boundary conditions [45–48]. The resultis a theory deﬁned 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 − δx / p min = − π/ δx and δp = 2 π/N δx ,Writing φ i ≡ φ ( x i ), the twisted boundary conditions cor-respond to the condition φ i + N = − φ i . The Hamiltonianbecomes [5] H = δx N − (cid:88) i =0 (cid:104) ˙ φ i − φ i [ ∇ φ ] j (cid:105) , (9)where the lattice operator ∇ is deﬁned through its ac-tion on a ﬁeld as [ ∇ φ ] i = (2 φ i − φ i − − φ i +1 ) / δx . Dueto the twisted boundary conditions [ ∇ φ ] = (2 φ + φ N − − φ ) / δx and [ ∇ φ ] N − = (2 φ N − − φ N − + φ ) / δx . The Wilson line operators can be written as Y n = P exp (cid:34) ig δx n (cid:88) i = n φ x i ( t = x i − n ) (cid:35) ,Y † ¯ n = P exp (cid:34) − ig δx n (cid:88) i =0 φ x i ( t = n − x i ) (cid:35) , (10)where n = ( N − / n Q qubits per lattice site allows for n φ ≡ n Q diﬀerent ﬁeld values. For each lattice point, the possibleﬁeld values are chosen to be by φ ( k ) i = − φ max + k δφ , with δφ = 2 φ max / ( n φ − φ max has to be chosento optimize the digitized description, which for free ﬁeldsis accomplished by φ max = 1 √ δx ¯ ω (cid:115) π n φ − n φ , (11)where ¯ ω = 1 N (cid:88) i ω i , ω i = 2 δx (cid:12)(cid:12)(cid:12)(cid:12) sin p i δx (cid:12)(cid:12)(cid:12)(cid:12) . (12)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 ﬁrst rewritethe time-ordered product of the two Wilson lines asT[ Y n Y † ¯ n ] = e − iH n δx exp (cid:2) ig δx (cid:0) φ x n − φ x (cid:1)(cid:3) (13) × e iHδx exp (cid:2) ig δx (cid:0) φ x n − − φ x (cid:1)(cid:3) × · · · × e iHδx exp (cid:2) ig δx (cid:0) φ x n − φ x n (cid:1)(cid:3) , where we have used the time translation operator tomake the time dependence on the ﬁeld operators explicit.Thus, the Wilson line operator consists of a sequence oftime-evolution operators for a time interval correspond-ing to the lattice spacing and exponentials of the ﬁeldoperator. The last time evolution evolves the state backfrom the largest time to which the Wilson lines can besensibly evolved, namely t max = n δx , to t = 0 at whichall states are deﬁned.Ultimately, to make contact with the continuum ﬁeldtheory any such simulation will have to be performed ona series of increasing lattices, and the result extrapolatedto the N → ∞ , δx → φ ( k ) i can be written in terms of sums of σ z opera-tors. This implies that the exponential of products ofﬁelds φ i can be implemented through combinations ofCNOT gates and R Z rotations [6–8]. The exponential ofthe conjugate operator ˙ φ can be implemented by tak-ing a quantum Fourier transformation of the exponentialof the operator φ . These can then be combined viathe Suzuki–Trotter formula [49–51]. For details, see theSupplemental Materials. The initial ground state of thescalar ﬁeld theory is a multi-variate Gaussian distribu-tion, which can be created using the approach of Kitaevand Webb (KW) [52]. To identify states of deﬁnite multi-plicity and momentum in | X (cid:105) one can follow the generalideas laid out in [1, 5].Our quantum circuit has been implemented in Qiskit [53] and is available from the authors upon re-quest. In this ﬁrst exploratory paper we compute thefoundational quantities, namely Y X = (cid:12)(cid:12)(cid:12) (cid:104) X | T[ Y n Y † ¯ n ] | Ω (cid:105) (cid:12)(cid:12)(cid:12) (14) for | X (cid:105) = | Ω (cid:105) and | X (cid:105) = | p i (cid:105) , the one-particle momentumeigenstates of the theory. It should be noted that thesequantities are not infrared (IR) safe, and will thereforedepend on the IR scale in the problem, the lattice size L .However, as discussed in more detail in the SupplementalMaterials, there is no non-trivial IR safe observable thatcan be deﬁned in (1+1) dimensions, and these transitionrates are therefore representative quantities of what canbe computed in this theory.The quantum circuit for this measurement can be rep-resented as | l (cid:105) / U Ω U Y U † X | . . . (cid:105) / | l N − (cid:105) / where | l n (cid:105) denotes the register of qubits for the n th lat-tice site. This creates the multivariate Gaussian vacuumstate from the initial state with all qubits zero using U Ω ,acts on this vacuum with the time ordered product ofthe two Wilson lines using U Y , and ﬁnally applies the in-verse of the state preparation of state | X (cid:105) . The details ofthese various circuits can be found in the SupplementaryMaterial.For our numerical results, we work with an N = 3 sitelattice with spacing δx = 1. With only 3 lattice sites, theWilson line operator simpliﬁes to Y X = (cid:12)(cid:12)(cid:12) (cid:104) X | T[ Y n Y † ¯ n ] | Ω (cid:105) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) (cid:104) X | e igδx ( φ x − φ x ) | Ω (cid:105) (cid:12)(cid:12)(cid:12) , (15)since all time evolution operators act on the initial or ﬁ-nal eigenstate of the Hamiltonian only and can thereforebe neglected as contributing an overall phase. In Fig. 1we show the dependence of the expectation values Y X on the coupling g for n Q = 2 qubits per lattice site fordiﬀerent ﬁnal states, and compare them against the ana-lytical results, shown by black lines. Results are given forboth a quantum simulation and from the 65-qubit IMBQManhattan quantum computer. The operators for imple-menting all states are exact, as the resources for doing soon a small lattice are modest. On a larger lattice approx-imate methods, such as KW ground state approximationand the excited state preparation techniques of [5] will benecessary; the eﬀect of such approximations is presentedin the Supplementary Material.Errors in the quantum circuits, especially readout er-rors and CNOT gate errors are quite large on existinghardware. As discussed in the Supplemental Materials,the exponential of the ﬁeld operator at a given positionrequires only n Q single qubit gates, such that the opera-tor in Eq. (15) requires no CNOT gates. For n Q = 2, hestate preparation requires 6 CNOT gates for gates. Notethat for more than 3 lattice sites the time evolution op-erator is required, which requires a much larger amount g T r a n s i t i o n r a t e | X | T [ Y n Y n ] || X = | X = | p Analytic CalculationDigitized Calculation/Noiseless Simulation (Qiskit)Raw Quantum (IBMQ)Corrected Quantum (IBMQ)

FIG. 1. Result of transition rates from the valuum of theWilson line for 3 lattice sites and n Q = 2 qubits per siteto the vacuum and the lowest-lying single excited state. Thesolid lines shows the analytical result with no ﬁeld digitizationwhile the dashed lines represents the result from a quantumsimulator of our circuit. The black data points show the re-sult from the 65-qubit IBMQ Manhattan quantum computer,corrected both for readout errors and CNOT gate errors. Weonly show result from the Manhattan computer for X = Ω,since the circuit to measure the excited state was too deep togive reliable results. of gates, although the resulting circuits are known. Forexample, even for 3 lattice sites the standard implemen-tation of a single Trotter step of our Hamiltonian re-quires 60 CNOT gates. We have applied both readouterror mitigation as described in [54] as well as CNOTgate noise mitigation [55]. For more details, includingReferences [56–77], see the Supplemental Materials. Onecan see that the digitized result with 2 qubits per latticesite diﬀers from the analytic calculation by up to 10 %.This would be reduced to at most 1 . ACKNOWLEDGEMENTS

We would like to thank Dorota Grabowska, Bert deJong, Michael Kreshchuk, Pier Monni, John Preskill,Martin Savage and Miro Urbanek for useful discus-sions. CWB and BN are supported by the U.S. Depart-ment of Energy (DOE), Oﬃce of Science under contractDE-AC02-05CH11231. In particular, support comesfrom Quantum Information Science Enabled Discovery(QuantISED) for High Energy Physics (KA2401032) andthe Oﬃce of Advanced Scientiﬁc Computing Research(ASCR) through the Accelerated Research for QuantumComputing Program. MF is supported by the DOE un-der grant DE-SC0010008. This research used resourcesof the Oak Ridge Leadership Computing Facility, whichis a DOE Oﬃce of Science User Facility supported underContract DE-AC05-00OR22725. ∗ [email protected] † [email protected] ‡ [email protected][1] S. P. Jordan, K. S. M. Lee and J. Preskill, Science ,1130-1133 (2012) [arXiv:1111.3633 [quant-ph]].[2] J. B. Kogut and L. Susskind, Phys. Rev. D , 395-408(1975)[3] J. B. Kogut, Rev. Mod. Phys. , 659 (1979)[4] J. B. Bronzan, Phys. Rev. D , 2020-2028 (1985)[5] S. P. Jordan, K. S. M. Lee and J. Preskill, Quant. Inf.Comput. , 1014-1080 (2014) [arXiv:1112.4833 [hep-th]].[6] A. Macridin, P. Spentzouris, J. Amundson andR. Harnik, Phys. Rev. Lett. , no.11, 110504 (2018)[arXiv:1802.07347 [quant-ph]].[7] A. Macridin, P. Spentzouris, J. Amundson andR. Harnik, Phys. Rev. A , no.4, 042312 (2018)[arXiv:1805.09928 [quant-ph]].[8] N. Klco and M. J. Savage, Phys. Rev. A , no.5, 052335(2019) [arXiv:1808.10378 [quant-ph]].[9] D. C. Hackett, K. Howe, C. Hughes, W. Jay, E. T. Neiland J. N. Simone, Phys. Rev. A , 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 , no.3, 032306 (2019) [arXiv:1811.12332 [quant-ph]].[11] M. Kreshchuk, W. M. Kirby, G. Goldstein, H. Beau-chemin 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. , 1343-1346 (1973)[15] H. D. Politzer, Phys. Rev. Lett. , 1346-1349 (1973)[16] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys.B , 199-224 (1985)[17] S. Catani, L. Trentadue, G. Turnock and B. R. Webber,Nucl. Phys. B , 3-42 (1993)[18] R. Bonciani, S. Catani, M. L. Mangano and P. Na-son, Phys. Lett. B , 268-278 (2003) [arXiv:hep-ph/0307035 [hep-ph]].[19] G. Marchesini and B. R. Webber, Nucl. Phys. B ,1-29 (1984)[20] M. Bengtsson and T. Sjostrand, Nucl. Phys. B , 810-846 (1987)[21] C. W. Bauer, W. A. De Jong, B. Nachman and D. Prova-soli, [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. ,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¨om, G. Hagen, T. Papenbrock,A. M. Weiner and M. J. Savage, et al. Phys. Rev. A ,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 , no.8, 083001 (2019) [arXiv:1908.03511 [hep-ph]].[25] A. Roggero, A. C. Y. Li, J. Carlson, R. Gupta andG. N. Perdue, Phys. Rev. D , 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 andM. Saﬀman, [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 , 014006 (2000) [arXiv:hep-ph/0005275 [hep-ph]].[29] C. W. Bauer, S. Fleming, D. Pirjol and I. W. Stewart,Phys. Rev. D , 114020 (2001) [arXiv:hep-ph/0011336[hep-ph]].[30] C. W. Bauer and I. W. Stewart, Phys. Lett. B , 134-142 (2001) [arXiv:hep-ph/0107001 [hep-ph]].[31] C. W. Bauer, D. Pirjol and I. W. Stewart, Phys. Rev. D , 054022 (2002) [arXiv:hep-ph/0109045 [hep-ph]].[32] C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothsteinand I. W. Stewart, Phys. Rev. D , 014017 (2002)[arXiv:hep-ph/0202088 [hep-ph]].[33] C. W. Bauer, A. V. Manohar and M. B. Wise, Phys. Rev.Lett. , 122001 (2003) [arXiv:hep-ph/0212255 [hep-ph]].[34] C. W. Bauer, C. Lee, A. V. Manohar and M. B. Wise,Phys. Rev. D , 034014 (2004) [arXiv:hep-ph/0309278[hep-ph]].[35] A. V. Manohar, Phys. Rev. D , 114019 (2003)[arXiv:hep-ph/0309176 [hep-ph]].[36] T. Becher and M. Neubert, Phys. Lett. B , 251-259(2006) [arXiv:hep-ph/0603140 [hep-ph]].[37] R. Goerke and M. Luke, JHEP , 147 (2018)[arXiv:1711.09136 [hep-ph]]. [38] T. Becher and M. Neubert, Phys. Lett. B , 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 andW. J. Waalewijn, Phys. Rev. D , 114030 (2011)[arXiv:1102.4344 [hep-ph]].[41] R. Kelley, M. D. Schwartz, R. M. Schabinger andH. X. Zhu, Phys. Rev. D , 045022 (2011)[arXiv:1105.3676 [hep-ph]].[42] P. F. Monni, T. Gehrmann and G. Luisoni, JHEP ,010 (2011) [arXiv:1105.4560 [hep-ph]].[43] R. Boughezal, X. Liu and F. Petriello, Phys. Rev. D ,no.9, 094035 (2015) [arXiv:1504.02540 [hep-ph]].[44] I. Moult and H. X. Zhu, JHEP , 160 (2018)[arXiv:1801.02627 [hep-ph]].[45] C. Lin, F. H. Zong and D. M. Ceperley, Phys. Rev. E ,016702 (2001) [arXiv:cond-mat/0101339 [cond-mat.stat-mech]].[46] C. T. Sachrajda and G. Villadoro, Phys. Lett. B ,73-85 (2005) [arXiv:hep-lat/0411033 [hep-lat]].[47] P. F. Bedaque, Phys. Lett. B , 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 , no.7, 074509 (2014) [arXiv:1311.7686[hep-lat]].[49] H. F. Trotter, Proceedings of the American MathematicalSociety 10 (1959) 545.[50] Masuo Suzuki, Communications in MathematicalPhysics 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 , 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 (1995) 487.[58] L. B. Lucy, Astron. J. (1974) 745.[59] W. H. Richardson, J. Opt. Soc. Am. (1972) 55.[60] R. Hicks, C. Bauer, and B. Nachman, Phys. Rev. A ,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. (2017) 180511.[63] M. Gong et al., Phys. Rev. Lett (2019) 110501.[64] K. E. Hamilton et al., [arXiv:2006.01805 [quant-ph]].[65] K. X. Wei et al., Phys. Rev. A (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 (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. (2018)210501. [76] S. Endo, S. C. Benjamin, and Y. Li, Phys. Rev. X (2018) 031027.[77] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev.Lett. (2017) 180509. Simulating eﬀective ﬁeld 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 ﬁeld theory under considerations that our quantumalgorithms can be compared against. The lattice ﬁeld theory and its quantization agree with the analogous treatmentin [5], while the calculations involving the Wilson lines have not appeared anywhere else to our knowledge.

Deﬁnition of the lattice

We consider a lattice of d dimensions, N lattice sites for each dimension and lattice spacing δx . Points on the laticeare indexed by d integers ranging from 0 to N −

1. We denote this set of integers by r = ( r , . . . , r d ) , r i ∈ { , . . . , N − } (S1)with the the positions of the lattice in each dimension given by x r = − x max + r δx , x max = ( n δx , . . . , n δx ) using n = (cid:22) N − (cid:23) . (S2)The corresponding momentum values in each dimension are then given by ( s denotes the same set of integers as r ) p s = − p max + ( s − ∆ ) δp , p max = ( n δp , . . . , n δp ) using δp = 2 πN δx ≡ πL . (S3)Periodic boundary conditions correspond to ∆ = (0 , . . . , −

1) for each wrap around the lattice), correspond to ∆ =(1 / , . . . , / (cid:88) r ≡ d (cid:88) u =1 N − (cid:88) r u =0 . (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 p m = − p n − m +1 for m > , (S5)except for an odd number of lattice sites, in which case the momentum p does not have corresponding partner. (Notethat one can of course translate the lattice and conjugate lattice such that x r = r δx and p s = ( s + ∆ ) δp .)To go between ﬁelds deﬁned on the position spaced lattice to those deﬁned on the momentum spaced lattice oneuses the discrete Fourier transform which is given by φ x = 1(2 π ) d (cid:88) p e − i p · x φ p , φ p = (cid:88) x e i p · x φ x . (S6)Note that we have introduced a short-hand form for the sum over lattice sites as (cid:88) x ≡ ( δx ) d (cid:88) r , (cid:88) p ≡ ( δp ) d (cid:88) s , (S7)where r and s denote the positions on the regular and conjugate lattice, respectively.One has the usual completeness relation1(2 π ) d (cid:88) p e i p · ( x − x (cid:48) ) = δ xx (cid:48) , (cid:88) x e i x · ( p − p (cid:48) ) = (2 π ) d δ pp (cid:48) , (S8)were we have again used a short hand notation δ xx (cid:48) ≡ δx ) d δ rr (cid:48) , δ pp (cid:48) ≡ δp ) d δ ss (cid:48) , (S9)and we have chosen the factor of 2 π to reﬂect the standard choice made in continuum ﬁeld theory.Before we conclude this section, we brieﬂy discuss the relationship between the boundary conditions and the valueof ∆ . From the deﬁnition of the Fourier transformation in Eq. (S6), one can immediately see that transforming aﬁeld by the lattice size in some dimension yields φ x + L n i = e i (2 π ∆ i ) φ x . (S10)Thus, non-integer values of ∆ i give an extra phase in the boundary condition, with half-integer values adding arelative minus sign, giving twisted boundary conditions. The free ﬁeld theory

The free scalar ﬁeld theory on the lattice is given by the Hamiltonian H = 12 (cid:88) x (cid:104) ˙ φ x − φ x [ ∇ φ ] x (cid:105) , (S11)where the sums run over the n possible integers (lattice sites) for each dimension. The d’Alembertian operator isdeﬁned through its action on a ﬁeld, analogously to the second derivative operator in the main text. Performing aFourier transform, one can write this expression as H = 12 1(2 π ) d (cid:88) p (cid:104) ˙ φ p ˙ φ − p + ω p φ p φ − p (cid:105) , (S12)where frequencies are given by ω p = 2 δx (cid:118)(cid:117)(cid:117)(cid:116) d (cid:88) i =1 sin (cid:18) p i δx (cid:19) . (S13)Introducing the creation and annihilation operators as˙ φ p = − i (cid:114) ω p (cid:16) a p − a †− p (cid:17) , φ p = 1 (cid:112) ω p (cid:16) a p + a †− p (cid:17) , (S14)which satisfy the commutation relations (cid:2) a p , a † q (cid:3) = (cid:18) πδp (cid:19) d δ pq , (S15)one ﬁnds for the Hamiltonian H = 1(2 π ) d (cid:88) p ω p (cid:34) a p a † p + 12 (cid:18) πδp (cid:19) d (cid:35) . (S16)Thus, in momentum space the free ﬁeld 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 − , where the occupation number at each lattice site is determined by d integers ≥

0. Thus, a general state is given by | Ψ l k (cid:105) = (cid:79) l k | l k (cid:105) , (S17)with the ground state given by | Ω (cid:105) = (cid:79) l k | (cid:105) . (S18)The energy of the ground state | Ω (cid:105) is given by H | Ω (cid:105) = E | Ω (cid:105) , E = 12 (cid:88) s ω s , (S19)while the energy of excited states are easily determined in terms of their occupation numbers as H | Ψ l k (cid:105) = E l k | Ψ l k (cid:105) , E l k = E + (cid:88) s l s ω s (S20) 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 linesoriginate from the point x = 0, which is labeled by the lattice position n , deﬁned in Eq. (S2). Note that with thisnotation x k with k < n are negative, while x k with k > n are positive. Note that in this section we assume that wehave an odd number of lattice sites, such that we have as many lattice points with x > x <

0. This meansthat 2 n = N − Y n = P exp (cid:34) ig δx N − (cid:88) i = n φ n x i ( t = x i − n ) (cid:35) ,Y † ¯ n = P exp (cid:34) − ig δx n (cid:88) i =0 φ n x i ( t = n − x i ) (cid:35) . (S21)The Wilson line operator is therefore given byT[ Y n Y † ¯ n ] = e − iHδx n exp (cid:2) ig δx (cid:0) φ n x n − φ n x (cid:1)(cid:3) × e iHδx exp (cid:2) ig δx (cid:0) φ n x n − − φ n x (cid:1)(cid:3) × · · · × e iHδx exp (cid:2) ig δx (cid:0) φ δnx n − φ δnx n − (cid:1)(cid:3) × e iHδx exp (cid:2) ig δx (cid:0) φ δnx n − φ δnx n (cid:1)(cid:3) ≡ U n ( − t ) U n − ( δt ) · · · U ( δt ) U ( δt ) . (S22)where we have used that ¯n = − n and x n + m = − x n − m . In the last line we deﬁned U m ( t ) ≡ e iHt exp (cid:2) ig δx (cid:0) φ n x n m − φ n x n − m (cid:1)(cid:3) , (S23)as well as t ≡ n δx , δt ≡ δx . Note that since x n = 0 one has U ( t ) = e − iHt .Using Eq. (S14) we can write U m ( t ) = e iHt (cid:89) s exp (cid:34) ig δx (cid:18) δp π (cid:19) d (cid:0) e i n · p s mδx − e − i n · δp s mδx (cid:1) (cid:112) ω p s (cid:0) a † p s − a p s (cid:1)(cid:35) = e iHt (cid:89) s exp (cid:34) − g δx (cid:18) δp π (cid:19) d (cid:115) ω p s sin( n · p s m δx ) (cid:0) a † p s − a p s (cid:1)(cid:35) = e iHt (cid:89) s D p s ( α p s ( m )) , (S24)with α p s ( m ) = 2 g δx (cid:18) δp π (cid:19) d (cid:115) ω p s sin( n · p s m δx ) . (S25)In the last line we have written the result in terms of the displacement operator that is well know from the theory ofcoherent states D p ( α p ) = e α p a † p − α ∗ p a p , (S26)which satisﬁes the relation D p ( α p ) D p ( β p ) = (2 π ) d D p ( α p + β p ) e i Im( α p β ∗ p ) ( πδp ) d , (S27)and displacement operators acting on diﬀerent momentum modes commute[ D p ( α p ) , D q ( β q )] = 0 . (S28)The action of the time evolution on the displacement operator is e iHt D p [ α p ] = D [ α p ( t )] where α p ( t ) = α p e iω p t . (S29)Combining all this information, one ﬁndsT[ Y n Y † ¯ n ] = (cid:89) s D p s (cid:34) n (cid:88) m =0 α p s ( m ) e i ( n − m ) ω s δt (cid:35) e iφ ps , (S30)with φ p s = (cid:18) πδp (cid:19) d (cid:88) n>m Im (cid:2) α p s ( n ) α ∗ p s ( m ) (cid:3) . (S31)To compute expectation values of this operator, we use (cid:104) n p | D p [ α p ] | Ω (cid:105) = exp (cid:34) − (cid:18) πδp (cid:19) d | α p | (cid:35) α n p p (cid:112) n p ! . (S32)In particular, the magnitude square of the vacuum expectation value is given by (cid:12)(cid:12)(cid:12) (cid:104) Ω | T[ Y n Y † ¯ n | Ω (cid:105) (cid:12)(cid:12)(cid:12) = (cid:89) s exp − (cid:18) πδp (cid:19) d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n (cid:88) m =0 α p s ( m ) e i ( n − m ) ω ps δt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:89) s exp − g ( δx ) (cid:18) πδp (cid:19) d ω p s (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n (cid:88) m =0 e i ( n − m ) ω ps δt sin ( n · p s m δx ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = exp − g (2 π ) d (cid:88) p ω p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) x e − iω p x sin( n · p x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = exp − g (2 π ) d (cid:88) p ω p (cid:88) x ≥ y cos( ω p ( x − y )) sin( n · p x ) sin( n · p y ) , (S33)where the sum over x and y runs over the lattice position and we have deﬁned (cid:88) x ≥ y f ( x, y ) ≡ (cid:88) x f ( x, x ) + (cid:88) x>y f ( x, y ) . (S34)One can verify that this expression reproduces the perturbative result to order g . To that order, one can write (cid:104) Ω | T[ Y n Y † ¯ n ] | Ω (cid:105) = 1 − g (cid:90) L d s (cid:90) L d t [ D F ( n s − n t ) + D F ( ¯n s − ¯n t ) − D F ( n s − ¯n t ) − D F ( ¯n s − n t )]= 1 − g (cid:90) L d s (cid:90) s d t [ D ( n s − n t ) + D ( ¯n s − ¯n t ) − D ( n s − ¯n t ) − D ( ¯n s − n t )]= 1 − g (cid:90) L d s (cid:90) s d t (cid:90) d d p (2 π ) d ω p e − iω p ( s − t ) (cid:104) e i n · p ( s − t ) + e − i n · p ( s − t ) − e i n · p ( s + t ) − e − i n · p ( s + t ) (cid:105) = 1 − g (cid:90) L d s (cid:90) s d t (cid:90) d d p (2 π ) d ω p e − iω p ( s − t ) [cos( n · p ( s − t )) − cos( n · p ( s + t ))]= 1 − g (cid:90) L d s (cid:90) s d t (cid:90) d d p (2 π ) d ω p e − iω p ( s − t ) sin( n · p s ) sin( n · p t ) (S35)Taking the square of this result one ﬁnds (cid:12)(cid:12)(cid:12) (cid:104) Ω | T[ Y n Y † ¯ n ] | Ω (cid:105) (cid:12)(cid:12)(cid:12) = 1 − g (cid:90) L d s (cid:90) s d t (cid:90) d d p (2 π ) d ω p cos( ω p ( s − t )) sin( n · p s ) sin( n · p t ) (S36)which is the continuum limit of Eq. (S33). THE EFFECT OF DIGITIZATION

In this section we provide some details of the eﬀect of digitization on the Hamiltonian of a free massless ﬁeld theory.For this, we consider a (1+1) dimensional ﬁeld theory with 3 lattice sites at points x = − δx , x = 0 , x = δx , (S37)with Hamiltonian H = δx (cid:88) i =0 (cid:104) ˙ φ i − φ i [ ∇ φ ] j (cid:105) , (S38)The momenta of the 3-site lattice are given by p = − π δx , p = − πδx , p = πδx . (S39)The corresponding frequencies are obtained from Eq. (12) ω = 2 δx , ω = ω = 1 δx . (S40)The exact eigenvalues of this Hamiltonian are now easily obtained. The energy of the ground state is given by halfof the sum of frequencies, and therefore E = 2 / δx . The excited states are simply given by the ground state energyincreased by up to n φ additions of each frequency. Thus, the lowest lying states are1 × δx , × δx , × δx , × δx , . . . (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 arereproduced with an accuracy of 10 − , for n q = 3 with an accuracy of 10 − and for n q = 2 with an accuracy of about10 %. This is in agreement with the ﬁndings of [8].In Fig. 2 we show the dependence of the vacuum expectation value Y X on the coupling g on the number of qubitsper lattice site. One can see that already for n Q = 2 qubits per lattice site one gets a very good approximation tothe exact lattice result, and for n Q = 3 the results including digitization are indistinguishable from the exact resultsif 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 starsfor n Q = 2 and 3. One can see that using the KW approximation introduces an additional error, with the resultingemission probability reduced compared to the exact result, but that the resulting diﬀerence is sub-percent level bythe time n Q = 3. For these plots only, we compare the KW approximation originally presented in [52] and used inthe rest of the text with a modiﬁed form discussed in Section . Eigenvalue F r a c t i o n a l E rr o r i n E i g e n v a l u e n Q = 2 n Q = 3 n Q = 4 FIG. 1. The fractional error in the 15 lowest-lying eigenstates on an N = 3 lattice. The dependence on the number of qubitsper 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 displayseigenvalues greater (smaller) than their undigitized values. For n Q = 3, all 2 particles states are already reproduced at betterthan 1 % accuracy. (If modes with (cid:38)

10 % deviations from their continuum values are excluded, this increases to all 5 particlestates.) The exponential improvement with number of qubits per site is clearly visible. g | X | T [ Y n Y n ] || X = X = p Analytic calculationDigitized calculations n Q = 1 n Q = 2 n Q = 3 FIG. 2. The eﬀect of ﬁeld digitization on transition rates in the presence of the Wilson line for N = 3. The blue solid lines arethe analytical result for transitions to the vacuum (Ω) and the lowest-lying one one-particle state ( p ). Digitization limits thenumber of states the ﬁeld 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, theexact ground state is digitized and evolved. CONSEQUENCES OF THE NON-LOCAL NATURE OF WILSON LINES SPECIAL CONSIDERATIONSIN (1+1) DIMENSIONS

The two Wilson lines Y n and Y † ¯ n describe how the dynamical soft degrees of freedom in the eﬀective theory (bosonsin our case) interact with the fermions that move at the speed of light through the system. These Wilson lines containan integral over the scalar ﬁelds along the world lines of the fermions, going from point 0 to ∞ along the directions n µ = (1 , n ) and ¯ n µ = (1 , − n ), and are therefore non-local in nature, with the non-locality extending all the wayto inﬁnity. This non-locality along the light-like directions gives rise to the well known collinear singularities in thematrix elements. An important consequence of the non-local nature of the operators and their inﬁnite extent is thattheir lattice implementation necessarily makes them directly sensitive to the lattice volume L and separation δx . Thefact that the Wilson lines can only extend up to the edge of the lattice provides a regulator for the angle between themomentum of bosons and the direction 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 deﬁnition of the Wilson line g | X | T [ Y n Y n ] || X = X = p n Q = 2 Analytic calculationDigitized calculations exact correlations Kitaev Webb Kitaev Webb (modified) g | X | T [ Y n Y n ] || X = X = p n Q = 3 Analytic calculationDigitized calculations exact correlations Kitaev Webb Kitaev Webb (modified)

FIG. 3. The eﬀect of approximated state preparations on transition rates in the presence of the Wilson line for N = 3 latticesites. The blue solid lines are the analytical result for transitions to the vacuum (Ω) and the lowest-lying one one-particle state( p ). We compare the digitization on the exact state with the original and modiﬁed 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). operators. The presence of the ﬁnite lattice produces terms L n · p and L ¯ n · p , regulating divergences as n · p and ¯ n · p tend to zero. On the other hand, for ﬁnite lattice spacing one has ω − | p | ∼ | p | δx , such that the UV regulator isalso prohibiting the momenta of the bosons to go on-shell, which in turn also regulates the collinear divergence. Thisimplies that collinear divergences in matrix elements give rise to mixed UV-IR divergences, which would be absent inthe continuum theory (even though the continuum theory also contains mixed UV-IR divergences coming from thethe limit n · p → p → ∞ ). Once infrared and collinear (IRC) safe physical observables are being calculated, inwhich all collinear divergences cancel, these lattice artifacts will also cancel, leaving only pure UV divergences relatedto the normalization of operators.There is an additional subtlety in a (1+1) dimensional ﬁeld theory as is used in the numerical results of this paper.With only a single spatial dimension one can of course only deﬁne a single direction. This means that in such a theoryall momenta have to be aligned with the direction n or ¯ n , which means that every momentum satisﬁes n · p = 0 or¯ n · p = 0 and therefore gives rise to a collinear divergence. In other words, each state of the quantum ﬁeld theory isactually collinear to the directions n µ and ¯ n µ . It should be recalled that the underlying short distance process we aredescribing is the production of two energetic back-to-back particles. The diﬀerential distribution in the total energyof this process is described by short distance physics captured by the hard function, while the jet function and thesoft function discussed in this paper describe the distribution of collinear and soft radiation originating from theseenergetic 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 inclusiveof any additional radiation. Therefore, the only IRC safe observables are those deﬁned by the hard kinematics, whichare completely inclusive over the soft radiation calculated here. Since the Wilson line operator is unitary, the squareof the completely inclusive operator is unity and therefore trivial. In other words, any non-trivial observable in thistheory is by construction IRC unsafe, and we therefore choose the simplest such observables, namely the expectationvalue of the Wilson line operator between the vacuum states, and the transition from the vacuum to a state with asingle momentum mode. CONSTRUCTING THE QUANTUM CIRCUITSRepresentation of the Hilbert space and relationship of qubit representation to ﬁeld values

To represent the Hilbert space of the massless scalar ﬁeld theory we discretize the spatial coordinates on a latticeof size L and N sites per dimension. The value of the scalar ﬁeld on each lattice site is digitized and represented by n φ = 2 n Q distinct values. One can write the ﬁeld as φ ( k ) i = − φ max + k δφ for k = 0 , . . . , n φ − , δφ = 2 φ max n φ − . (S42)Representing the integer | k (cid:105) i through its binary representation of the n Q qubits at lattice site i , we can deﬁneˆ φ i | k (cid:105) i = φ ( k ) i | k (cid:105) i . (S43)The value of φ max should be chosen to minimize the minimize the error due to the digitization and depends on theHamiltonian implemented on the lattice. The integer state | k (cid:105) i is represented as usual through the bitstring of the n Q qubits at the given lattice site | k (cid:105) i = | q (cid:105) i · · · (cid:12)(cid:12) q n Q − (cid:11) i . (S44)The full Hilbert space is then represented through the tensor product of the states | k (cid:105) on each lattice site | Ψ (cid:105) = | k (cid:105) · · · | k (cid:105) N d . (S45)Our explicit circuit constructions in this paper will only use a single spatial direction, such that we use | Ψ (cid:105) = | k (cid:105) · · · | k (cid:105) N . (S46) The free Hamiltonian

The construction of the Hamiltonian of free massless scalar ﬁeld theory follows previous work [6–8]. One can easilyconvince oneself that the operator ˆ φ i can be written through its action on the n Q qubits at each lattice siteˆ φ i = n Q − (cid:88) j =0 j ˆ σ ( j ) z,i , (S47)where the operator ˆ σ ( j ) z,i is a single σ 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 ﬁeldis the conjugate ﬁeld π i ≡ ˙ φ i , and one can write H = H π + H φ (S48)with H π = δx (cid:88) i π i , H φ = δx N − (cid:88) i =0 φ i [ ∇ φ ] i . (S49)As discussed before, the φ i [ ∇ φ ] i operator only requires nearest neighbor interactions on the lattice. The timeevolution is then written in terms of the Suzuki–Trotter formula (we give the ﬁrst order expression here) (cid:2) e − iHt (cid:3) n = (cid:104) e iH π t/n e iH φ t/n (cid:105) n . (S50)To construct the exponential of the H π operator, we use that the φ and π are related by a Fourier transform onthe φ register. Thus, we can write e iH π t = QFT − e iδx tφ i QFT , (S51)where QFT denotes the (symmetrized) Quantum Fourier Transform, which was discussed in [8]. We do not repeatits circuit here.Given this, on needs to ﬁnd a circuit representation exp[ iθφ i φ j ] for general i and j , from which one can constructboth the circuits for exp[ iH π t ] and exp[ iH φ t ].Using Eq. (S47), one can writeˆ φ i ˆ φ j = (cid:34) n Q − (cid:88) l =0 l ˆ σ ( l ) z,i (cid:35)(cid:34) n Q − (cid:88) k =0 l ˆ σ ( k ) z,j (cid:35) = n Q − (cid:88) l =0 n Q − (cid:88) k =0 ( l + k ) σ ( l ) z,i σ ( k ) z,j . (S52)This allows us to write exp (cid:104) iθ ˆ φ i ˆ φ j (cid:105) = n Q − (cid:89) l =0 n Q − (cid:89) k =0 exp (cid:104) i ( l + k ) θ σ ( l ) z,i σ ( k ) z,j (cid:105) . (S53)The action exp (cid:104) i ( k + l ) θ σ ( l ) z,i σ ( k ) z,j (cid:105) is equal to exp (cid:2) i ( k + l ) θ (cid:3) if qubits q ( l ) i and q ( k ) i are equal and exp (cid:2) − i ( k + l ) θ (cid:3) if theyare opposite. Thus, it can be implemented by the circuit | l (cid:105) i • •| k (cid:105) i e − i ( k + l θ ) Z 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 − / Ground state preparation

The ground state of a massless scalar ﬁeld theory is given by a multivariate Gaussian | Ψ (cid:105) = exp (cid:20) −

12 ˆ φ i G ij ˆ φ j (cid:21) | k (cid:105) · · · | k (cid:105) n (S54)While Qiskit provides a function to generate an aribtrary Gaussian multivariate distribution, the number of gatesin the resulting circuit scale exponentially with the number of qubits in the system. However, an algorithm withpolynomial scaling was derived by Kitaev and Webb [52]. While this algorithm does not produce the exact multivariatedistribution 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 correlationmatrix in terms of a diagonal matrix D and a lower unit-triangular matrix L (a lower triangular matrix with 1’s onthe diagonal) G = LDL † (S55)An arbitrary multi-dimensional Gaussian distribution can then be created by generating a series of uncorrelatedGaussian distributions according to the diagonal matrix D , and then applying shear matrices L through a remappingof the basis states of the Hilbert space. For details of this algorithm, we refer the reader to the original paper. InSection we mentioned a modiﬁed version of the KW procedure, which only rounds the shear matrices applied to theuncorrelated states to the nearest digitized ﬁeld value after applying the the full shear matrix rather than after everyindividual shearing operation, as was originally proposed in [52]. While requiring more ancilla qubits for memory,this does not aﬀect the polynomial scaling of the approximation and results in exponentially greater ﬁdelity with theexact ground state for the n Q ≥ 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 functionalityof 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 determinedthrough the successive application of the operator exp[ ig δx φ n ] and Hamiltonian evolution. Given the circuit for0Hamiltonian evolution derived above, one therefore only needs the circuit for the exponential of a single ﬁeld. This iseasily derived using a similar derivation to the Hamiltonian case, and one writesexp[ iθ ˆ φ i ] = n Q − (cid:89) j =0 exp (cid:104) i j θσ ( j ) z,i (cid:105) , (S56)which can easily be written in circuit form | (cid:105) i e − iθZ ... · · · ... | n Q − (cid:105) i e − i ( nq − θZ VALIDATION OF THE CIRCUITS FOR THE HAMILTONIAN

As a validation of the quantum circuits we ﬁrst check the implementation of the free ﬁeld theory (ground statepreparation and time evolution). In particular, we check to what degree the ground state is an eigenstate of theHamiltonian, and how well the energy of the ground state agrees with the known analytical value. We begin bycomputing the overlap f ( t ) = (cid:12)(cid:12) (cid:104) Ω | (cid:2) e − iHt (cid:3) n | Ω (cid:105) (cid:12)(cid:12) , (S57)where (cid:2) e − iHt (cid:3) n is the Trotterized Hamiltonian with n steps given in Eq. (S50). This is implemented with the circuit / U state (cid:2) e − iHt (cid:3) n U † state and the function f ( t ) is obtained by the fraction of measurements where all qubits are back the initial | (cid:105) state. Ifthe ground state is indeed an eigenstate of the Hamiltonian, and the Trotterized Hamiltonian is equal to the fullHamiltonian, this function should be identically 1, and any deviation from this value should be due to Trotterizationerrors or because the state | Ω (cid:105) not being equal to the true ground state. In Fig. 5 we show this overlap for the exactground state on the left and for the KW ground state on the right. The result on the left conﬁrms that with moreTrotter steps we approach unity, while the result on the right shows that even for a large amount of Trotter steps theKW approximation leads to deviations from unity.While this measurements tests to what degree the ground state is an eigenstate of the Trotterized Hamiltonian, itcan 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 / U state (cid:2) e − iHt (cid:3) n U † state H • H can measure the energy of the ground state. The fraction of measurements with all qubits in the | (cid:105) state is given by f ctr ( t ) = 14 (cid:12)(cid:12) (cid:104) Ω | (cid:2) e − iHt (cid:3) n | Ω (cid:105) (cid:12)(cid:12) , (S58)which is sensitive to the energy E Ω . One can see that using the digitization of exact ground state and suﬃcient Trottersteps one produces the analytically expected dependence on the ground state energy up to the small diﬀerence inperiod due to the shift in ground state energy due to digitization. Conversely, with the KW states one also sees asmall 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-processingstep. We prepare all 2 n qubit possible states | i (cid:105) and measure the frequency of observing the state | f (cid:105) . These probabilities1 t || T [ Y n Y n ] || exact ground state n = 1 n = 5 t exact t || T [ Y n Y n ] || KW approximation n = 1 n = 5 t exact FIG. 4. The value of (cid:12)(cid:12) (cid:104) Ω | (cid:2) e − iHt (cid:3) n | Ω (cid:105) (cid:12)(cid:12) for two diﬀerent values of n . The blue line shows the result for n = 1, while the redline shows the result for n = (cid:100) t/ . (cid:101) . 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 exactground state, the deviation from unity is only due to the Trotterization of the Hamiltonian, and gets smaller with more Trottersteps (and more quantum gates). Since the KW approximation is not an eigenstate of the full Hamiltonian, the deviation fromunity is due to both this approximation and the Trotterization. t || + T [ Y n Y n ] || / exact ground state n = 1 n = 5 t exact t || + T [ Y n Y n ] || / KW approximation n = 1 n = 5 t exact FIG. 5. The value of (cid:12)(cid:12) (cid:104) Ω | (cid:2) e − iHt (cid:3) n | Ω (cid:105) (cid:12)(cid:12) / n . The blue line shows the result for n = 1, while thered line shows the result for n = (cid:100) t/ . (cid:101) . Measurements from a noiseless simulation of 512 shots/point are overlayed. On theleft, 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. are tabulated in a 2 n qubit × n qubit matrix called the response matrix R . The observations from the main experimentare represented as a vector m i and then readout error mitigation proceeds iteratively [54, 57–59]: t n +1 i = (cid:88) j Pr(truth is i | measure j ) × m j = (cid:88) j R ji t ni (cid:80) k R jk t nk × m j , (S59)where t i = 1 / n qubit is the uniform prior and the number of iterations n is chosen to be 20. The result is notsensitive to small variations in these choices. The iterative procedure in Eq. S59 avoids pathologies from otherforms of regularized matrix inversion and can be further improved with variations like readout rebalancing [60] andsubexponential approximations [61–65]. Matrix inversion approaches from the quantum simulators pyQuil [66] (byRigetti), Cirq [67, 68] (by Google), and

XACC [69–71] and the least squares method from

Qiskit by IBM [53, 72] (seealso Ref. [73, 74]) produce similar results [54]. In current IBMQ machines, the dominant gate noise is from multiqubit2operations, which are only the CNOT gates in our experiments. A common model for CNOT errors is the quantumdepolarizing noise channel. This error is mitigated by systematically increasing the noise and then extrapolating tozero noise. In particular, we remove the leading order depolarizing noise by repeating the nominal measurement witha circuit where every CNOT gate is replaced with 3 CNOT gates. Following Ref. [55], we then set the mitigatedresult as 3 / //