General Methods for Digital Quantum Simulation of Gauge Theories

TL;DR

This paper introduces a general framework for simulating gauge theories with matter fields on digital quantum computers, ensuring gauge invariance and demonstrating the approach on small lattice models including non-Abelian theories.

Contribution

It develops a gauge-symmetry-preserving Trotterized evolution scheme and provides a method to obtain gauge-invariant operators for digital quantum simulation.

Findings

Successfully simulated a 2+1D non-Abelian gauge theory on small lattices.

Presented a gauge-invariant Trotterized time-evolution operator.

Outlined procedures applicable to a broad class of gauge theories.

Abstract

A general scheme is presented for simulating gauge theories, with matter fields, on a digital quantum computer. A Trotterized time-evolution operator that respects gauge symmetry is constructed, and a procedure for obtaining time-separated, gauge-invariant operators is detailed. We demonstrate the procedure on small lattices, including the simulation of a 2+1D non-Abelian gauge theory.

Figures

Tables

Table 1. Table 1: Gate requirements for the propagation of a lattice with N P subscript 𝑁 𝑃 N_{P} plaquettes and L 𝐿 L links, for a time T 𝑇 T with time-steps of size Δ ​ t Δ 𝑡 \Delta t

Gate	Number
${𝔘}_F$	$2L\frac{T}{\Delta t}$
${𝔘}_{\text{phase}}$	$L\frac{T}{\Delta t}$
${𝔘}_{-1}$	$6N_P\frac{T}{\Delta t}$
${𝔘}_{\times }$	$6N_P\frac{T}{\Delta t}$
${𝔘}_{\mathrm{Tr}}$	$N_P\frac{T}{\Delta t}$

Equations

$${\mathfrak{U}}_{-1}|g⟩=|g^{-1}⟩\text{.}$$

$${\mathfrak{U}}_{\times }|g⟩|h⟩=|g⟩|gh⟩\text{.}$$

$${\mathfrak{U}}_{\mathrm{Tr}}(\theta )|g⟩=e^{i\theta \mathrm{Re}\mathrm{Tr}g}|g⟩.$$

$${\mathfrak{U}}_F\sum _{g\in G}f(g)|g⟩=\sum _{\rho \in \hat{G}}\hat{f}(\rho {)}_{ij}|\rho ,i,j⟩$$

$${\mathfrak{U}}_{\times ,f}|g⟩|\varphi ⟩=|g⟩|g\varphi ⟩\text{.}$$

$$⟨{\tilde{\varphi }}_1{\tilde{\varphi }}_2|{\mathfrak{U}}_{⟨\cdot ,\cdot ⟩}(\theta )|{\varphi }_1{\varphi }_2⟩={\delta }_{{\varphi }_1}^{{\tilde{\varphi }}_1}{\delta }_{{\varphi }_2}^{{\tilde{\varphi }}_2}{e}^{i\theta \left[{\varphi }_2^{†}{\varphi }_1+{\varphi }_1^{†}{\varphi }_2\right]}\text{.}$$

$$\varphi (V)|U_{ij}\dots ⟩=|V_j{U}_{ij}{V}_i^{†}\dots ⟩$$

$$P|U_{12}\dots ⟩$$

$$=\frac{1}{|G{|}^N}{\int }_{G^N}\mathrm{d}V\varphi (V)|U⟩$$

$$H=$$

$$=$$

$$\mathcal{U}(t)=\prod _{t/\Delta t}{e}^{-i{H}_K\Delta t}{e}^{-i{H}_V\Delta t}\text{.}$$

$$e^{-i{H}_K\Delta t}=\prod _{⟨ij⟩}{\mathcal{U}}_K^{(1)}(i,j)\text{ and }{e}^{-i{H}_V\Delta t}=\prod _p{\mathcal{U}}_V^{(1)}(p)\text{,}$$

$${\mathcal{U}}_K^{(1)}(i,j)={\mathfrak{U}}_F{\mathfrak{U}}_{\mathrm{phase}}{\mathfrak{U}}_F^{†}.$$

$$E(t)=⟨\Omega |{\tilde{\mathcal{U}}}^{†}(t)\hat{E}\tilde{\mathcal{U}}(t)|\Omega ⟩$$

$$T=T_V^{1/2}{T}_K{T}_V^{1/2}$$

$$⟨{\tilde{U}}_{12}\dots |T_V|U_{12}\dots ⟩={\delta }_{U_{12}\cdots }^{{\tilde{U}}_{12}\cdots }\exp \left(\frac{a_0}{a{g}^2}\sum _x{W}_{\mu \nu }(x)\right)\text{.}$$

$$W_{\mu \nu }(x)=\mathrm{Re}\mathrm{Tr}[U_{x,x+\hat{\nu }}^{†}{U}_{x+\hat{\nu },x+\hat{\mu }+\hat{\nu }}^{†}{U}_{x+\hat{\mu },x+\hat{\mu }+\hat{\nu }}{U}_{x,x+\hat{\mu }}],$$

$$⟨{\tilde{U}}_{12}\dots |T_K|U_{12}\dots ⟩=\prod _{⟨ij⟩}{e}^{\frac{a}{g^2{a}_0}\mathrm{Re}\mathrm{Tr}\left[{\tilde{U}}_{ij}^{†}{U}_{ij}\right]}$$

$$Z$$

$$=\left({\int }_G\mathrm{d}{U}_{12}\cdots \right)\exp \left[\frac{1}{g^2}\sum \mathrm{Re}\mathrm{Tr}\prod U_{ij}\right]$$

$$H=H_V$$

$$\text{where }{H}_V=-\log T_V$$

$$U_C|\Psi ⟩|0⟩=|\Psi ⟩|0⟩\text{ and }{U}_C|\Psi ⟩|1⟩=U|\Psi ⟩|1⟩\text{.}$$

$$\left(⟨\Psi |\otimes ⟨+|\right)U_C^{†}\left(\text{mathbbm}1\otimes {\sigma }_x\right)U_C\left(|\Psi ⟩\otimes |+⟩\right)=\mathrm{Re}⟨\Psi |U|\Psi ⟩$$

$$⟨\Psi |\mathcal{U}(-t)W_{{\mu }^{\prime }{\nu }^{\prime }}(x^{\prime })\mathcal{U}(t)W_{\mu \nu }(x)|\Psi ⟩\text{.}$$

$$H_{{ϵ}_1,{ϵ}_2}(\tau )=H_0+{ϵ}_2\delta (\tau -t)W_{{\mu }^{\prime }{\nu }^{\prime }}(x^{\prime })+{ϵ}_1\delta (\tau )W_{\mu \nu }(x)$$

$$-{\frac{{\partial }^{2}C({ϵ}_1,{ϵ}_2)}{\partial {ϵ}_1\partial {ϵ}_2}|}_{{ϵ}_1={ϵ}_2=0}=⟨\mathcal{U}(-t)W_{{\mu }^{\prime }{\nu }^{\prime }}(x^{\prime })\mathcal{U}(t)W_{\mu \nu }(x)⟩$$

$$H_{ϵ}^{ab}(\tau )=H_0$$

