Accelerating lattice quantum field theory calculations via interpolator optimization using NISQ-era quantum computing

TL;DR

This paper introduces a hybrid quantum-classical method that uses NISQ-era quantum computations to optimize interpolators, significantly accelerating lattice quantum field theory calculations, demonstrated on the Schwinger model.

Contribution

It presents a novel approach combining small-scale quantum computations with classical calculations to improve efficiency in non-perturbative quantum field theory simulations.

Findings

Quantum-assisted interpolator optimization accelerates calculations.

Method successfully applied to the 1+1-dimensional Schwinger model.

Potential to extend to more complex quantum field theories.

Abstract

The only known way to study quantum field theories in non-perturbative 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.

Equations

$$\mathcal{L}$$

$${\varphi }_n=\{\begin{array}{l}\sqrt{a}{\psi }_{e^{-}}(x_n),\text{ n even,}\\ \sqrt{a}{\psi }_{e^{+}}(x_n),\text{ n odd}.\end{array}$$

$$\hat{U}(x_n,x_{n+1})=e^{iag{\hat{A}}_1(x_n)}\equiv e^{i{\hat{\theta }}_n}.$$

$${\hat{\ell }}_n=-\frac{1}{g}\frac{d{\hat{A}}_1}{dt},$$

$${\hat{\ell }}_n|{\ell }_n⟩={\ell }_n|{\ell }_n⟩,{\ell }_n\in \mathbb{Z}\forall n.$$

$${\hat{L}}_n^{\pm }\equiv e^{\pm i{\hat{\theta }}_n},{\hat{L}}_n^{\pm }|{\ell }_n⟩=|{\ell }_n\pm 1⟩.$$

$${\hat{H}}_{\mathrm{lat}}=$$

$$-iw({\hat{\varphi }}_{N-1}^{†}{\hat{L}}_{N-1}^{+}{\hat{\varphi }}_0-h.c.)$$

$$+m\sum _{n=0}^{N-1}(-1{)}^n{\hat{\varphi }}_n^{†}{\hat{\varphi }}_n+J\sum _{n=0}^{N-1}{\hat{\ell }}_n^2.$$

$$G(\tau )=\sum _x\left[⟨\hat{\mathcal{O}}(x,\tau ){\hat{\mathcal{O}}}^{†}(0,0))⟩-⟨\hat{\mathcal{O}}(x,\tau )⟩⟨{\hat{\mathcal{O}}}^{†}(0,0)⟩\right].$$

$$G(\tau )=\sum _n{\left|⟨n|\mathcal{O}|\Omega ⟩\right|}^2{e}^{-(E_n-E_{\Omega })\tau }.$$

$$M_{\mathrm{eff}}(\tau )=\frac{1}{a}\log \left(\frac{G(\tau )}{G(\tau +a)}\right)\underset{\tau \to \infty }{\to }(E_{\mathcal{O}}-E_{\Omega }).$$

$${\hat{O}}_V(x_n,{\tau }_j)$$

$$|G⟩\approx |G,\vec{\theta }⟩=U^{(m)}({\vec{\theta }}_m)U^{(m-1)}({\vec{\theta }}_{m-1})\dots U^{(1)}({\vec{\theta }}_1)|0⟩.$$

$${\ell }_i^2\le {\Lambda }^2,\forall i=0,\dots ,N-1,\sum _{i=0}^{N-1}{\ell }_i^2\le {\tilde{\Lambda }}^2,$$

$$|{\tilde{V}}^{-}⟩=\sum _k{\alpha }_k{\hat{O}}_k|\Omega ,\vec{\theta }⟩.$$

$${\hat{O}}_{\mathrm{VQE}}=\sum _k{\alpha }_k{\hat{O}}_k,$$

$$({\hat{O}}_1{)}_{\mathrm{src}}(x_i)$$

$$({\hat{O}}_2{)}_{\mathrm{src}}(x_i)$$

$$({\hat{O}}_3{)}_{\mathrm{src}}(x_i)$$

$$({\hat{O}}_4{)}_{\mathrm{src}}(x_i)$$

$$({\hat{O}}_5{)}_{\mathrm{src}}(x_i)$$

$$({\hat{O}}_6{)}_{\mathrm{src}}(x_i)$$

Full text

††preprint:

MIT-CTP/5138

ADP-19-17/T1097

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), can not 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 Wilson (1974); Gattringer and Lang (2010). 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 Aoki et al. (2019); Lehner et al. (2019). 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 Joó et al. (2019); Carlson et al. (2018). 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 Preskill (2018).

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 wavefunction 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+1-dimensional 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 an significant acceleration of LQFT computations for nuclear physics.

The Schwinger model: The Schwinger model Schwinger (1962), 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, including confinement. 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 $\psi$, with mass $m$, coupled via charge $g$ to an electromagnetic field $A_{\mu}$, with the Lagrangian defined as

[TABLE]

where the covariant derivative $\not{D}$ and field strength tensor $F_{\mu\nu}$ are defined in terms of $A_{\mu}$. A discretized formulation of the 1+1D Schwinger model can be defined on a staggered space-time lattice via the Kogut-Susskind prescription Kogut and Susskind (1975); the staggered fermion field $\phi(x_{n})=\phi_{n}$ is defined with the two components of the fermion field occupying even and odd sites by

