A resource efficient approach for quantum and classical simulations of gauge theories in particle physics

Gauge theories are the basis of high energy physics and the foundation of the standard model (SM). They describe the elementary interactions between particles, which are mediated by the electroweak and strong forces [1, 2, 3], making the SM one of the most successful theories with tremendous predictive power [4]. Still, there are numerous phenomena which cannot be explained by the SM. Examples include the nature of dark matter, the hierarchy of forces and quark masses, the matter antimatter asymmetry and the amount of CP violation [5]. Answering these questions and accessing physics beyond the SM, though, often requires the study of non-perturbative effects.

A very successful approach to address non-pertubative phenomena is lattice gauge theory (LGT) [6, 7, 8], as proposed by Kenneth Wilson in 1974 [9]. In LGT, Feynman’s path integral formulation of quantum field theories (QFTs) is employed on an Euclidean space-time grid. Such a discretised form of the path integral allows for numerical simulations utilizing Markov Chain Monte Carlo (MCMC) methods. The prime target of LGT is quantum chromodynamics (QCD), i.e. the theory of strong interactions between quarks and gluons. In this field, LGT has been extremely successful, allowing for example the computation of the the low-lying baryon spectrum [10], the structure of hadrons, fundamental parameters of the theory and many more [11, 12, 13, 14].

However, many of the aforementioned open questions in modern physics cannot be addressed within the standard approach, due to the sign-problem [15, 16, 17] that renders MCMC methods ineffective. A possible solution is to employ a Hamiltonian formulation of the underlying model. Classical Hamiltonian-based simulations using tensor network states (TNS), including fermionic projected entangled-pair states, have been successful [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], but are so far restricted to mostly one spatial dimension (for link model 2D calculations with DMRG and tree tensor network see e.g [29, 30]). Consequently, there is a necessity for new approaches to both access higher dimensions and address problems where standard MCMC methods fail. It is presently not known whether efficient classical methods can be developed to overcome this problem.

Hamiltonian-based simulations on quantum hardware provide an alternative route, since there is no such fundamental obstacle to simulating QFTs in higher dimensions [32, 33, 34, 35]. Therefore, this approach holds the potential to address questions that cannot be answered with current and even future classical computers. The rapidly evolving experimental capabilities of quantum technologies [36, 37] have led to proof-of-concept demonstrations of simulators tackling one-dimensional theories [38, 39, 40, 41, 42, 43, 44]. Extending these results to higher dimensions is a lively area of research [32, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54], since it represents a crucial step for this field, and realisations on ‘Noisy Intermediate-Scale Quantum’ devices [55, 56], i.e. current quantum hardware, require novel approaches to make this leap.

To meet this challenge, we provide a resource-efficient approach that facilitates the quantum simulation of higher dimensional LGTs that would otherwise be out of reach for current and near-term quantum hardware, which is exemplified Table 1. In addition, purely classical simulations based on the Hamiltonian formalism also benefit from our resource-optimised approach. Hence, we bring both quantum and classical calculations closer to developing computational strategies that do not rely on Monte Carlo methods, and thus circumvent their fundamental limitations.

Our new approach addresses the important problem of reaching the continuum limit (in which the lattice spacing approaches zero) with finite computational resources. Since QFTs are continuous in their time and space variables, the need to take a controlled continuum limit is inherent to any lattice approach and necessary to extract physically relevant results from a lattice simulation.

Taking QCD as a concrete example, we require an accurate description for particles interacting at both short and long distances. Lattice QCD and other LGTs offer the unique tool to investigate both regimes. At long distances, e.g. the bound state spectrum can be computed. At short distances, and after taking the continuum limit, it is possible to connect the perturbative results derived with QFTs with non-perturbative simulations, thus assessing the range in which perturbation theory is valid. However, taking the continuum limit is in general computationally expensive. MCMC methods, for instance, have the intrinsic problem of autocorrelations, that become more and more severe when decreasing the lattice spacing. This drawback in turn leads to a significant increase in the computational cost, and fixes the smallest value of the lattice spacing that can be reached. On one hand, Hamiltonian approaches circumvent this problem. On the other, however, Hamiltonian-based formulations face the challenge that for continuous (Abelian and non-Abelian) gauge groups, local gauge degrees of freedom are defined in an infinite dimensional Hilbert space. As a consequence, any simulation – classical or quantum – requires a truncation of the gauge fields, which is inherently at conflict with the required continuum limit.