$$+{\widetilde{ϵ}}_2\delta (\tau -t)\mathrm{Im}[U_{ij}^{†}{]}^{ba}+{\widetilde{ϵ}}_1\delta (\tau )\mathrm{Im}[U_{ij}{]}^{ab}\text{,}$$

$$-{\left[\frac{{\partial }^2}{\partial {ϵ}_1\partial {ϵ}_2}-\frac{{\partial }^2}{\partial {\widetilde{ϵ}}_1\partial {\widetilde{ϵ}}_2}\right]}_{ϵ=0}$$

$$\mathrm{Re}⟨\mathcal{U}(-t)[U_{ij}^{†}{]}^{ba}\mathcal{U}(t)[U_{ij}{]}^{ab}⟩\text{.}$$

$$|U_{12}\cdots ⟩\otimes |{\varphi }_1\cdots ⟩\mapsto |(V_2{U}_{12}{V}_1^{†})\cdots ⟩\otimes |(V_1{\varphi }_1)\cdots ⟩\text{.}$$

Full text

I Introduction

Quantum simulations are motivated by inherent obstacles to the classical, nonperturbative simulation of quantum field theories Feynman (1982). Deterministic methods struggle with exponential state spaces. Sign problems stymie Monte Carlo methods both at finite fermion density Troyer and Wiese (2005); Alexandru et al. (2018) and in real-time evolution Alexandru et al. (2016). While large-scale quantum computers will greatly enhance calculations in quantum field theory, for the foreseeable future quantum computers are limited to tens or hundreds of non-error-corrected qubits with circuit depths fewer than a 1000 gates — the so-called Noisy Intermediate-Scale Quantum (NISQ) era. Along with hardware development, theoretical issues impede full use of quantum computers. Despite these limits, algorithms have been demonstrated in toy field theories Dumitrescu et al. (2018); Roggero and Carlson (2018); Lu et al. (2018); Martinez et al. (2016); Kokail et al. (2018); Klco et al. (2018); Lamm and Lawrence (2018); Macridin et al. (2018). Viable quantum simulations require addressing four interconnected issues: representation, preparation, evaluation, and propagation. Proposals in the literature address one or more of these topics. For gauge theories, additional complications arise which we discuss here.

The first hurdle is the representation of a gauge field in a set of qubits, or quantum register. While natural matchings exist for fermionic fields Jordan and Wigner (1928); Verstraete and Cirac (2005); Zohar and Cirac (2018); Whitfield et al. (2016), the digitization of a bosonic field is nontrivial. This is reminiscent of the early days of classical lattice field theory where memory resources limited calculations. Ideas include approximating by: finite subsets Hackett et al. (2018); Alexandru et al. (2019), Fock-state truncation Macridin et al. (2018); Yeter-Aydeniz et al. (2018); Klco and Savage (2018), dual variables Bazavov et al. (2015); Zhang et al. (2018); Unmuth-Yockey (2018); Unmuth-Yockey et al. (2018); Kaplan and Stryker (2018); Meurice (2019), optical representations Zache et al. (2018), and the prepotential formulation Raychowdhury and Stryker (2018). A common suggestion is to limit the register to physical states by gauge fixing, at the cost of increased circuit depth in the time-evolution, although a practical method for nonabelian theories remains undescribed.

Given a register, the next concern is the preparation of desired quantum states. In strongly-coupled theories, how the asymptotic states depend upon the fundamental fields is unknown. Most methods focus on ground state construction. Examples include using quantum variational methods Peruzzo et al. (2014); Kokail et al. (2018), Quantum Phase Estimation (QPE) Abrams and Lloyd (1999); Nielsen and Chuang (2000); Wiebe and Granade (2016), Quantum Adiabatic Algorithm (QAA) Farhi et al. (2000, 2001), and spectral combing Kaplan et al. (2017). Some preparation methods require efficient real-time evolution. Classical-quantum hybrid methods Lamm and Lawrence (2018) and dimensional reduction Bilgin and Boixo (2010) have been proposed to initializing thermal states. Simulations of scattering are hindered by state preparation Jordan et al. (2012, 2011); García-Álvarez et al. (2015); Jordan et al. (2014, 2017); Hamed Moosavian and Jordan (2018).

Evaluation of observables represented by a single, hermitian operator at a single time is straight forward. This includes most static properties. Other observables, like time-separated composite operators (e.g. parton distribution functions) are more complicated. Naively, the first nonunitary operator collapses the quantum state. One resolution introduces ancillary probe and control qubits Pedernales et al. (2014), which we use here. QPE Abrams and Lloyd (1999) has been used to compute linear response Roggero and Carlson (2018). Another intriguing idea uses quantum sensors to implement generating functionals Bermudez et al. (2017). A second problem arises from state contamination. In Euclidean lattice field theory, we can use imaginary-time evolution to separate states overlapping with the same operator. Real-time evolution lacks this separation, so we need novel methods to project onto desired states, perhaps through a distillation-like technique Peardon et al. (2009).

In this work, we construct a method for mapping the Hilbert space of a gauge theory onto the Hilbert space of a quantum computer, and simulating its propagation. Propagating a system for a time $t$ requires a unitary operator $\mathcal{U}(t)=e^{-iHt}$, which cannot be efficiently constructed on a quantum computer. Instead, one Trotterizes the evolution by $\mathcal{U}(t)\approx(e^{-iH\frac{t}{N}})^{N}$ which allows efficient simulation Jordan et al. (2012, 2011, 2017); García-Álvarez et al. (2015); Jordan et al. (2014); Hamed Moosavian and Jordan (2018); Haah et al. (2018). Trotterized versions of $\mathcal{U}(t)$ exist for electron-phonon systems Macridin et al. (2018), the abelian theories Banerjee et al. (2012); Hauke et al. (2013); Kuno et al. (2015, 2017); Martinez et al. (2016); Marcos et al. (2014); Klco et al. (2018); Zache et al. (2018); Muschik et al. (2017); Bazavov et al. (2015); Zhang et al. (2018); Unmuth-Yockey (2018); Unmuth-Yockey et al. (2018); Gustafson et al. (2019); Zohar et al. (2012), and nonabelian theories Byrnes and Yamamoto (2006); Banerjee et al. (2013); Wiese (2013); Zohar et al. (2013); Tagliacozzo et al. (2013); Zohar et al. (2016, 2017); Bender et al. (2018). Each depends on the register used.

Violations of gauge invariance in the time-evolution are potentially introduced through noisy gates, digitization, and Trotterization. To remedy these, Stryker (2018) proposed using oracles to project into the physical subspace. In this work, we construct a Trotterized time-evolution without gauge-fixing, for arbitrary gauge theories with coupled matter fields. Gauge invariance is exactly preserved after Trotterization. We connect this procedure to the transfer matrix formalism Luscher (1977); Mitrjushkin (2002); Creutz (2000); Caracciolo and Palumbo (2013). The effect of noisy gates on gauge-invariance are hardware-dependent, and we leave this for future work. The error introduced by field digitization can be treated without impacting the algorithm.

