Simulating collider physics on quantum computers using effective field theories
SSimulating collider physics on quantum computers using effective field 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 field theory over a wide range of energies requiresexceptionally large quantum computing resources. Yet for many observables in particle physics,perturbative techniques are sufficient to accurately model all but a constrained range of energieswithin the validity of the theory. We demonstrate that effective field theories (EFTs) provide anefficient 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 first 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 field 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 field theo-ries [1]. The main technique involves disretizing the spa-tial degrees of freedom by introducing a lattice [2–4], anddigitizing the field values at a given lattice point [5–13].This turns the uncountably infinite dimensional Hilbertspace of standard quantum field theories into a finite 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 finite volume of space introduceupper and lower cutoffs to the energies E over which theresulting lattice field theory is a good approximation forthe continuum. In particular, one finds1 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 fixed order per-turbation theory, and high precision can be reached withexisting techniques. Physics at lower energies introducesnew significant 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 fixed order per-turbative calculations are of little use to describe lowerenergy effects 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 effects 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 first principles.In this paper we show that this is indeed possible byusing effective field 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 first 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 simplifies the description of theinitial and final state and significantly 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 final 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 Effective 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 coefficient 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 thedifferent 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 final state, with the dynamics described by the EFTHamiltonian. The final 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 effects arising from the high multiplicity final 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 different 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 field theory (with any degrees of freedomthat only contribute to short distance effects 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 effect is captured by a so-called Wilsonline, which describes the interaction of a charge movingalong a fixed 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 fixed 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 fields [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 field and n µ describes the di-rection of the light-like direction n µ = (1 , n ) with n = 1such that n µ n µ = 0. Thus, the gauge field is evaluatedon a path going from the origin to spatial infinity alongthe world-line characterized by the direction n µ . Thedynamics of the soft gauge field 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 final state. If re-producing multi-particle weakly coupled final states, | X (cid:105) will contain a fixed number of gauge momentum modeswith given momenta k µi . To compute the soft function fora given observable one needs to sum over all possible finalstates that can contribute to this observable. Traditionalperturbative calculations [38–44] can only compute thesematrix elements for final 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 field theory wouldprovide the full non-perturbative result for the soft ma-trix element. To evaluate the matrix element on a quan-tum computer one first needs to define 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 simplification of working within an EFT.To be precise, we consider a massless field theory in(1+1) dimensions with Hamiltonian and Wilson lines de-fined 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 effects 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 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 − δ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 defined through its ac-tion on a field 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 different field values. For each lattice point, the possiblefield 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 fieldsis 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 first 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 field 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 fieldoperator. 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 defined.Ultimately, to make contact with the continuum fieldtheory 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 offields φ 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 field theory is a multi-variate Gaussian distribu-tion, which can be created using the approach of Kitaevand Webb (KW) [52]. To identify states of definite 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 first 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 defined 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 finally 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 simplifies 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 fi-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 fordifferent final 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 effect 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 field 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 field 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 differs 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), Office of Science under contractDE-AC02-05CH11231. In particular, support comesfrom Quantum Information Science Enabled Discovery(QuantISED) for High Energy Physics (KA2401032) andthe Office of Advanced Scientific 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 Office 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. Saffman, [arXiv:2011.04097 [hep-th]].[27] T. F. Stetina, A. Ciavarella, X. Li and N. Wiebe,[arXiv:2101.00111 [quant-ph]].[28] C. W. Bauer, S. Fleming and M. E. Luke, Phys. Rev. D , 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 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 quantumalgorithms can be compared against. The lattice field 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.
Definition 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 fields defined on the position spaced lattice to those defined 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 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 valueof ∆ . From the definition of the Fourier transformation in Eq. (S6), one can immediately see that transforming afield 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 field theory
The free scalar field 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 isdefined through its action on a field, 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 finds 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 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 − , 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 , defined 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 defined 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 satisfies 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 different 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 findsT[ 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 defined (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 finds (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 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 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 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 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 difference 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 modified 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 effect of field 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 field can take on for n Q = 1 (red) to 2, n Q = 2 (orange) to 4 and n Q = 3 green to 8. In each case, 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 effective 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 fields 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 infinity. 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 infinite 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 definition 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 effect 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 modified KW approximation (see text for details)which can be prepared with polynomial resources in number of qubits. The plots show results for n Q = 2 (left) and n Q = 3(right). operators. The presence of the finite 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 finite 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 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 theoryall momenta have to be aligned with the direction n or ¯ n , 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 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 differential 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 defined 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 field values
To represent the Hilbert space of the massless scalar field theory we discretize the spatial coordinates on a latticeof 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 φ ( 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 defineˆ φ 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 field 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 fieldis the conjugate field π 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 first 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 find 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 field 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 modified version of the KW procedure, which only rounds the shear matrices applied to theuncorrelated states to the nearest digitized field 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 affect the polynomial scaling of the approximation and results in exponentially greater fidelity 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 field. 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 first check the implementation of the free field 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 confirms 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 sufficient Trottersteps one produces the analytically expected dependence on the ground state energy up to the small difference 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 different 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 / //