In this work, we present a practical solution to overcome this crucial bottleneck and to allow for resource-efficient Hamiltonian simulations of LGTs. Although our approach is general and applicable to LGTs of any dimension, we consider two-dimensional quantum electrodynamics (QED) as a benchmark example.

Truncation of the gauge fields is typically performed in the ‘electric basis’, i.e. the basis in which the electric Hamiltonian and Gauss’ law are diagonal. As such, truncation preserves the gauge symmetry, and the resulting models are known as gauge magnets or link models [57, 58, 59], which are of direct relevance in condensed matter physics [60, 61, 62, 63]. As recently shown in Ref. [54], spin-$1/2$ truncations are within the reach of current quantum simulators. From the perspective of fundamental particle interactions, electric truncations can result in an accurate description of the system in the strong coupling regime. However, by decreasing the value of the coupling or equivalently the lattice spacing, the magnetic contributions to the energy become increasingly important and the number of states that have to be included in the electric basis grows dramatically (a similar increase can be realised by adding an auxiliary spatial dimension to the lattice [64]). An alternative approach to describe the gauge degrees of freedom is to approximate continuous gauge groups with discrete ones, for instance, to approximate $U(1)$ with $\mathbb{Z}_{2L+1}$ ($L\in\mathbb{N}$) [65, 66, 67]. Such approaches also face similar limitations as the ones described above, as $L$ has to be progressively increased with decreasing coupling.

A natural solution to simulate the weak coupling regime consists of exploiting the self-duality [68] of the electric and magnetic terms by Fourier transforming the Hamiltonian and working in the ‘magnetic basis’, i.e. the basis in which magnetic interactions are diagonal, as suggested in [69]. However, the fact that the magnetic degrees of freedom are continuous variables with a gapless spectrum, poses intricate challenges for a resource-efficient truncation scheme, that have yet (to the best of our knowledge) to be addressed. In this work, we provide a practical solution by combining state truncation with a gauge group discretisation that is dynamically adjusted to the value of the coupling. This approach allows for controlled simulations at all values of the bare coupling, smoothly connecting the weak, strong and intermediate coupling regimes. As a proof-of-principle of this new approach and its ability to faithfully simulate non-perturbative phenomena, we target the renormalised coupling in QED in $2+1$ dimensions.

To observe non-perturbative phenomena such as confinement, the simulated physical length scale needs to be larger than the scale at which confinement sets in. As a result, large lattice sizes are required and the number of lattice points grows rapidly when approaching the continuum limit of the theory. This results in computational requirements that cannot be satisfied using current classical and quantum computers. Still, as previously done in the pioneering work by Creutz [70], we can study the bare coupling dependence of the local plaquette operator. This quantity allows us to benchmark our formalism and to show that a smooth connection between the weak and the strong coupling regimes can be established. In addition, our method allows for estimating the precision with which a given truncation approximates the untruncated results.

The paper is organised as follows. In Sec. 2, we review lattice QED in $2+1$ dimensions as an example of Abelian and non-Abelian gauge theories with continuous groups and magnetic interactions. We consider lattices with periodic boundary conditions and reformulate the lattice Hamiltonian in terms of gauge-invariant degrees of freedom. By eliminating redundant variables, we obtain an effective Hamiltonian description that allows for simulations at a low computational cost. In Sec. 3, we introduce a new magnetic representation of lattice QED that is equipped with a regularisation in terms of a $\mathbb{Z}_{2L+1}$ group and an efficient truncation scheme. In Sec. 4, we study the performance of our method and benchmark its precision by calculating the expectation value of the plaquette operator on a periodic plaquette in the static charge limit. We show that both the truncation cut-off parameter, i.e. the maximum number of gauge basis states included in the simulation, and $L$, the dimension of the $\mathbb{Z}_{2L+1}$ group, can be used as adjustable variational parameters. Both are used to optimise the simulation and estimate its accuracy. In the following Sec. 5, we present the generalisation to an arbitrary, two-dimensional periodic lattice with dynamical matter. Finally, we outline the prospects of this method for classical and quantum simulations is in Sec. 6