A final, overarching concern is renormalization. Work in quantum computing Haegeman et al. (2013); Verstraete and Cirac (2010) and elsewhere Rossi and Testa (2018); Briceño et al. (2017) suggests that the renormalization factors may be cheaply computed classically, allowing for the matching of the bare parameters on the quantum computer to renormalized and physical ones. In the context of the method considered here, a mapping between the Trotterization and a Euclidean action may allow renormalized couplings to be computed efficiently on a classical computer with standard lattice methods. The transfer matrix formalism naturally provides such a mapping. Other quantum-computing based proposals exist as well Cotler et al. (2018a, b).

This paper is structured as follows: in Sec. II we describe the prerequisite registers and fundamental gates for our construction. With these, Sec. III formulates gauge-invariant $\mathcal{U}(t)$ for the gauge fields. Sec. IV links this procedure to the transfer matrix formalism, and outlines how renormalization may be performed classically. Sec. V explains how nontrivial correlators may be computed. Sec. VI and VII extend the formulation to include scalar and spinor fields respectively. Sec. VIII demonstrates of our method in the non-Abelian $D_{4}$ theory and the $\mathbb{Z}_{2}$ theory with staggered fermions. We conclude in Sec. IX.

II Prerequisites

Developers of quantum algorithms assume idealized qubits, and certain unitary gates operating on those qubits (a sufficient universal gateset, for example, is the Hadamard gate, the $\frac{\pi}{8}$ gate ($T$), and the controlled-not gate CNOT). This follows a pattern in classical computing, where certain constructs (e.g., floating-point numbers and arrays) are assumed to specify an algorithm. The implication of this is that any correct implementation of the primitive constructs will do: the correctness of the quantum algorithm does not depend on whether superconducting qubits or trapped ions are used.

In this section, we describe a set of higher-level constructs that allows the simulation of a general nonabelian gauge theory. We outline how they might be implemented in practice, from the lower-level qubit/gateset primitives. Hopefully these specifications are helpful to focus the implementation of those tools that are most critical.

We will simulate a theory with a gauge group $G$. The central construct is the $G$-register, analogous to a classical variable storing an element of $G$. The Hilbert space of an ideal $G$-register is $\mathcal{H}_{G}=\mathbb{C}G$, the complex vector space with one basis element $\left|g\right>$ for every element $g\in G$. (Equivalently, it is the space of square-integrable functions from $G$ to $\mathbb{C}$.) The Hilbert space of multiple $G$-registers is constructed as a tensor product, so a computer with $N$ $G$-registers has Hilbert space $\mathcal{H}_{G}^{\otimes N}$.

For discrete groups, ideal $G$-registers are easy to implement by virtue of their finite number of elements; for instance, an ideal $\mathbb{Z}_{2}$-register requires one qubit. For continuous Lie groups (e.g. $SU(3)$), ideal $G$-registers cannot be made with any finite number of qubits – we must settle for approximate $G$-registers (typically with Hilbert space $\tilde{\mathcal{H}}_{G}=\operatorname{span}\{\left|g\right>\}_{g\in\tilde{G}}$ for some subset $\tilde{G}\in G$) which leads to intrinsic discretization error. The tradeoff between number of qubits and the size of discretization error will be crucial to constructing efficient simulations near-term.

In this work, we consider only the case of ideal $G$-registers. Complications arising from the discretization, which may include the breaking of gauge invariance Stryker (2018) and difficulties with the definition of the kinetic and Fourier gates are not considered here. A proposal for approximate $SU(N)$ registers that does not have these problems may be found in Hackett et al. (2018); Alexandru et al. (2019).

For the Hamiltonians we consider, a set of useful primitive gates defined on the $G$-register is: inversion, matrix multiplication, trace, and the quantum Fourier transform. The inclusion of matter fields requires additional multiplication and Fourier transform gates, and an inner product gate. A sample implementation for these operations is specified in Appendix A for the $D_{4}$-register case.

The inversion gate acts on a single register to produce the register of its inverse group element. This is defined in the fiducial basis by

[TABLE]

The other group operation is matrix multiplication, which requires a two-register gate, defined by

[TABLE]

Here we have $\mathfrak{U}_{\times}$ implementing left multiplication, in the sense that the content of the target register was multiplied on the left. This is all that is required here. In general, the combination of $\mathfrak{U}_{-1}$ and left-acting $\mathfrak{U}_{\times}$ permits right multiplication: $\mathfrak{U}_{\times,R}(1,2)=\mathfrak{U}_{-1}(1)\mathfrak{U}_{-1}(2)\mathfrak{U}_{\times}(2,1)\mathfrak{U}_{-1}(2)\mathfrak{U}_{\times}(1,2)$.

Most quantum gatesets — certainly all gatesets proposed for practical use — contain the inverse of every gate as a primitive. This makes it easy, given a circuit implementing one unitary, to obtain a circuit implementing its inverse. Therefore, if $\mathfrak{U}_{\times}^{\dagger}$ is available, $\mathfrak{U}_{-1}$ is readily obtained with the assistance of an ancillary register.

The trace of a plaquette appears in our gauge Hamiltonian, and so to perform this operation we combine with the matrix multiplication gate with a single-register trace gate:

[TABLE]

In our construction, the final gate required on the $G$-register is the Fourier transform gate Coppersmith (2002); Ekert and Jozsa (1996); Jozsa (1998) $\mathfrak{U}_{F}$. This gate acts on a $G$-register to rotate into Fourier space. It is defined by

[TABLE]

The second sum is taken over $\rho$, the representations of $G$, and $\hat{f}$ denotes the Fourier transform of $f$. This gate diagonalizes what will be the ‘kinetic’ part of the Trotterized time-evolution operator. After application of the gate, the register is no longer a $G$-register but a $\hat{G}$-register.

The $G$-register suffices for the representation of the gauge field. In order to represent scalar matter fields, we introduce the $\mathbb{C}^{N}$-register, where $N$ is the dimension of the fundamental representation111For matter in the adjoint representation, $N$ is the dimension of the adjoint representation. of $G$. This register stores a component of $\mathbb{C}^{N}$, and we will define gates allowing transformation under elements of $G$. As with the $G$-register, it is only possible to implement a ‘truncated’ $\mathbb{C}^{N}$ register. (For fermionic matter fields, the $\mathbb{C}^{N}$-register is unneeded, and no truncation is necessary.)

Three gates are needed to make the $\mathbb{C}^{N}$-register useful in a quantum simulation. The fundamental (or adjoint) multiplication gate $\mathfrak{U}_{\times,f}$ is defined by

[TABLE]

The inner product gate, necessary for the field-theory kinetic term (sometimes called the ‘hopping’ term, to distinguish it from the quantum mechanical kinetic term), is defined by

[TABLE]

The inner product is generally a complex number, but this gate extracts only the real part while remaining unitary. Happily, this is all required for simulation of field theory, where only the linear combination $\phi^{\dagger}_{2}\phi_{1}+\phi^{\dagger}_{1}\phi_{2}$ appear.

Lastly, the scalar Fourier transform gate is defined analogously to $\mathfrak{U}_{F}$ above, except for $\mathbb{C}^{N}$: we notate this $\mathfrak{U}^{\mathbb{C}^{N}}_{F}$. The resulting register is a $\mathbb{C}^{N}$-register. Since the quantum Fourier transform on $\mathbb{C}^{N}$ factorizes into $N$ separate transforms on $\mathbb{C}$, the implementation of $\mathfrak{U}^{\mathbb{C}^{N}}_{F}$ consists of $N$ separate calls to $\mathfrak{U}^{\mathbb{C}}_{F}$.

