Accelerating lattice quantum field theory calculations via interpolator optimization using noisy intermediate-scale quantum computing

The only known way to study quantum field theories in nonperturbative regimes is using numerical calculations regulated on discrete space-time lattices. Such computations, however, are often faced with exponential signal-to-noise challenges that render key physics studies untenable even with next generation classical computing. Here, a method is presented by which the output of small-scale quantum computations on noisy intermediate-scale quantum era hardware can be used to accelerate larger-scale classical field theory calculations through the construction of optimized interpolating operators . The method is implemented and studied in the context of the 1 þ 1 -dimensional Schwinger model, a simple field theory which shares key features with the standard model of nuclear and particle physics.

Numerical approaches to quantum field theory are the only known way to make predictions for a wide range of physical quantities from the standard model of particle physics, our best current theory of nature at the smallest scales.Standard model calculations of nuclear physics processes-such as those needed to interpret experiments using nuclei as targets-are particularly challenging.In particular, the strong-interaction component of the standard model, which is encoded in the theory of quantum chromodynamics (QCD), cannot be approached analytically at the relevant energy scales.The only first-principles approach to QCD at these scales is numerical: a discretized form of the QCD equations can be solved using supercomputers through Monte Carlo integration on a finite four-dimensional grid representing space-time [1,2].This technique, named lattice quantum field theory (LQFT), plays an important role in modern particle and nuclear physics and has been essential in testing the standard model against precise measurements of the decays and interactions of particles at frontier machines such as the Large Hadron Collider [3,4].Calculations of nuclei, however, are limited by exponentially bad scaling of computational cost with the atomic number of the system being studied.Using current methods, direct studies of nuclei with tens of nucleons, as relevant to diverse physics programs from direct searches for dark matter to neutrino physics, will remain intractable, even with the advent of exascale classical computing in the next years; progress on this front will require a revolutionary approach, and there is great interest in the potential applications of quantum computing to overcome this challenge [5,6].Hybrid methods coupling classical and quantum computing offer a natural pathway to exploit quantum computation despite the small number of qubits, sparse qubit connectivity, lack of error correction, and noisy quantum gates that are hallmarks of current and near-term quantum computing in the noisy intermediate-scale quantum (NISQ) era [7].
A significant contribution to the computational cost of LQFT studies could be eliminated by the construction of optimized interpolating operators, corresponding in broad terms to approximations to the quantum wave function of the desired state.Precisely, to determine matrix elements of interest in some state in a LQFT computation, such as those describing an interaction or decay process, correlation functions are calculated which encode the creation, interaction, and annihilation, of the state in question.These correlation functions, however, receive contaminating contributions from the many other states with the quantum numbers of the state of interest.In order to reliably extract the desired piece, the contributions from all of these unwanted higher-energy states must be suppressed.Typically, this is achieved via an evolution in the Euclidean time of the calculation; the unwanted states are exponentially suppressed by the energy gap to the ground state at large times, but at the cost of an exponential growth in the statistical noise of the Monte Carlo sampling used in the computation (and thus computational cost).By using optimized interpolating operators for state creation and annihilation, constructed to have significant overlap onto the state of interest, this Euclidean time evolution, and thus exponential growth in noise, can be reduced.In this Letter, it is demonstrated for the 1 þ 1dimensional Schwinger model how one can construct such interpolating operators for classical LQFT calculations using small-scale quantum computation.Ultimately, the extension of this approach to the more complex theory of QCD, along with advancement in quantum hardware, could enable a significant acceleration of LQFT computations for nuclear physics.
The Schwinger model.-TheSchwinger model [8], which describes the theory of quantum electrodynamics in one space and one time dimension, is a prototypical lattice gauge theory that shares a number of key features with QCD.This model thus provides a simplified framework to test new algorithms and approaches to LQFT studies.The theory describes fermions as a two-component spinor field, with mass m, coupled via charge g to an electromagnetic field, A μ .A discretized formulation of the 1 þ 1D Schwinger model can be defined on a staggered space-time lattice via the Kogut-Susskind prescription [9], with the staggered fermion field operators denoted by φðx n Þ ¼ φn .
In temporal gauge (A 0 ¼ 0), the single spatial component of the gauge field is encoded on links connecting adjacent staggered sites x n and x nþ1 .The associated electric flux operators can be defined in terms of operators ln , which, together with corresponding raising and lower operators, act on the space of links connecting sites: The eigenvalue l n describes the value of the electric flux at the link connecting sites n and n þ 1.
Combining the link space with fermionic occupation numbers, a complete Fock space of states in this theory can be expressed as fj⃗ n; ⃗ lig.On a lattice with N staggered sites (N=2 spatial sites), the Hamiltonian of this theory can be expressed in terms of these operators as In terms of the spatial lattice spacing a, the parameters w ¼ 1=2a and J ¼ g 2 a=2 are identified with the temporal scale and coupling strength respectively.This theory has a simple spectrum of low-lying states of conserved parity quantum number; the first excited state is odd-parity, interpreted as the massive photon (the lightest vector meson), while the second excited state is the evenparity scalar e þ e − "meson." Classical computations of ground-state energies.-Usingclassical computation, energy levels of the Schwinger model can be obtained using standard Monte Carlo (MC) methods.Here a local Hamiltonian MC method is studied [10]; details of the application of this approach to the 1 þ 1D Schwinger model are given in Ref. [11].In this formalism, ground-state energy levels can be determined by the analysis of correlation functions GðτÞ, defined in terms of the expectation values of interpolating operators Ôðx; τÞ, which are constructed to create and annihilate states with the quantum numbers of a target state of interest at some Euclidean position (x, τ): Here, a state is created at some initial spatial position and time ðx ¼ 0; t ¼ 0Þ, and annihilated τ Euclidean time steps later.The sum over x projects onto the zero-momentum sector.This correlation depends on the energy gaps between the ground state (vacuum) of the system jΩi and the tower of excitations coupled to the vacuum through Ô: Numerically, the energy gap to the lowest state of interest is determined from the asymptotic value of the effective mass function: Interpolating operators for lattice field theories can be constructed by inspection, and often the simplest operators which have the quantum numbers of the state of interest are chosen.For the Schwinger model, the lowest-energy excitation is described by the lightest vector meson, V − , a massive photon.An interpolating operator for this state can be constructed as for odd-parity meson states in staggered lattice formulations of QCD [12]: Physically, the operator ÔV creates an e þ e − pair on sites x n−1 and x n , and a second e − e þ pair on sites x n and x nþ1 .
The relative minus sign between the two terms ensures that the construction has odd parity.Improved interpolating operators for LQFT calculations can be constructed classically via a variational approach: rather than a single interpolating operator, a set of operators with the same quantum numbers is chosen, and the resulting system is diagonalized via a generalized PHYSICAL REVIEW LETTERS 124, 080501 (2020) eigenvalue problem (GEVP) to achieve an optimized ground-state energy extraction in that sector [13][14][15].Despite tremendous success [16], the classical variational method is computationally expensive, scaling quadratically with the size of the basis.In particular, using a large variational basis requires high-statistics numerical calculations to ensure a nonsingular covariance matrix.This method thus remains intractable for many studies, such as QCD calculations of nuclei [5].In this work, an alternative variational approach to interpolating operator construction is explored, in which the expense of a classical variational method is replaced with small-scale computations on quantum hardware.
Quantum approaches to Schwinger model dynamics.-Severalyears ago, the first experimental demonstration of a digital quantum simulation of a lattice gauge theory was achieved by realizing the Schwinger model on a few-qubit trapped-ion quantum device [17][18][19][20].Recently, the twospatial-site Schwinger model has also been studied on IBM's superconducting quantum hardware [21].These quantum calculations use an equivalent formulation of the Schwinger model based on bosonic degrees of freedom, related by a Jordan-Wigner transformation [22] to the formulation described here.In both demonstrations, the ground state energy level of the theory was extracted using the variational quantum eigensolver (VQE) method [23].In the VQE approach, a sequence of unitary operators U ðiÞ ð ⃗ θÞ, implemented as a sequence of one-and two-qubit gates, is tuned using variational parameters θ to transform an initial, easy-to-prepare state j0i into an approximation of the ground state of the system in a given symmetry sector, jGi: From this construction, an approximate value of the ground-state energy of the system can be calculated.In the study of Klco et al. [21], which presents a formulation of the Schwinger model on quantum hardware which has a natural relation to classical approaches to the theory, explicit electric flux degrees of freedom are retained in the basis of states studied using the VQE approach.To render the Schwinger model Hilbert space in this description finite thus requires truncating the possible values of flux on each link.This can be achieved by enforcing for some choice of truncation fΛ 2 ; Λ2 g; harsher truncations result in larger systematic uncertainties in the determined energy level, but require fewer qubits for simulation.Naturally, the small system sizes accommodated by NISQ-era quantum-computing hardware result in additional finite-size systematic uncertainties in numerical studies.
Given the present status of quantum computation, scaling these studies to determine ground-state energies of systems of physical interest is a long-term challenge.Using VQE calculations in hybrid approaches to accelerate classical LQFT computations, which can be more easily scaled at the present time, however, offers the potential of near-term exploitation of these new tools.
Quantum-improved interpolating operators.-Here,it is proposed to use the information encoded in a variationally obtained approximation to a ground-state wave function to construct an improved interpolating operator for use in classical LQFT computations of that state.This is achieved via a two-step approximation process: first, VQE calculations are used to yield approximate representations of the wave functions of both the dynamical vacuum and the ground state of the symmetry sector of interest; second, a linear combination of operators is optimized to maximize the transition matrix element between the vacuum and the state of interest.This procedure can be considered as analogous to the classical variational approach to interpolating operator construction, and may have significant advantages in a near-future era of efficient small-scale quantum computation.In particular, as is demonstrated here, one can undertake a quantum computation in a truncated Hilbert space on a small lattice volume, and use this to determine a reduced basis of operators to compute classically on the full state space.
As an explicit example, for the Schwinger model one might study the first excited state of the theory (the odd-parity massive photon).As outlined in Eq. ( 8), VQE computations can provide approximate descriptions of both this state, jV − ; ⃗ θ 0 i, and the vacuum of the theory jΩ; ⃗ θi.Acting with a fixed set of interpolating operators f Ôk g on the ground state provides a variational basis that can be used to approximate the massive photon: The operators Ôk can be defined, for example, in terms of field and link operators φn , ln , LAE n , and should be constructed to transition between the vacuum and the symmetry sector of interest.With an estimate for the matrix elements hV − ; ⃗ θ 0 j Ôk jΩ; ⃗ θi, a classical computer can be used to optimize the overlap jhV − ; ⃗ θj Ṽ− ij (with appropriate normalization) with respect to ⃗ α.This optimization defines an improved interpolating operator which has suppressed overlap onto excited states and can be implemented in a classical Euclidean MC calculation.Within the formalism of Ref. [21], the relevant transition operators can be evaluated by projective measurements PHYSICAL REVIEW LETTERS 124, 080501 (2020) [20,[22][23][24] for operators spanning the two symmetry sectors of interest.Implementation for the 1 þ 1D Schwinger model.-Theproposed approach to interpolating operator construction is implemented in the context of the 1 þ 1D Schwinger model.The parameters of the Hamiltonian are chosen to match those in recent studies of this model on quantum devices [21]: J=ω ¼ 5=3 and m=ω ¼ 1=6, and the temporal lattice spacing is set as Δt ¼ a=5 based on previous Monte Carlo studies [11].
For direct comparison, the mass of the lightest vector meson in this model can be extracted using the local Hamiltonian MC method [10,11].Here, a system with 4 spatial sites (8 staggered sites), and 40 temporal (80 staggered sites) is studied, with the standard interpolating operator defined in Eq. ( 7) used as a benchmark; the effective mass from a numerical computation with 10 7 configurations sampled is shown in Fig. 1.The effective mass with this operator is also reconstructed exactly from the numerical diagonalization of Eq.This diagonalization is performed with respect to a truncation of the electric fluxes, fΛ 2 ; Λ2 g ¼ f1; 8g, that is chosen to encode the same physics that is being sampled within the Euclidean MC setup.
A quantum-improved interpolating operator, constructed via the approach proposed here, can also be investigated.First, the even and odd-parity ground states are obtained from the corresponding exact solution, which is used as a proxy for a VQE in this proof-of-principle demonstration.A linear combination of up to 6 operators defined in terms of electric flux operators l [Eq.( 1)], detailed in the Supplemental Material [25], is then optimized to maximize the overlap with the target state.Using the exact solution as a proxy for a VQE, in this proof-of-principle work the matrix elements hV − ; ⃗ θ 0 j Ôk jΩ; ⃗ θi were calculated on a classical computer.Of course, this approach is not scalable to even modestly larger systems; ultimately these matrix elements must be evaluated on quantum hardware.This improved operator, of the form Eq. ( 11), is then used to construct a correlator [Eq.( 4)] which has significantly improved overlap with the lowest-lying negative-parity state.Figure 1 displays the corresponding effective masses obtained for both the Monte Carlo ensemble and those constructed from exact diagonalization.The effective mass obtained from a MC computation with the quantumimproved interpolating operator leads to a far better constrained mass extraction at the same statistics than the benchmark operator.Moreover, the quantum-improved operator produces an exact effective mass curve which is indistinguishable, on the scale of Fig. 1, from that obtained via an exact version of the classical generalized eigenvalue program.
Naturally, in a true quantum computation the VQE wave functions will be only approximately determined, with statistical uncertainties on the variational parameters f ⃗ θ; ⃗ θ 0 g, and the corresponding matrix elements hV − ; ⃗ θ 0 j Ôk jΩ; ⃗ θi will be similarly limited in measurement precision.The effect of such uncertainties on the definition of the quantum-improved interpolating operator ÔVQE is investigated by taking a 0.05 radian error on the variational parameters, corresponding to approximately 15% error on the coefficients of expansion.These uncertainties are roughly equivalent to the fidelities obtained in recent studies of the Schwinger model using modern quantum hardware [20].Figure 1 shows the corresponding uncertainty in the VQE-improved correlation function.Importantly, even with an imperfect quantum computation, the quantum-improved operator offers improved isolation of the ground state compared with the benchmark operator.
As described, it is evident that a quantum calculation has the potential to improve the results of a conventional MC computation.Nevertheless, the scaling of quantum hardware to encode larger Hilbert spaces is expected to be an enduring problem.It would therefore be of tremendous value if quantum simulations in a (significantly) truncated Hilbert space can still produce improved operators for a Euclidean MC calculation-which can more readily be scaled to larger systems.Restricting the gauge space, by truncating the gauge link variables, and limiting the spatial volume, are two natural ways to reduce the Hilbert space for simulations on quantum hardware.6)] constructed exactly via diagonalization (solid lines), and obtained from MC computations (open points), for the lightest vector meson.The green circles and red squares denote results obtained using the benchmark [Eq.( 7)] and quantum-improved [Eq.(11)] operators, respectively.The black dashed line shows the result obtained solving an exact form of the classical GEVP using the same basis of interpolating operators, while the purple solid band shows the exact result for the optimized operator, given some uncertainty on the VQE input into its construction as described in the text.Note that the uncertainty is suppressed as the Euclidean separation between operators is increased.The large time behavior is a consequence of the backward propagating states around the boundary at finite Euclidean time.
PHYSICAL REVIEW LETTERS 124, 080501 (2020) 080501-4 Figure 2 displays an exact calculation of masses, for different operator constructions, for a Schwinger model with 6 spatial sites.Each contour shows results obtained using quantum-improved operators of the form Eq. ( 11) optimized on smaller Hilbert spaces than the space they are applied to.In particular, optimizations are undertaken on the 4 spatial-site system, and with different levels of truncation on the total link-square parameter Λ2 .It is clear that even on this smaller system, all operator constructions with truncations Λ2 > 2 perform significantly better than the naive operator presented in Fig. 1.For comparison, note that the dimensionality of the truncated Hilbert space with L ¼ 4 and Λ2 ¼ 4 is 10 in the zero-momentum, odd-parity subspace, compared to 100 for the same subspace of the Λ2 ¼ 12 L ¼ 6 system.It should also be noted that the improved operators from the VQE on truncated spaces approach the effective mass of the GEVP computed on the exact system for the same operator basis.
The potential advantages of this approach are clear.In particular, the scaling of quantum simulations is challenging, and the approach presented here exploits the strengths of both quantum and classical computation; an imprecise interpolating operator extracted from a noisy quantum simulation can still outperform a more naive interpolating operator in a classical calculation, while the classical computation using that operator can be much more easily scaled to large spacetime volumes than the quantum calculation.
Discussion.-Here it is demonstrated how improved interpolating operators for lattice field theory calculations can be constructed using information from VQE computations on NISQ-era quantum hardware.A key feature of this method is that it does not require the large quantum systems that would be needed for a direct calculation of the physics of interest on the quantum machine, but can still accelerate the classical calculation analogously to a classical variational approach, with a potentially significant reduction in classical computing resources.Moreover, while complicated properties of a ground-state system are challenging to access using quantum hardware and have not yet been directly computed on quantum hardware even for the simple Schwinger model system, an optimal interpolating operator obtained via the procedure proposed here can be used in classical computations to accelerate the evaluation of many other properties of the state, including its interactions with external probes.
Extensions of this approach to field theories of phenomenological interest such as QCD could proceed via methods similar to those investigated in Refs.[26,27].Ultimately, there is tremendous opportunity for NISQ-era quantum devices to improve classical field theory calculations.The approach described here represents a clear step towards realizing this goal.

FIG. 1 .
FIG.1.Effective mass functions [Eq.(6)] constructed exactly via diagonalization (solid lines), and obtained from MC computations (open points), for the lightest vector meson.The green circles and red squares denote results obtained using the benchmark [Eq.(7)] and quantum-improved [Eq.(11)] operators, respectively.The black dashed line shows the result obtained solving an exact form of the classical GEVP using the same basis of interpolating operators, while the purple solid band shows the exact result for the optimized operator, given some uncertainty on the VQE input into its construction as described in the text.Note that the uncertainty is suppressed as the Euclidean separation between operators is increased.The large time behavior is a consequence of the backward propagating states around the boundary at finite Euclidean time.

FIG. 2 .
FIG. 2. Exact effective mass functions in an L ¼ 6 spatial-site Schwinger model computation, calculated using quantumimproved interpolating operators constructed on smaller Hilbert spaces.Green, orange, and purple curves (moving vertically down the figure) show results obtained via optimization of systems with flux truncations Λ2 ¼ f2; 4; 8g on an L ¼ 4 system, respectively.The grey dashed and dotted curves indicate the exact classical GEVP solved on the full Hilbert space, and the exact result corresponding to the naive interpolating operator [Eq.(7)], respectively.