In this chapter, we provide a Hamiltonian formulation for LGTs with continuous gauge groups that allows for resource-efficient classical and quantum simulations. First, we review the standard Kogut-Susskind Hamiltonian subject to Gauss’ law (the local constraints ensuring gauge invariance) in Sec. 2.1, considering QED on a square lattice as a paradigmatic example. In Sec. 2.2, we proceed to provide a minimal formulation of the QED lattice Hamiltonian, in which redundant degrees of freedom have been removed.

In the following, we describe a scheme that allows switching from the so-called electric representation, where $\hat{H}_{E}$ is diagonal, to the magnetic one, where $\hat{H}_{B}$ is diagonal. Our method is based on the replacement of the $U(1)$ gauge group with the group $\mathbb{Z}_{2L+1}$, and an accompanying transition from the compact formulation to a completely compact formulation, where both field degrees of freedom are treated as compact variables. While this procedure is general, we illustrate it for the minimal periodic system introduced in Sec. 2.2 and consider generalisations in Sec. 5.

Before presenting the scheme, we discuss the following observations about the considered Hamiltonian, that is now reduced to the sum of Eqs. (12) and (2.2), while all charges $\hat{q}_{\boldsymbol{n}}$ in Eq. (2.2) are set to zero. In particular, the lowering (raising) operators $\hat{P}_{x}^{({\dagger})}$ and $\hat{P}_{y}^{({\dagger})}$ acting on the strings are solely contained in the now absent kinetic Hamiltonian in Eq. (22). The total Hamiltonian thus commutes with $\hat{R}_{x}$ and $\hat{R}_{y}$, i.e. $[\hat{H}_{\textrm{gauge}},\hat{R}_{x}]=[\hat{H}_{\textrm{gauge}},\hat{R}_{y}]=0$, and as a consequence the strings become constants of motion. The dynamics induced by the pure-gauge Hamiltonian are thus restricted to different subspaces defined by $\hat{R}_{\mu}|{r_{\mu}}\rangle=r_{\mu}|{r_{\mu}}\rangle$, for $\mu=x,y$. Starting in Sec. 4, we will be interested in a ground state property, therefore we restrict ourselves to the subspace where both strings are confined to the vacuum. The effective Hamiltonian of this subspace can be readily obtained by setting $\hat{R}_{x}=\hat{R}_{y}=0$ in Eqs. (12) and (2.2) which yields

where we introduced the superscript $(\mathrm{e})$ to emphasise is the electric representation.

Since the three rotators possess discrete but infinite spectra, any numerical approach for simulating the Hamiltonian in Eq. (3) requires a truncation of the Hilbert space. In the following, $l$ denotes the cut-off value which is identified by

Thus, the action of the truncated lowering operators is given as

Note that the total dimension of the Hilbert space is reduced to $(2l+1)^{3}$, which is still challenging to simulate even for relatively small values of $l$. In particular, calculations in the weak coupling regime suffer from this severe limitation and until now, no practical methods to solve this issue have been available.

Let us now introduce a formulation that allows for an efficient representation of the Hamiltonian’s eigenstates in the weak coupling regime, where $g\ll 1$. It is based on the exchange of the continuous $U(1)$ group with the discrete group $\mathbb{Z}_{2L+1}$, which provides a discrete basis for the vector potential operators $\hat{A}_{\boldsymbol{n},\boldsymbol{e}_{\mu}}$ and enables a direct transformation into this dual basis via a Fourier transform. The approach is motivated by the key observation that, in the electric representation, the Hamiltonians of the continuous $U(1)$ group and the discrete $\mathbb{Z}_{2L+1}$ group after truncation ($l<L$) are equivalent. The group $\mathbb{Z}_{2L+1}$ consists of $2L+1$ elements, thus the parameter $L$ indicates the size of the Hilbert space. In particular, the rotators $\hat{R}_{j}$ and lowering (raising) operators $\hat{P}_{j}^{({\dagger})}$ ($j=1,2,3$) take the form