III Pure Gauge

In this section, we construct a time-evolution operator suitable for the simulation of pure gauge theories on a quantum computer. We take care to exactly preserve gauge invariance (up to unavoidable gate errors and noise). We take the standard approach, by discussing the Hamiltonian formulation of the theory, and then Trotterizing the Hamiltonian into a sum of pieces easily implemented via the primitive constructs of Sec. II. This approach will be extended in Sec. VI and VII to include scalar and fermionic matter fields, respectively.

A similar approach to nonabelian gauge simulation has been discussed previously Bender et al. (2018); Zohar (2018); Zohar et al. (2012); Tagliacozzo et al. (2013); Zohar et al. (2017, 2016). In this section, we take care to work only with logical qubits and gates, keeping the algorithm independent from the underlying geometry of the quantum computer. Additionally, we avoid the need for ancillary $G$-registers.

In the present work, we use the basis diagonalizing the gauge field operators $U$, rather than their conjugate momenta, unlike e.g. Byrnes and Yamamoto (2006) in which the fiducial basis is chosen to diagonalize the kinetic term of the Hamiltonian. Our choice simplifies the calculation of plaquettes and Wilson lines, as well as the inclusion of matter fields. Additionally, simulations done in this ‘position’ basis correspond more closely to calculations performed in classical lattice field theory.

We begin by reviewing the Hamiltonian formulation of lattice gauge theory Kogut and Susskind (1975). Each gauge link has a local Hilbert space $\mathbb{C}G$. Combined, the $L$ links on our lattice have Hilbert space $\mathcal{H}=\mathbb{C}G^{\otimes L}$.

On a lattice with $N$ sites, the space of gauge transformations is $G^{N}$, where one element of $G$ must be specified at each site. The transformation rule for a link $U_{ij}$ from site $i$ to site $j$ is $U_{ij}\mapsto V_{j}U_{ij}V^{\dagger}_{i}$. The linear action of $V\in G^{N}$ on the Hilbert space is given by $\phi(V)$:

[TABLE]

Gauge invariance demands that two states differing only by a gauge orbit be considered physically equivalent. The physical Hilbert space222The Gribov ambiguity prevents us from unambiguously assigning a single representative configuration to each gauge orbit; however, we do not take that approach. Physical states in the Hilbert space are not configurations, but formal linear combinations of configurations, and these can be constructed unambiguously., then, is $\mathcal{H}_{P}=\mathbb{C}G^{\otimes L}/\phi(G^{\otimes N})$.

$\mathcal{H}_{P}$ can be obtained from the larger space $\mathcal{H}$ by a gauge-symmetrization operator, which acts as a projection operator:

[TABLE]

This discussion of the gauge-invariance of physical states is equivalent to the imposition of Gauss’s law. An implementation of the operator $P$ for $U(1)$ is given by Stryker (2018); however our method does not require implementing $P$.

We define a Hamiltonian which acts on the entire space $\mathcal{H}$ of the gauge theory, implicitly defining a ‘physical’ Hamiltonian on the subspace $\mathcal{H}_{P}$. As will be made clear in later sections, this choice of Hamiltonian is motivated by a desire to recover the Euclidean Wilson action.

[TABLE]

The first term, $H_{V}$, is a sum over spatial plaquettes $p$ and can be thought of as a potential term. The trace is taken in a suitable representation of the gauge group $G$ — for matrix groups, this is typically the fundamental representation. The second term, $H_{K}$, is the quantum-mechanical kinetic term describing the mechanics of a free particle moving on the surface of $G$. Therefore, $\pi_{ij}^{2}$ is the Laplace-Beltrami operator on the surface $G$, and $\pi_{in}$ is the momentum operator conjugate to $U_{ij}$. This prescription gives us no guidance for how to form the kinetic term of a discrete gauge theory — this is addressed via the transfer matrix formalism in Section IV.

Not only is the Hamiltonian gauge-invariant, but each term is gauge-invariant. This is crucial, as it allows us to Trotterize without breaking gauge invariance at any order in $\Delta t$. A time evolution operator suitable for implementation on a quantum computer is then

[TABLE]

Below, we will use the notation $\mathcal{U}(\Delta t)$ to indicate a single $\Delta t$ step. We now describe how the operators $e^{-iH_{K}\Delta t}$ and $e^{-iH_{V}\Delta t}$ may be obtained. First, each part (kinetic and potential) of the Hamiltonian consists of mutually commuting, local terms, and each such term may be treated individually and sequentially without approximation. That is,

[TABLE]

Here, again because the factors commute, the order of application of the $\mathcal{U}_{V}^{(1)}$ makes no difference to the final result. The plaquette evolution $\mathcal{U}_{V}^{(1)}$ is constructed by first preparing the product of the plaquette in an ancillary register, and then operating on the ancillary with $\mathfrak{U}_{\operatorname{Tr}}(\frac{1}{g^{2}a})$. $\mathcal{U}_{K}^{(1)}$ is implemented by diagonalizing via $\mathfrak{U}_{F}$, and then applying a diagonal unitary:

[TABLE]

Circuits implementing $\mathcal{U}_{K}^{(1)}$ and $\mathcal{U}_{V}^{(1)}$ are shown in Fig. (1). The implementation of $\mathfrak{U}_{\mathrm{phase}}$ depends strongly on the group and the implementation of the $G$-register, but in all cases it is a diagonal operator. The total gate requirements for the propagation are shown in Table 1.

Note that $\mathcal{U}_{V}$ involves multiple $\mathcal{U}_{V}^{(1)}$, some of which affect overlapping links. Because all potential operators commute (indeed all $\mathcal{U}_{V}^{(1)}$ are diagonal in the fiducial basis), the order they are performed in does not matter.

When preparing an initial state, it is important to ensure that it lies within the physical subspace $\mathcal{H}_{P}$. One procedure for projecting onto the gauge-invariant states for an Abelian gauge theory is provided in Stryker (2018).

Although the preparation of the ground state $\left|\Omega\right>$ is beyond the scope of this paper, it is worth discussing that common proposals typically require the time-evolution operator $\mathcal{U}(t)$ as constructed in this section. Assuming such methods are used, nonequilibrium physics can be computed by preparing the ground state with $\mathcal{U}(t)$, and then evolving with a different time-evolution operator $\mathcal{\tilde{U}}(t)$ Lamm and Lawrence (2018). A typical expectation value of interest, accessible in such a scheme, is

[TABLE]

where $\hat{E}$ measures the expectation of a plaquette (i.e. the energy density of the gauge field). Other expectation values of interest are discussed in Sec. V below.

IV The Lattice Path Integral

Substantial insight into QCD has been achieved via Euclidean-spacetime lattice methods on classical computers. Accordingly, the Euclidean lattice path integral occupies a dominant place in the QCD literature and the minds of practitioners. This is in contrast to quantum simulations where the Hamiltonian formulation appears more natural. In this section, we connect the Hamiltonian formalism to the Euclidean lattice path integral via a transfer matrix. This implies an important practical consequence: the determination of appropriate bare coupling constants may be performed on a classical computer. This procedure also guides us in constructing a Hamiltonian for a discrete gauge theory (discussed in Harlow and Ooguri (2018)).