[TABLE]

In temporal gauge ($A_{0}=0$), the remaining spatial component of the gauge field is encoded on links connecting adjacent staggered sites $x_{n}$ and $x_{n+1}$ by

[TABLE]

Electric flux operators can be defined in terms of $\hat{A}_{1}$ as

[TABLE]

and can be interpreted as acting on the Fock space of links connecting sites: with $\ell_{n}$ denoting the value of the electric flux at the link connecting sites $n$ and $n+1$,

[TABLE]

Ladder operators in this space can be defined as

[TABLE]

Combining the link space with fermionic occupation numbers, a complete Fock space of states in this theory can be expressed as \{\big{|}\vec{n},\,\vec{\ell}\big{\rangle}\}. On a lattice with $N$ staggered sites ($N/2$ spatial sites), and with coupling constants $w=\frac{1}{2a}$ and $J=\frac{g^{2}a}{2}$, the Hamiltonian of this theory can be expressed in terms of these operators as

[TABLE]

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 even-parity scalar $e^{+}e^{-}$ ‘meson’.

Classical computations of ground-state energies: Using classical computation, energy levels of the Schwinger model can be obtained using standard Monte-Carlo (MC) methods. Here a local Hamiltonian MC method is studied Hirsch et al. (1982); details of the application of this approach to the 1+1D Schwinger model are given in Ref. Schiller and Ranft (1983). In this formalism, ground-state energy levels can be determined by the analysis of correlation functions $G(\tau)$, defined in terms of the expectation values of interpolating operators $\hat{\mathcal{O}}(x,\tau)$, which are constructed to create/annihilate states with the quantum numbers of a target state of interest at some Euclidean position $(x,\tau)$:

[TABLE]

Here, a state is created at some initial spatial position and time $(x=0,t=0)$, and annihilated $\tau$ 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 $\ket{\Omega}$ and the tower of excitations coupled to the vacuum through $\hat{\mathcal{O}}$:

[TABLE]

Numerically, the energy gap to the lowest state of interest is determined from the asymptotic value of the effective mass function:

[TABLE]

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 by following the same prescription as for odd-parity meson states in staggered lattice formulations of QCD Golterman (1986):

[TABLE]