The only difference between the truncated $U(1)$ group and untruncated $\mathbb{Z}_{2L+1}$ group is the cyclic property of the lowering (raising) operator, which transforms $|{L}\rangle$ into $|{-L}\rangle$ (and vice versa). However, after a truncation of $\mathbb{Z}_{2L+1}$ with $l<L$, this property is lost, meaning that Eqs. (26) and (27) correspond to each other and the two truncated groups become equivalent.

For now, consider the Hamiltonian which employs the complete $\mathbb{Z}_{2L+1}$ group. Importantly, the relations in Eq. (27) resort to a compact description of the electric field since the spectra of the rotators and strings are constrained to the compact interval $[-L,L]$. We now introduce the following replacement rules for these operators,

which reassemble Fourier series expansions. Crucially, this replacement is exact, i.e. there is no truncation of the Fourier series. Employing the fact that the spectrum of $\hat{R}$ is discrete and takes integer values, the periodicity of the trigonometric functions can be exploited, which allows one to perform a summation over all coefficients where the sine (cosine) is equivalent. Hence, a finite number of $2L$ coefficients remain, which take the form

Here, $\psi_{k}(\bullet)$ is the $k$-th polygamma function. Let us further remark that these rules can be extended to higher powers in the variables $\hat{R}$ than considered in (3).

This replacement turns out to be convenient for the basis transformation explained below. Introducing the convention $|{\boldsymbol{r}}\rangle=|{r_{1}}\rangle|{r_{2}}\rangle|{r_{3}}\rangle$ and recalling Eq. (27), the electric contribution from Eq. (3) reads

Note that we use the notation $\sum_{\boldsymbol{r}=-\boldsymbol{L}}^{\boldsymbol{L}}$ to indicate that the sum collects all combinations of $r_{j}$, where $r_{j}\in[-L,L]$, $j=1,2,3$ and we neglected the constant energy shifts introduced by Eq. (3). The $\mathbb{Z}_{2L+1}$ magnetic Hamiltonian $\hat{H}_{B}^{(\mathrm{e})}$ can be obtained by substituting the cyclic $\hat{P}_{j}$ of Eq. (27) into Eq. (3).

We are now in a position to perform the switch to the dual basis. As shown in App. C, for any $\gamma\in\mathbb{N}$, the discrete Fourier transform $\hat{\mathcal{F}}_{2L+1}$ diagonalises the lowering operators as

Hence, by applying the discrete Fourier transform to the total Hamiltonian we diagonalise the magnetic contributions, while sacrificing the diagonal structure in the electric part, i.e.

and

Note that we introduced the superscript $(\mathrm{b})$, which refers to the magnetic representation of the Hamiltonian. Using this representation, computations in the weak coupling regime $g\ll 1$ can be performed efficiently, as a truncation $l$ now chooses a cut-off for the magnetic field energy. We emphasize that although we employed the rotator formulation of the Hamiltonian, the just presented procedure is likewise valid for the link formulation utilizing the electric field operators. Indeed, the replacement rules in Eq. (3) are then formulated in terms of $\hat{E}$ instead of $\hat{R}$ and inserted into the Hamiltonian in Eq. (2.1). The corresponding magnetic representation is analogously obtained via an application of the Fourier transform.

The parameter $L$ now affects the accuracy of the simulation. In fact, while $L$ is completely irrelevant in the electric representation (truncated $U(1)$ and truncated $\mathbb{Z}_{2L+1}$ are equivalent), it strongly influences the results derived in the magnetic representation. While examining the relationship between $L$ and $l$ in more detail in Sec. 4, we qualitatively discuss our procedure to simulate the $U(1)$ group with the two representations of $\mathbb{Z}_{2L+1}$ in the following. To be more precise, for any $g$ we might always formulate a sequence of approximating representations for any quantum state of the system in the computational basis defined by $|{\boldsymbol{r}}\rangle$, i.e.,

Here, $p^{(\mathrm{e})}$ denotes the expansion coefficients in the electric representation, with the subscript indicating the group to which they are referring to (no subscript stands for the truncated $\mathbb{Z}_{2L+1}$). The first approximation in Eq. (35) is due to the transition from $U(1)$ to the $\mathbb{Z}_{2L+1}$ group, while the second approximation represents the truncation from $(2L+1)^{3}$ down to $(2l+1)^{3}$ states.