Our formalism for the transfer matrix, $T$, slightly differs from the usual one Luscher (1977); Mitrjushkin (2002); Creutz (2000); Caracciolo and Palumbo (2013); Creutz (1977) in that $T$ is defined on the entire space $\mathcal{H}$. The usual transfer matrix is defined only on the physical Hilbert space, and may be obtained by projection with $P$. In the course of performing this projection, the lattice path integral with the Wilson action is recovered.

The transfer matrix is an approximation to imaginary-time evolution by one unit of time: $T\approx e^{-H}$. We construct the transfer matrix in the fiducial basis of $\mathcal{H}$ from separate kinetic $T_{K}$ and potential $T_{V}$ contributions via the split-operator approximation

[TABLE]

where the potential term resembles the product of spatial plaquettes as appears in the Wilson action

[TABLE]

Here $a$ is the spatial lattice spacing, and $a_{0}$ the temporal spacing; we do not assume an isotropic lattice. We have borrowed from lattice field theory the Wilson plaquette

[TABLE]

and $\mu,\nu$ are restricted to space-like directions.

For the kinetic term there is no coupling between different links, but it is related in lattice field theory to plaquettes containing timelike links:

[TABLE]

Note that, like the Trotterized real-time evolution, $[T,P]=0$ due to the fact that $T_{V}$ and $T_{K}$ individually commute with $P$.

When the transfer matrix is restricted to physical states, it exactly reproduces the usual lattice path integral, with the Wilson action.

[TABLE]

Here $\operatorname{Tr}_{P}$ denotes the trace over only the physical subspace $\mathcal{H}_{P}$. The inverse temperature is denoted $\beta$, and $Z$ is the Euclidean lattice approximation to the partition function.

Now that we have a transfer matrix corresponding to the Wilson action, we may connect it to the Hamiltonian. This is done in detail by Creutz Creutz (1977) for the gauge-fixed transfer matrix, so we will only summarize the gauge-free result (which is essentially the same) here. The potential part $T_{V}$ of the transfer matrix is exactly $e^{-H_{V}}$. The relationship between the kinetic part $H_{K}$ of the Hamiltonian and $T_{K}$ is more subtle. It is shown in Creutz (1977) that, in the limit where the temporal lattice spacing $a_{0}$ is taken to [math], the kinetic part of the Hamiltonian is recovered.

With a mapping between the Hamiltonian and lattice action through the transfer matrix, the bare constants in Eq. (III) for a desired renormalized, physical set of couplings are the same as the lattice action ones. This is particularly important for ensuring that Lorentz invariance is recovered in the continuum limit, but also allows us to perform simulations at known values of physical parameters. In particular, in a lattice QCD simulation, the task of determining appropriate bare parameters for a desired physical pion mass may be performed cheaply on a classical computer, reserving the quantum processor for the difficult real-time evolution.

Clearly, the Hamiltonian evolution on a quantum computer is closely related to standard classical lattice simulations. Indeed, the two differ only in: the discretization of the gauge group (quantum computers have smaller $G$-registers), whether the Trotterization is done in real or imaginary time, and the exact form of the kinetic term. The first difference can be removed by performing a classical simulations with reduced precision, and the others are reduced at small temporal lattice spacing.

In the case of the Trotterization $\Delta t$ being equal to $a$, the quantum evolution is simply an analytic continuation of the usual isotropic lattice calculations performed in lattice QCD. This implies that Lorentz invariance is recovered in the continuum limit.

We conclude by exploiting this procedure to produce a Hamiltonian for a discrete gauge theory — that is, a gauge theory where the gauge group is a discrete group (such as $D_{4}$, demonstrated below), rather than a connected manifold. In this case, the kinetic part of the Hamiltonian is determined not by a Laplacian on a manifold, but can be derived from relevant portion of the transfer matrix by taking a logarithm.

[TABLE]

where $T_{V}$ and $T_{K}$ are defined exactly as in Eq. (14). Note that $H_{K}$ is only well-defined if $T_{K}$ is positive-definite: this fact is established for a general discrete group (with the Wilson action) in Appendix B.

As a consequence of Schur’s lemma, the kinetic transfer gate defined here is diagonal in the Fourier basis. This implies that, where $H_{K}$ is obtained in this fashion, $\mathcal{U}_{K}$ may be easily implemented from $\mathfrak{U}_{F}$.

V Correlators

Having constructed time-evolution for a pure gauge theory, we now must show how observables of physical interest may be obtained. In this section, we show how to obtain a plaquette-plaquette correlator and the expectation value of an arbitrary, temporally extended Wilson loop. Such correlation functions serve as the nonperturbative input to larger calculations, such as the determination of parton distribution functions.

We accomplish both correlators by the technique set out in Pedernales et al. (2014). To begin, we recall how one can compute an expectation value of a unitary operator $U$ Ortiz et al. (2001) in a given state $\left|\Psi\right>$. Introducing a single ancillary qubit, we construct a controlled unitary operator $U_{C}$, defined by

[TABLE]

(Standard quantum gatesets make the construction of this operator straightforward.)

Generally, the expectation value of $U$ has both real and imaginary parts. These must be measured separately. A measurement of $\operatorname{Re}\left<\Psi\right|U\left|\Psi\right>$ is given by initializing the ancillary qubit to $|+\rangle\equiv\frac{1}{\sqrt{2}}\left(\left|0\right>+\left|1\right>\right)$, performing evolution via $U_{C}$, and then measuring $\sigma_{x}$ on the ancillary qubit.

[TABLE]

For $\operatorname{Im}\left<\Psi\right|U\left|\Psi\right>$ one performs the same process except ending with a measurement of $\sigma_{y}$.

With this procedure in mind, we can show how to compute a correlator of the form

[TABLE]

In QCD this Wilson plaquette correlator corresponds to a glueball propagator. The operator in Eq. (22) is clearly not unitary, so it cannot be directly evaluated by means of the procedure described above. Instead, using the Hermiticity of $W_{\mu\nu}(x)$, we construct a parameterized family of unitary operators. The derivatives of this family give the operator of interest. Explicitly, this is done via a time-dependent perturbation of the Hamiltonian:

[TABLE]

Time evolving forward in time with $H_{\epsilon_{1},\epsilon_{2}}$, and back with $H_{0}$, allows us to measure the expectation value $C(\epsilon_{1},\epsilon_{2})\equiv\left<\Psi\right|\mathcal{U}(-t)\mathcal{U}_{\epsilon_{1},\epsilon_{2}}(t)\left|\Psi\right>$. Differentiating twice with respect to the perturbation yields the plaquette-plaquette correlator as desired.

[TABLE]

In practice, this derivative must be obtained by finite differencing, after computing $C(\epsilon_{1},\epsilon_{2})$ for several small values of $\epsilon$. This approach is naturally extended to the computation of $k$-time correlation functions, at the cost of requiring the $k$th numerical derivative.

The terms $W_{\mu\nu}(x)$ added in the perturbation also appear in the original Hamiltonian, and are readily implemented with the primitive gates defined in Section II.