Physically, the operator $\hat{O}_{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 generalised eigenvalue problem (GEVP) to achieve an optimized ground-state energy extraction in that sector Michael (1985); Luscher and Wolff (1990); Blossier et al. (2009). This approach has had tremendous success, particularly in the area of meson spectroscopy studies in LQFT Wilson et al. (2015). Nevertheless, the classical variational method is computationally expensive, scaling quadratically with the size of the basis, and in particular using a large variational basis requires high-statistics numerical calculations to ensure a non-singular covariance matrix. This method thus remains intractable for many studies, such as state-of-the-art QCD calculations of nuclear systems Joó et al. (2019). 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: Several years 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 Hauke et al. (2013); Muschik et al. (2017); Martinez et al. (2016); Kokail et al. (2019). Recently, the two-spatial-site Schwinger model has also been studied on IBM’s superconducting quantum hardware111These quantum calculations use an equivalent formulation of the Schwinger model based on bosonic degrees of freedom, related by a Jordan-Wigner transformation Ortiz et al. (2001) to the formulation described here. Klco et al. (2018). In that work, the ground state energy level of the theory in was extracted using the variational quantum eigensolver (VQE) method Peruzzo et al. (2014). In the VQE approach, a sequence of unitary operators $U^{(i)}(\vec{\theta})$, implemented as a sequence of one- and two-qubit gates, are tuned using variational parameters $\theta$ to transform an initial, easy-to-prepare state $\ket{0}$ into an approximation of the ground state of the system in a given symmetry sector, $\ket{G}$:

[TABLE]

From this construction, an approximate value of the ground-state energy of the system can be calculated. In the study of Klco et al. Klco et al. (2018), 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

[TABLE]

for some choice of truncation $\{\Lambda^{2},\tilde{\Lambda}^{2}\}$; 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 wavefunction 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 wavefunctions 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 significantly reduced basis of operators to compute classically on the full state space.

As an explicit example, for the Schwinger model one might aim to study the first excited state of the theory (the odd-parity massive photon). As outlined in Eq. (12), VQE computations can provide approximate descriptions of both this state, $|{V^{-},\vec{\theta}^{\prime}}\rangle$, and the vacuum of the theory $|{\Omega,\vec{\theta}}\rangle$. Acting with a set of interpolating operators $\{\hat{O}_{k}\}$ on the ground state provides a variational basis that can be used to approximate the massive photon:

[TABLE]

The operators $\hat{O}_{k}$ can be defined, for example, in terms of field and link operators $\hat{\phi}_{n}$, $\hat{\ell}_{n}$, $\hat{L}_{n}^{\pm}$, and should be constructed to transition between the vacuum and the symmetry sector of interest. With an estimate for the matrix elements $\langle V^{-},\vec{\theta}^{\prime}|\hat{O}_{k}|\Omega,\vec{\theta}\rangle$, a classical computer can be used to optimize the overlap $|\langle V^{-},\vec{\theta}|\widetilde{V}^{-}\rangle|$ (with appropriate normalization) with respect to $\vec{\alpha}$. This optimization defines an improved interpolating operator

[TABLE]

which has suppressed overlap onto excited states and can be implemented in a classical Euclidean MC calculation. Ultimately, the matrix elements $\langle V^{-},\vec{\theta}^{\prime}|\hat{O}_{k}|\Omega,\vec{\theta}\rangle$ needed for this approach will be evaluated on quantum hardware. Within the present formalism, the Hamiltonian has been block-diagonalized into distinct symmetry sectors, whereas the operators of interest transition between these sectors. Nevertheless, in a space concatenating the two sectors, the operators can be written as products of Pauli operators such that projective measurements can be made in an approach similar to that in Refs. Reiner et al. (2019); Kokail et al. (2019); Peruzzo et al. (2014); Ortiz et al. (2001). In the near term, these matrix elements can also be calculated directly on a classical computer given the quantum gate definitions and the angles {$\vec{\theta}$,$\vec{\theta}^{\prime}$} which encode the approximate wavefunctions; for example, given the 2-layer depth of the variational circuits implemented for the Schwinger model on NISQ devices in Ref. Klco et al. (2018), this can be undertaken in $\mathcal{O}(N^{2})$ operations for an $N$-dimensional VQE basis, which is still computationally tractable on intermediate system sizes where exact diagonalization, at $\mathcal{O}(N^{3})$, is not.

Implementation for the 1+1D Schwinger model: The proposed approach to interpolating operator construction is implemented in the context of the 1+1D Schwinger model defined by Eq. (1). The parameters of the Hamiltonian are chosen to match those in recent studies of this model on quantum devices Klco et al. (2018): $J/\omega=5/3$ and $m/\omega=1/6$, and the temporal lattice spacing is set as $\Delta t=a/5$ based on previous Monte-Carlo studies Schiller and Ranft (1983).

For direct comparison, the mass of the lightest vector meson in this model can be computed exactly and extracted using the local Hamiltonian MC method Hirsch et al. (1982); Schiller and Ranft (1983). 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. (11) 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. (Accelerating lattice quantum field theory calculations via interpolator optimization using NISQ-era quantum computing). This diagonalization is performed with respect to a truncation of the electric fluxes, $\{\Lambda^{2},\tilde{\Lambda}^{2}\}=\{1,8\}$, 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 $\hat{\ell}$ (Eq. (5)), detailed in the supplementary material, is then optimized to maximise the overlap with the target state. This improved operator, of the form Eq. (15), is then used to construct a correlator (Eq. (8)) 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 diagonalisation. The effective mass obtained from a MC computation with the quantum-improved 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 wavefunctions will be only approximately determined, with statistical uncertainties on the variational parameters {$\vec{\theta}$,$\vec{\theta}^{\prime}$}, and the corresponding matrix elements $\langle V^{-},\vec{\theta}^{\prime}|\hat{O}_{k}|\Omega,\vec{\theta}\rangle$ will be similarly limited in measurement precision. The effect of such uncertainties on the definition of the quantum-improved interpolating operator $\hat{O}_{\rm 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 Kokail et al. (2019). Fig. 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 in comparison 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.

Figure 2 displays an exact calculation of effective 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. (15) 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 $\tilde{\Lambda}^{2}$. It is clear that even on this smaller system, all operator constructions with truncations $\tilde{\Lambda}^{2}>2$ perform significantly better than the naive operator presented in Fig. 1. For comparison, we note that the dimensionality of the truncated Hilbert space with $L=4$ and $\tilde{\Lambda}^{2}=4$ is 10 in the zero-momentum, odd-parity subspace, compared to 100 for the same subspace of the $\tilde{\Lambda}^{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 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 space-time 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 is 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. Brower et al. (1999); Wiese (2014). 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.

Acknowledgements

We thank W. Detmold, D. B. Kaplan, B. Kulchytskyy, C. Muschik, M. J. Savage and J. M. Zanotti for helpful discussions. This work is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-SC0011090. PES is supported by the National Science Foundation under CAREER Award 1841699, and PES and AA are supported in part by Perimeter Institute for Theoretical Physics. Research at the Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research and Innovation. RDY is supported in part by the Australian Research Council, grant DP19010029.

I Supplementary Material

This supplementary material makes explicit the variational basis of interpolating operators used for the demonstration of the quantum-improved interpolating operator optimization described in the main text. A basis of 6 operators is chosen, constructed to have the appropriate parity:

[TABLE]

Note that the operator in Eq. 16 may act non-trivially on configurations with one or more electron-positron pairs; operators in Eq. (17)–(18), two or more; operators in Eq. (19)–(21), three or more. For the more harsh truncations in $\tilde{\Lambda}^{2}$, some of these operators annihilate the entire space and hence the corresponding optimisation (Eq. (15)) is performed with respect to a smaller operator basis.