The same scheme exists for the magnetic representation, where the weights $p^{(\mathrm{b})}(g,\boldsymbol{r},L)$ are used for the state $|{\psi^{(\mathrm{b})}(g,L)}\rangle$. In this case, however, the choice of $L$ is important. While the truncated electric representation directly corresponds to the truncated and compact $U(1)$ formulation, the completely compact formulation employed in the magnetic representation is crucially affected by the level of discretisation $L$. This relation is examined further in Sec. 4.2, where we study the convergence of the two representations to $U(1)$ for intermediate values of the coupling $g$. Hence, in the remainder of this manuscript we consider the completely compact formulation for the magnetic representation only and resort to the compact formulation for the electric representation, i.e. to Eq. (3) for the case of pure gauge.

The interplay of the parameters $L$ and $l$ can be intuitively understood by employing a geometrical illustration. In Fig. 2, the black circles represent the continuous $U(1)$ group, which is approximated by $2L+1$ possible states (blue lines) of the $\mathbb{Z}_{2L+1}$ group. For $l=L$, we faithfully describe the untruncated $\mathbb{Z}_{2L+1}$ group, and use both the solid and the dashed blue lines in the figure. By choosing $l<L$, we select the states marked with solid blue lines that lie symmetrically around $|{\boldsymbol{r}=\boldsymbol{0}}\rangle$ and achieve a binned approximation of any continuous $p_{U(1)}$ lying in the grey area. Furthermore, for any fixed $l$, the parameter $L$ controls the spread of the available basis states (or bins) around $|{\boldsymbol{0}}\rangle$. Since we are interested in the convergence of the truncated $\mathbb{Z}_{2L+1}$ to $U(1)$ which occurs for $L\to\infty$, we only consider the $2l+1$ states that are important for the dynamics. In particular, we disregard cyclic effects from the lowering operator $\hat{P}$ that are a distinctive feature of $\mathbb{Z}_{2L+1}$ with respect to $U(1)$ [see Eq. (27)].

As an example of this relationship between $l$ and $L$, consider the two distributions drawn from $U(1)$, represented by the red and green dashed lines in Fig. 2. Clearly, the combination $L=3$, $l=1$ is insufficient to approximate the broad red distribution. Hence, we increase $l$ to completely cover the target distribution within the grey shaded area. A reduction of $L$ is also a practicable option, especially given the situation where $l$ could not be increased further because of a lack of computational resources. By raising $L$ instead, our binned approximation has a higher resolution around the zero state. For the more localised green state, therefore, it is advantageous to choose the higher value $L=6$ when $l=2$ is an available option. In fact, the combination $L=6$, $l=1$ leads to a worse approximation of the green state than the choice $L=3$, $l=1$, which is therefore preferable when $l$ is limited to unity.

To that end, it is clear that the interplay of $L$ and $l$ represents a crucial point when estimating results and their error due to the performed discretisation. Note that the spread of the distribution determines the value of the free parameter $L$, while $l$ will, for all practical purposes, be limited by the amount of physical resources, e.g. by the memory of a classical computer or the number of available qubits in a quantum simulator.

In the previous section we outlined the transformation from the electric to the magnetic representation, suited to describe the strong and the weak coupling limits, respectively. Here, we develop a protocol which allows for assessing convergence of the truncated representations.

First, we qualitatively describe the system’s behaviour at different values of the bare coupling, which will be useful for motivating the protocol. Later, we consider the non-asymptotic cases, where it is not known whether a given truncation is sufficient to describe the considered system with the desired precision. In our example dealing with a pure gauge $U(1)$ theory, this happens for $g\approx 1$, where both representations might be inaccurate. Finally, we discuss convergence in the weak and strong coupling regimes, respectively and apply our protocol to estimate the plaquette expectation value. In particular, we consider the interplay between the parameters $l$ and $L$, which plays an important role when $g\gg 1$.

From now on, we work within a unit lattice spacing, i.e. $a=1$, but emphasise that this represents no restriction for the following results.