We can use the same method to compute the expectation values of temporally-extended Wilson loops, $\left<\operatorname{Re}\operatorname{Tr}U^{\dagger}_{ij}(t)U_{ij}(0)\right>$. To do so, we decompose it into terms $\left<\operatorname{Re}([U^{\dagger}_{ij}(t)]^{ba}[U^{ab}_{ij}(0)]^{ab})\right>$ where $a,b$ are color indices. We can use this to derive a perturbed Hamiltonian

[TABLE]

which can be utilized in a correlator $C^{ab}(\epsilon_{1},\tilde{\epsilon}_{1},\epsilon_{2},\tilde{\epsilon}_{2})\equiv\left<\Psi\right|\mathcal{U}(-t)\mathcal{U}^{ab}_{\epsilon}(t)\left|\Psi\right>$ We now take the difference of two second derivatives of the resulting correlator:

[TABLE]

Summing over the $a,b$ yields the desired, gauge-invariant trace. The imaginary part can be obtained similarly.

For many gauge groups, the different correlators $\left<(U^{\dagger}_{ij}(t))^{ba}(U_{ij}(0))^{ab}\right>$ are related by gauge symmetry. Consequently, for these gauge groups, the correlator in Eq. (26) need only be evaluated for one particular selection of $a$ and $b$, as all the other terms will be equal. In particular, this is true for $D_{4}$ (simulated below) and $SU(N)$ in the fundamental representation.

The perturbation introduced, unlike the one used for plaquette-plaquette correlators above, is not gauge-invariant, meaning that the state during this time evolution does not lie within the physical subspace $\mathcal{H}_{P}$. Also unlike the plaquette-plaquette perturbation, this perturbation contains terms not present in the original Hamiltonian, and requires gates not specified in Section II. The gates required are diagonal (in the fiducial basis) phase gates on individual $G$-registers.

This method can also compute a non-gauge-invariant expectation value, e.g. $\left<\operatorname{Tr}U_{ik}(t)U_{ij}(0)\right>$ where $j\neq k$. Provided the initial state is correctly prepared, this expectation value must be [math].

VI Scalar fields

Here we extend our Hamiltonian construction to include scalar matter fields which transform in a $D$-dimensional representation of $G$. The Hilbert space defined in Section III is extended via tensor product to include $D$ scalar degrees of freedom at each site. The new Hilbert space is $\mathcal{H}=(\mathcal{H}_{G})^{\otimes L}\otimes(\mathcal{H}_{S})^{\otimes N}$, where $\mathcal{H}_{G}$ is the one-link Hilbert space defined above, $N$ is the number of sites, and $\mathcal{H}_{S}$ is the one-site Hilbert space. For $\mathcal{H}_{S}$, it is convenient to work in the basis of occupation number, writing $\mathcal{H}_{S}=\mathop{\mathrm{span}}\{\left|n_{1}\cdots n_{D}\right>:n_{i}=0,1,\ldots\}$. This is different from classical lattice field theory where scalars are typically represented in a basis diagonalizing the field operator or integrated out, but is more natural for quantum simulations.

As before, we first consider how gauge transformations act on $\mathcal{H}$. A gauge transformation is still specified by $V_{i}\in G$ at each site $i$. The gauge links still transform as $U_{ij}\mapsto V_{j}U_{ij}V_{i}^{\dagger}$, and the scalar fields transform (in the fundamental) via $\phi_{i}\mapsto V_{i}\phi_{i}$. In this way, $\phi_{j}^{\dagger}U_{ij}\phi_{i}$ is gauge-invariant. The action of a gauge transformation $V$ on the Hilbert space is

[TABLE]

The gauge-symmetrization operator is still defined by Eq. (III), but acts on the expanded Hilbert space.

To represent the physical Hilbert space on a quantum computer, we now make use of the $\mathbb{C}^{N}$-registers transforming under $G$. For the scalar fields we consider, the Hamiltonian is

[TABLE]

Although we assume scalar interactions only to order $\phi^{2}$, the generalization to higher-order couplings is trivial and introduces no complications into what follows. We construct a Trotterized time evolution operator, just as in Section III, by decomposing into potential and kinetic terms. The new gates required for the scalar fields are $\mathcal{U}_{M}^{(1)}$ (implementing the mass term), $\mathcal{U}_{H}^{(1)}$ (implementing the hopping term), and $\mathcal{U}^{(1)}_{K,\mathrm{scalar}}$. These are defined by

[TABLE]

Note that the mass and hopping gates commute, as both are diagonal in the fiducial basis. The same-site $\phi_{i}^{\dagger}\phi_{i}$ terms generated by the hopping term of the Hamiltonian have been absorbed into the mass gate, assuming a square lattice of dimension $d$. These gates, together with the pure-gauge $\mathcal{U}_{V,K}^{(1)}$, are combined to form

[TABLE]

With these, the full time-evolution operator for a single time-step is

[TABLE]

A single time step, therefore, consists of five steps on the quantum computer:

For each plaquette, compute the plaquette product and apply $\mathfrak{U}_{\operatorname{Tr}}(\Delta t/g^{2}a)$. This implements $\mathcal{U}_{V,\mathrm{gauge}}$. 2. 2.

For each link, apply the scalar hopping gate $\mathcal{U}_{H}^{(1)}$. 3. 3.

For each site, apply the scalar mass gate $\mathcal{U}_{M}^{(1)}$. 4. 4.

For each link, apply the kinetic gate $\mathcal{U}_{K,\mathrm{gauge}}^{(1)}$. 5. 5.

For each site, apply the scalar kinetic gate $\mathcal{U}_{K,\mathrm{scalar}}^{(1)}$.

This is the general algorithm for time-evolving gauge fields in the presence of matter. Before moving on to discuss scalar field correlators, we note that the transfer matrix discussion of Sec. IV could be repeated here, including the scalar fields, to connect the evolution of the quantum computer to a classical, Euclidean spacetime path integral. Thus, renormalization calculations again may be performed classically.

Steps 1 and 4 above are the same as for a pure-gauge simulation. Step 2 requires two calls to $\mathfrak{U}_{\times,f}$, two calls to $\mathfrak{U}_{-1}$, and one call to $\mathfrak{U}_{(\cdot,\cdot)}$. Step 3 is a register-implementation dependent diagonal phase gate. Step 5, analogous to the implementation of $\mathcal{U}_{K,\mathrm{gauge}}^{(1)}$, entails $2N$ calls to $\mathfrak{U}_{F}^{\mathbb{C}}$ and $N$ diagonal phase gates per site, where $N$ is the dimension of the scalar’s representation.

The procedure for computing plaquette correlators is the same as in Section V.

The construction of Section V may be generalized to obtain the gauge-invariant two-point function of a fermion, connected by a Wilson line: $\left<\bar{\psi}(y,t)W\psi(x,0)\right>$. However, for a general Wilson line, this construction requires one derivative to be taken at every timeslice. As the derivatives must be taken via finite differencing, this is not a plausible method for computing these correlators, even in the far future. Constructing a more efficient method remains an open problem.

VII Fermion fields

Here we consider the inclusion of spinor matter fields, which can be done in the same manner as was discussed above for scalar fields. Also as before, the representation of fermionic fields on the quantum computer will not match that on a classical computer, where they are typically integrated out. In the implementation of spinor fields on a quantum computer, while the register is trivial to specify and exact, a significant technical difference emerges: the fundamental operators provided on most digital quantum computers are bosonic, meaning that operators acting on different qubits will commute with each other. For spinor fields, we want fermionic operators, which anticommute at different physical sites. We will translate commuting operators into anticommuting operators via the Jordan-Wigner transformation, paying the price that local fermionic operators are generated by nonlocal bosonic operators. More efficient, albeit complex, schemes, relating local fermionic operators to more local bosonic operators have been devised Verstraete and Cirac (2005); Zohar and Cirac (2018); Whitfield et al. (2016).

The Hilbert space of a gauge theory with $d$-component fermions can be written $\mathcal{H}=\left(\mathcal{H}_{G}\right)^{\otimes L}\otimes\left(\mathcal{H}_{F}\right)^{\otimes dN}$, where there are $L$ links, $N$ sites, and $\mathcal{H}_{F}=\operatorname{span}\{\left|0\right>,\left|1\right>\}$ is the Hilbert space of a single fermion on a site333As a consequence of writing the fermionic space as a tensor product, the fermionic operators will be nonlocal.. Note that, as a consequence of the exclusion principle, $\mathcal{H}_{F}$ can be put on the quantum computer without approximation or truncation.

Writing the Hamiltonian using operators that obey the standard fermionic anticommutation relations is easily accomplished in the absence of gauge fields via the Jordan-Wigner transformations Jordan and Wigner (1928): each fermionic operator pair is given an index from $1\cdots dN$, and the $i$th operators are constructed from bosonic operators $a_{i}$ and $a^{\dagger}_{i}$:

[TABLE]

In the presence of gauge fields, we must ensure that the fermionic operators obey the appropriate gauge transformation law. To be concrete, consider when the fermions transform in the fundamental of $SU(3)$. Then at each site $i$, for each spinor index $\alpha$, there are three annihilation and creation operators, which must transform under gauge transformations as

[TABLE]

where $a,b$ are color indices. There are now three times as many fermion operators, as for each site and each spinor, there must be three different fermionic degrees of freedom.

This may be accomplished if the bosonic operators $a,a^{\dagger}$ from which the fermions are constructed transform in the fundamental of $SU(3)$. Thus each operator has a color index, and $a^{a}\mapsto g_{ab}a^{b}$. Crucially, no change to the Jordan-Wigner procedure is needed since the $\sigma$-matrices are independent of the gauge group (however the $\sigma$ matrices do inherit the color index).

With fermionic operators constructed, the time-evolution operator is constructed in a similar fashion to that from Section VI above.

[TABLE]

The key difference is that the spinor kinetic and potential operators are derived from the Dirac equation. The governing Hamiltonian is now

[TABLE]

and the single-site and single-link evolution operators are defined as before:

[TABLE]

An implementation of $\mathbb{Z}_{2}$ fermions is demonstrated in Section VIII.

VIII Demonstration

In this section we present two concrete demonstrations of the method described above: first, a two-plaquette theory with the discrete, non-Abelian gauge group $D_{4}$, and second, a $\mathbb{Z}_{2}$ gauge field coupled to fermions. For all our simulations, we work in units where the lattice spacing is $a=1$. We perform classical simulations of a quantum computer using qiskit Santos (2017); Cross et al. (2017) in order to demonstrate the correctness of the implementation. These calculations are performed without modeling of realistic noise sources, corresponding to a perfect, error-free quantum computer. As will be discussed below, the number of entangling gates required for a single Trotterization step is comparable to the total number used in state-of-the-art quantum simulations Nam et al. (2019), and attempting to run on a current quantum computer or using a realistic noise model would be impractical.

We simulate the $D_{4}$ gauge field on the two-dimensional lattice shown in Fig. (2). This lattice represents the smallest two-dimensional lattice which cannot be reduced to a one-dimensional theory. The simulation requires a four $D_{4}$ registers, and uses a total of 14 qubits: 12 for physical degrees of freedom, and 2 ancillary qubits. Note that, for brevity, we have broken with the notation of previous sections, in referring to a link not by the source and sink sites, but instead with a single direct index $0\ldots 3$.

We define a trace on $D_{4}$ (not uniquely specified by the group structure) by embedding $D_{4}$ into $U(2)$, and defining the trace via the fundamental representation of that Lie group. The embedding of $D_{4}<U(2)$ is generated by the elements $\sigma_{x}$ and $i\sigma_{z}$. The lattice action for this model is

[TABLE]

The last term is a sum over all temporal plaquettes on the $2+1$ lattice. The resulting Hamiltonian is

[TABLE]

where $\log T_{K}^{(1)}$ is the one-link kinetic term, determined as discussed in Section IV.

The starting state of our $D_{4}$ simulation is the gauge-invariant projection of the eigenstate of the operators $U_{i}$, where all links have been assigned to the identity matrix. We evolve the system for $t=10$ with two different Trotterization time steps, one of $\Delta t=0.2$ and one of $\Delta t=0.5$. With the larger Trotterization, errors build up more quickly, becoming noticeable at large times. The circuits used are detailed in Appendix A. In total, the quantum simulation of Fig. (3) entailed $\sim 200$ entangling gates per Trotterization time step. This is roughly the resources recently used to compute the ground state energy of the water molecule Nam et al. (2019). Fig. (3) shows the average plaquette energy of the left plaquette, as a function of time. The only systematic error is Trotterization, as there is no noise on the simulated quantum computer.

For the same lattice gauge theory, we also demonstrate the procedure described in Section V for determining the expectation value of a Wilson loop with temporal extent $t$: $\left<\operatorname{Re}\operatorname{Tr}U_{0}^{\dagger}(t)U_{0}(0)\right>$. Simulated results are shown in Fig. (4). For the finite differencing, we perform simulations with $(\epsilon_{1},\epsilon_{2})=(0.1,0.0),(0.0,0.1),(0.1,0.1)$, and the same for the parameters $\tilde{\epsilon}_{1},\tilde{\epsilon}_{2}$. As before, a Trotterization time step of $\Delta t=0.2$ is used. While smaller $\epsilon$ is preferred (smaller deformations improve the finite differencing calculation), there is a tradeoff because the resulting smaller finite difference requires a larger number of samples to be collected, scaling like $\epsilon^{-2}$.

In order to test our technique with matter fields, we simulated staggered fermions transforming in a $\mathbb{Z}_{2}$ gauge field, in $1+1$ dimensional spacetime. A similar model was also simulated on a quantum computer in Klco et al. (2018). The lattice Hamiltonian for this model is

[TABLE]

where the index $i$ denote lattice sites, the sum is taken over all pairs of adjacent sites, and $\chi$ is a one-component fermionic operator.

Using the Jordan-Wigner procedure, we translate from the fermionic operators $\chi_{i}$ into a set of bosonic operators $\sigma_{i}$: $\chi_{i}=\sigma^{(z)}_{1}\cdots\sigma^{(z)}_{i-1}\sigma_{i}$, where $\sigma_{i}$ is the lowering operator for site $i$. In one-dimension without periodic boundary conditions, this results in a local bosonic Hamiltonian.

[TABLE]

There are two distinct families of bosonic operators: the $\sigma_{i}$ which live at a site and correspond to fermions, and the $\sigma_{i(i+1)}$ on the link from site $i$ to site $i+1$, corresponding to the $\mathbb{Z}_{2}$ gauge field. The Trotterization of this Hamiltonian takes four steps, corresponding exactly to the four terms in Eq. (VIII). In this way full translational invariance (neglecting the boundary) is maintained even in the approximated evolution.

The starting state of this $\mathbb{Z}_{2}$ simulation is obtained by projecting the $\left|0000000\right>$ state to the gauge-invariant space, and then applying the gauge-invariant excitation operator $\chi_{1}^{\dagger}\chi_{2}^{\dagger}\sigma^{(z)}_{12}$. The number density of the left-most site as a function of $t$ is shown in Fig. (5) where the simulated results are compared to the exact values. Since we do not begin in an eigenstate of the Hamiltonian, this is not constant. The difference between the two is due to the Trotterization, and as with the $D_{4}$ simulation may be made systematically smaller by reducing the lattice spacing. This simulation uses $8$ qubits: $4$ fermion sites, $3$ gauge links, and a single ancillary. For this lattice, $36$ CNOT gates are used per timestep, putting the simulation (slightly) out of reach of present-day processors. A timestep of $\Delta t=0.6$ is used for $T=30$.

IX Discussion

In this work, we have given a general procedure to construct a unitary time-evolution operator and timelike separated correlation functions for use on an idealized universal quantum computer. This construction ensures that the Trotterization never mixes the physical and unphysical states within the Hilbert space, provided that the requirements of Sec. II are implemented exactly. Further, we show how gauge-invariant correlators (including a temporal Wilson loop) are constructed. We have demonstrated the efficacy of the procedure by the implementation of two gauge theories on a simulated quantum computer. Additionally, we have discussed the classical procedure for determining the bare couplings to be used in a quantum simulation.

Critically, we have defined a set of primitive constructs which may be used in the construction of gauge theory simulation. The practical and efficient implementation of such registers and gates is an important open question in the literature. Finally, future work must address the issue of how the behavior of different noisy gates affects the mixing of the physical and unphysical state, whether this can be mitigated by special constructions or error-mitigating codes within our construction, and indeed whether such mitigation is necessary.

Acknowledgements.

H.L., S.L., and Y.Y. are supported by the U.S. Department of Energy under Contract No. DE-FG02-93ER-40762. The authors thank P. Bedaque and A. Alexandru for conversations on the Hamiltonian formulation of gauge theories, Z. Davoudi for a quantum computing seminar, J. Stryker for helpful comments on this paper. Finally, the authors would like to thank C. Bonati, M. Cardinali, G. Clemente, L. Cosmai, M. D’Elia, A. Gabbana, S. F. Schifano, R. Tripiccione, and D. Vadacchino for pointing out errors in our manuscript.

Appendix A Implementation of a $D_{4}$-register

The group $D_{4}$ of isometries of the square is treated as a subgroup of $U(2)$ generated by the matrices

[TABLE]

The first matrix represents a $\frac{\pi}{2}$ rotation, and the second represents reflection.

The $D_{4}$ register is implemented with three qubits. The state $\left|abc\right>$ corresponds to the matrix

[TABLE]

Thus, the two least-significant bits specify an amount of rotation, to be followed by a flip if the most-significant bit is $1$. This establishes the isomorphism needed between $\mathbb{C}D_{4}$ and the three-qubit Hilbert space on a quantum computer.

It remains to describe the inversion, multiplication, trace, and Fourier transform circuits. The inversion and multiplication circuits are both equivalent to classical circuits (as they consist exclusively of controlled-not gates), and are shown in Fig. (6). The trace circuit is a three-qubit controlled phase gate, defined by $\mathfrak{U}_{\operatorname{Tr}}(\theta)\left|000\right>=e^{-2i\theta}$ and $\mathfrak{U}_{\operatorname{Tr}}(\theta)\left|010\right>=e^{2i\theta}$ (with all other states picking up no phase).

Finally, the Quantum Fourier Transform which diagonalizes the $D_{4}$ momentum operator is given by

[TABLE]

Rather than attempt to determine a quantum circuit by hand (although general methods do exist Moore et al. (2003); Pueschel et al. (1998); Beals (1997)), we perform a classical simulated annealing search to find a short circuit that implements exactly this unitary within the qiskit gate set using Hadamard, CNOT, and $\frac{\pi}{8}$-gates ($T$). The circuit found is shown in Fig. (7). The annealing algorithm will be detailed in a future work. Since the size of $F$ does not scale with volume, this method does not suffer from an exponentially large physical Hilbert space.

After being diagonalized by the Fourier transform, the kinetic gate corresponds simply to a phase gate on the most-significant qubit, along with a phase gate on the state $\left|000\right>$:

[TABLE]

The constants $\theta_{1}$ and $\theta_{2}$ are coupling-dependent. For $\beta=0.8$, the coupling used in this paper, we have $\theta_{1}\approx 1.263$ and $\theta_{2}\approx 0.409$.

Appendix B Positive-Definiteness of the Transfer Matrix

Here we provide a proof that the kinetic part $T_{K}$ of the transfer matrix for a discrete gauge theory is positive-definite to define a Hamiltonian, a fact used in Sec. IV. Unlike prior proofs, we are concerned with positive-definiteness on $\mathcal{H}$, not just $\mathcal{H}_{P}$.

$T_{K}$ is the tensor product of single-link operators

[TABLE]

and is therefore positive-definite as long as $T_{K}^{(1)}$ is. The single-link transfer matrix is defined by

[TABLE]

where $g,\tilde{g}$ are the elements of $G$, $\rho$ is a faithful representation, and $\beta$ is the inverse coupling. Therefore, our proof need only establish that $T^{(1)}_{K}$ is positive definite.

It is sufficient to show that $\langle\Psi|T_{K}^{(1)}|\Psi\rangle>0$ for all $|\Psi\rangle\neq 0$ in $\mathbb{C}G$. We must establish that for any function $f(g)$,

[TABLE]

We work order-by-order in $\beta$ to see that each term is bounded below by [math]:

[TABLE]

where ${}_{n}C_{i}$ are binomial coefficients. Each term indexed by $n,i$ is non-negative:

[TABLE]

where $\{a(b)\}=(a(b)_{1},\cdots,a(b)_{i})$ and $\{c(d)\}=(c(d)_{1},\cdots,c(d)_{n-i})$. It remains only to show that $\langle\Psi|T_{K}^{(1)}|\Psi\rangle=0$ is impossible. For the inequality to be saturated, the following must hold for any $n,i$ and any combination of matrix components in $\{a\},\{b\},\{c\},\{d\}$.

[TABLE]

As the representation $\rho$ is faithful, the space of polynomials on $G$ forms a complete basis, and thus Eq. (53) establishes that the inner product of $f(g)$ with any vector in $\mathbb{C}G$ vanishes, implying that $f(g)=0$.

Therefore, for faithful representations $\rho$ (such as the fundamental representation), $T^{(1)}_{K}$ is positive-definite.