Gluon Field Digitization for Quantum Computers

TL;DR

This paper introduces a digitization scheme for $SU(3)$ gauge theories suitable for quantum computers, enabling accurate classical and potential quantum simulations of QCD with minimal qubits.

Contribution

A novel digitization method using the subgroup $S(1080)$ for $SU(3)$ gauge theories, with a modified action allowing simulations in the scaling window.

Findings

Classical Monte Carlo results match full $SU(3)$ continuum extrapolations.

The digitization scheme achieves sub-percent precision at multiple lattice spacings.

The method is suitable for NISQ-era quantum chromodynamics simulations.

Abstract

Simulations of gauge theories on quantum computers require the digitization of continuous field variables. Digitization schemes that uses the minimum amount of qubits are desirable. We present a practical scheme for digitizing $SU(3)$ gauge theories via its discrete subgroup $S(1080)$. The $S(1080)$ standard Wilson action cannot be used since a phase transition occurs as the coupling is decreased, well before the scaling regime. We proposed a modified action that allows simulations in the scaling window and carry out classical Monte Carlo calculations down to lattice spacings of order $a\approx 0.08$ fm. We compute a set of observables with sub-percent precision at multiple lattice spacings and show that the continuum extrapolated value agrees with the full $SU(3)$ results. This suggests that this digitization scheme provides sufficient precision for NISQ-era QCD simulations.

Tables

Table 1. Table 1: Wilson flow parameters t 0 / a subscript 𝑡 0 𝑎 \sqrt{t_{0}}/a and t 0.2 / a subscript 𝑡 0.2 𝑎 \sqrt{t_{0.2}}/a found on lattices of size ( 3 ​ N t ) 4 superscript 3 subscript 𝑁 𝑡 4 (3N_{t})^{4} along our trajectory where N t subscript 𝑁 𝑡 N_{t} is temporal size used to determine T c subscript 𝑇 𝑐 T_{c} . In the last two columns, the first error is statistical, and the second is from the separatrix.

$N_t$	$N_s$	${\beta }_0^c$	${\beta }_1^c$	$\sqrt{t_0}/a$	$\sqrt{t_{0.2}}/a$
4	12	9.154(2)	-0.9061(3)	1.016(3)(3)	0.8316(12)(20)
6	18	12.795(9)	-1.3673(11)	1.508(3)(5)	1.2068(18)(42)
8	24	19.61(4)	-2.231(5)	2.000(4)(8)	1.595(3)(6)

Equations

$$S=-\sum _p\left(\frac{{\beta }_0}{3}\mathrm{Re}\mathrm{Tr}{U}_p+{\beta }_1\mathrm{Re}\mathrm{Tr}{U}_p^2\right),$$

$$P_{ij}({\varphi }_i,{\varphi }_j)=\text{min}\left(1,e^{({\beta }_{0,i}-{\beta }_{0,j})(\tilde{S}[{\varphi }_i]-\tilde{S}[{\varphi }_j])}\right).$$

$${\beta }_1=-0.1267{\beta }_0+0.253,$$

$${\dot{V}}_t(x,\mu )=-\frac{1}{{\beta }_0}$$

$$V_t(x,\mu ){|}_{t=0}$$

$${\left(t^2⟨E⟩\right)}_{t=t_X}=X$$

Full text

Quantum computers can attack problems in physics which appear intractable on classical computers Feynman (1982). Large-scale quantum computers would allow simulations of non-perturbative calculations of real time evolution and finite-density equations of state. For the foreseeable future, though, quantum computers will be limited to tens or hundreds of non-error-corrected qubits with circuit depths less than a thousand gates — the so-called Noisy Intermediate-Scale Quantum (NISQ) era. QCD simulations on quantum computers – especially in the NISQ era – depend upon formulating QCD in an efficient way.

Fermionic fields like quarks can be easily digitized as qubit registers by encoding their presence or absence in a given state Jordan and Wigner (1928); Verstraete and Cirac (2005); Zohar and Cirac (2018); Whitfield et al. (2016). This is evident from the few existing calculations performed on quantum computers Martinez et al. (2016); Klco et al. (2018); Lamm and Lawrence (2018); Shehab et al. (2019). The continuous nature of gauge fields preclude such exact digitization. Proposed solutions involve either eliminating the bosonic fields using some model-dependent properties or truncating in occupation number Hackett et al. (2018); Macridin et al. (2018); Yeter-Aydeniz et al. (2018); Klco and Savage (2018); Bazavov et al. (2015); Zhang et al. (2018); Unmuth-Yockey (2018); Unmuth-Yockey et al. (2018); Zache et al. (2018); Raychowdhury and Stryker (2018); Kaplan and Stryker (2018); Stryker (2019); Alexandru et al. (2019).

The situation is reminiscent of the pioneering days of lattice field theory when computer memory was limited and the cost of storing $SU(3)$ elements was prohibitive. Several attempts were made to replace the continuous gauge fields by a discrete set of values Petcher and Weingarten (1980); Jacobs and Rebbi (1981); Bhanot and Rebbi (1981); Grosse and Kuhnelt (1981); Bhanot (1982); Lisboa and Michael (1982); Flyvbjerg (1984a, b). Quantum computation is presently in a similar situation where every qubit comes at a high cost. Representing each gauge link by 9 complex numbers represented using double-precision format requires 1152 qubits. In contrast, the largest “crystal-like” discrete subgroup of $SU(3)$, $S({1080})$, contains 1080 elements and would require only 11 qubits to store each link value.

Digitization typically reduces the symmetry of the model. With this reduction, it is not a given that the original model is recovered in the continuum limit as the universality class of the lattice model may change Hasenfratz and Niedermayer (2001a); Caracciolo et al. (2001a); Hasenfratz and Niedermayer (2001b); Patrascioiu and Seiler (1998); Krcmar et al. (2016); Caracciolo et al. (2001b). For any discrete group, there is always a finite difference in the action, $\Delta S$, between the field configurations with the two smallest actions, as opposed to continuous groups where no such gap exists. This may lead to “freezing” at some critical $\beta_{f}$; that is, all field values except the identity (and gauge-equivalents) are exponentially suppressed. For values of $\beta$ beyond $\beta_{f}$ the theory with the discrete group differs drastically from the continuous group and is no longer a reasonable approximation. This is a particular problem for asymptotically-free theories like QCD. The spacetime continuum limit where the lattice spacing $a$ approaches zero is obtained by making $\beta$ large but that is where the continuous and discrete group theories differ. This is not necessarily fatal: realistic lattice calculations are performed on classical computers with a finite $a$ and extrapolated in a controlled manner to $a\rightarrow 0$. In these calculations, what is required is that $a$ is smaller than typical hadronic scales (e.g. the size of hadrons). Typically values used in state-of-art calculations are $\mathcal{O}(0.1$ fm). This corresponds $\beta\gtrsim 6$ when using the Wilson action, in the so-called scaling region. Our goal will be to set up a framework where discrete groups can be used to reproduce $SU(3)$ results in the scaling region, such that continuum extrapolations can be performed.

There were a number of early studies of the viability of crystal-like discrete subgroups of $U(1)$ Creutz et al. (1979); Creutz and Okawa (1983) and $SU(N)$ Bhanot and Rebbi (1981); Petcher and Weingarten (1980); Bhanot (1982) gauge theories. While the discrete subgroups all have freezing transitions for the Wilson action, $\beta_{f}$ increases with the size of the subgroup. For $U(1)$, the $\mathbb{Z}_{N}$ theories have a $\beta_{f}$ in the scaling regime (in this case $\beta\gtrsim 1$) for $N>4$. $SU(2)$ has only three crystal-like subgroups: the binary tetrahedral, $\mathbb{BT}$, the binary octahedral, $\mathbb{BO}$, and the binary icosohedral, $\mathbb{BI}$. Using the Metropolis algorithm with 100 measurements separated by 1000 updates, we refined the results of Petcher and Weingarten (1980), finding that while $\mathbb{BT}$ has $\beta_{f}=2.24(8)$, $\mathbb{BO}$ and $\mathbb{BI}$ have $\beta_{f}=3.26(8)$ and $\beta_{f}=5.82(8)$ respectively, both deep in the scaling regime $\beta\gtrsim 2.2$. Hence, these two last groups can be used in lieu of $SU(2)$ for practical calculations.

The story changes for $SU(3)$. There are five crystal-like subgroups: $S(60)$, $S(108)$, $S(216)$, $S(648)$, and $S({1080})$, designated by their number of elements. For all these, $\beta_{f}<6$ , with the largest, $S({1080})$, being reported to have $\beta_{f}=3.58(2)$ obtained on a $2^{4}$ lattice Bhanot and Rebbi (1981). Our own calculations on a $2^{4}$ volume with larger statistics for $S({1080})$ show a slightly larger value $\beta_{f}=3.935(5)$. In any case, it is evident that the $S({1080})$ theory with the Wilson action is inadequate to reach the scaling regime. Subsequent work Lisboa and Michael (1982) showed that extending the elements to include the midpoints between elements of $S({1080})$ was sufficient to push $\beta_{f}\approx 7$. However this requires more bits and sacrifices the group structure.

To overcome these limitations, attempts were made to approximate the $SU(3)$ Wilson action by a modified action based on a subgroup Edgar (1982); Bhanot (1982); Creutz and Okawa (1983); Fukugita et al. (1983); Horn et al. (1982); Flyvbjerg (1984a, b); Ayala and Baig (1990), although only in Bhanot (1982) were Monte Carlo calculations undertaken. There, simulations using $S(648)$ with a Wilson action modified by a $|\operatorname{{Tr}}U_{p}|^{2}$ term (equivalent to the trace in the adjoint representation) were performed. Even with this modified action $S(648)$ was inadequate to reach the scaling regime. Further, it was conjectured, based on small scale simulations and mean-field estimates, that $S({1080})$ would also fail to reach the scaling regime with that modified action. Calculations with modified actions of $SU(3)$ were also performed to study thermodynamics and reduce lattice spacing errors Blum et al. (1995); Heller (1995, 1996); Hasenbusch and Necco (2005, 2004).

With this history in mind, we study the viability of $S({1080})$ with a different action

[TABLE]

where $U_{p}\in S(1080)$ indicates a plaquette, and the first term has been normalized such that, for $\beta_{1}=0$, the action matches the $SU(3)$ Wilson action (with $\beta=\beta_{0}$). Simulations can be performed efficiently by employing precomputed multiplication and trace tables.

One could argue heuristically that the action in Eq. (1) will lead to the same continuum limit as, for instance, the Wilson action of $SU(3)$ ($S_{SU(3)}=-\tfrac{\beta_{SU(3)}}{3}\sum_{p}\operatorname{{Re}}\operatorname{{Tr}}U_{p}$) by noting that the continuum limit in asymptotically-free theories is obtained by setting the coupling to be small ($\beta$ large) and that in this limit only small field fluctuations are important. For small values of the field $S$ and $S_{SU(3)}$ agree as long as $\beta_{SU(3)}$ is a certain linear combination of $\beta_{0},\beta_{1}$. The flaw with this argument is that as $\beta_{SU(3)}$ is increased towards the continuum limit, so do $\beta_{0},\beta_{1}$ and, at some point, the fluctuations become smaller than the separation between the identity and the nearest element to it in $S({1080})$ and fields freeze. Beyond this point it makes no sense to consider small fluctuations in the $S({1080})$ theory. This is seen dramatically in Fig. 1 where the average plaquette of the Wilson action $SU(3)$ theory (black) and the $S({1080})$ action $S$ at $\beta_{1}=0$ (red) are shown to agree very well until they abruptly diverge when freezing occurs. This agreement is expected since the strong coupling expansion predicts the difference in the average plaquette for $SU(3)$ and $S(1080)$ with $\beta_{1}=0$ to be of order ${\cal O}(10^{-6}\beta^{5})$ Bhanot and Rebbi (1981). By setting $\beta_{1}<0$ the gap between the two lowest action values is reduced and freezing occurs only at large $\beta_{0}$ (see blue curve in Fig. 1.) Our proposal is to find a trajectory in the ($\beta_{0},\beta_{1}$) plane that avoids freezing and allows us to get closer to the continuum limit. The idea is that as we move on this trajectory towards larger values of $\beta_{0}$ we produce configurations with larger correlation lengths, ideally increasing all the way to infinity. We must then check that this action generates the same physics as $SU(3)$.

In order to guide our choice of trajectory we roughly mapping out the phase diagram of the $S({1080})$ theory in the ($\beta_{0},\beta_{1}$) plane to determine where the theory is frozen and not useful to approximate $SU(3)$. The identification of this boundary is complicated because each phase exists as a metastable state throughout the phase diagram. To deal with this metastability and the associated long autocorrelation times, we used a parallel tempering algorithm similar to Swendsen and Wang (1986); Earl and Deem (2005). We perform simulations with a set of $\{\beta_{0,i}\}$ for a fixed $\beta_{1}$. For every fifth local update, the configuration with $\beta_{0,i}$ is randomly selected and swapped with the ensemble with $\beta_{0,j}$ where $j=i\pm 1,i\pm 2$ with probability:

[TABLE]

With this algorithm, we were able to map out the full $(\beta_{0},\beta_{1})$ space by searching for peaks in the susceptibilities $\chi_{0}=\frac{\partial^{2}\langle S\rangle}{\partial\beta_{0}^{2}}$ and/or $\chi_{1}=\frac{\partial^{2}\langle S\rangle}{\partial\beta_{1}^{2}}$ (see Fig. (2)). Besides the freezing transition (shown in red), there are additional transitions shown in the upper left corner of the phase diagram delineating regions where the dynamics is partially frozen down to a subset of the group elements. This rich structure is qualitatively similar to the one found in the $SU(3)$ theory but it is of little concern to us. We focus instead on the lower ($\beta_{1}<0$) region of the phase diagram. At $\beta_{1}=0$ we measure $\beta_{c}=3.935(5)$, a larger value than the $\beta_{c}=3.58(2)$ found with smaller statistics in Bhanot and Rebbi (1981).

With the freezing transition mapped out, we can propose a way to approach the continuum limit by performing simulations along a trajectory that avoids the frozen phase. We choose the trajectory

[TABLE]

shown as a dotted line in Fig. 2. A pure gauge theory (without matter) can predict only ratios of observables. So, one observable should be used to set the scale, that is, to determine the dimensionful quantity $a$ at a fixed $\beta_{0},\beta_{1}$. After this, the theory can make predictions for all other observables. In order to quantify the approach to the continuum and compare it with the expected continuum result we then need a minimum of two observables: one to set the scale and another to be predicted. We choose the pseudocritical temperature $T_{c}$ for deconfinement and the scale $t_{0}$ defined by the Wilson flow Luscher (2010).

We determine $T_{c}$ by looking at the distribution of values of the Polyakov loop $P$. In the confined phase they concentrate around zero, while at higher temperatures cluster around the vertices of a triangle. We denote by $w_{c},w_{3}$ the number of configurations in which $P$ lies near zero or near one of the vertices of the triangle. Following the procedure outlined in Francis et al. (2015), we label as “center” configurations the ones that lie inside an equilateral triangle centered at the origin and rotated such that one side is perpendicular to the positive real axis. The separatrix is adjusted to the minimum of the histogram of the $\operatorname{{Re}}\operatorname{{Tr}}P$. This position is determined from simulations close to the transition point. There is some ambiguity in this definition but we include this variability into our error budget. We then look at the quantity $s(\beta_{0},\beta_{1})=\frac{3w_{c}-w_{3}}{3w_{c}+w_{3}}$: it is $+1$ deep in the confined phase and reaches -1 deep in the deconfined phase. When $s=0$, the theory is tuned to $T_{c}$. This definition of $T_{c}$ has been shown to have finite-volume effects that scale exponentially in the spatial volume for $SU(3)$, unlike the peak of $\chi_{0}$ which exhibits power law volume corrections Francis et al. (2015). We found that for $N_{t}=4$, the difference in $(\beta_{0}^{c},\beta_{1}^{c})$ for $N_{s}=12,16$ was negligible, confirming the same behavior in $S({1080})$. Therefore, we assumed $N_{s}=3N_{t}$ has negligible finite-volume effects at larger $N_{t}$ and use this volume to compute $T_{c}$. For each set of parameters we collected ${\cal O}(10^{6})$ measurements separated by 10 sweeps. To perform a sweep, we visit each link and update it using a multi-hit Metropolis step.

For $N_{t}=4,6,8$ we scan ($\beta_{0},\beta_{1}$) along the trajectory Eq. (3) to find ($\beta^{c}_{0},\beta^{c}_{1}$) for which $s=0$. The values obtained are listed on Table 1. For each of these the inverse physical temperature $1/T=N_{t}a(\beta_{0}^{c},\beta_{1}^{c})=1/T_{c}$ is the same.

We compute the scale $t_{0}$ defined by the Wilson flow for these sets of $(\beta^{c}_{0},\beta^{c}_{1})$ parameters. We first generate configurations on lattices of size $(3N_{t})^{4}$ where $N_{t}$ is the temporal size used to determine $T_{c}$. For each ensemble we generate 200 configurations, separated by 1000 sweeps. The configurations generated are represented as $SU(3)$ matrices and used as initial conditions for performing the Wilson flow Luscher (2010):

[TABLE]

where $S_{W}[V_{t}]$ is the Wilson action of $SU(3)$ fields $V_{t}$ at some Wilson-flow time $t$. Using the flow, we define two observables $t_{X=0.2}$ and $t_{X=0.3}$ implicitly by the expression

[TABLE]

where $X=0.2,0.3$, and $E$ is the lattice clover definition of the energy density. Following the convention, $t_{X=0.3}$ is called $t_{0}$. Both of these observables have been measured precisely for $SU(3)$ pure gauge theory, allowing for comparison. By also measuring $t_{0.2}$, we probe higher energy scales where larger discrepancies between $S({1080})$ and $SU(3)$ should be found. Our results are found in Table 1.

In the absence of discretization effects, the value of $t_{0}$ in physical units should be the same on all our lattices. We demonstrate that the variation of $t_{0}$ as we approach the continuum is mild and the extrapolated value agrees with the full $SU(3)$ result. With our data, it is possible to construct a dimensionless quantity, $T_{c}\sqrt{t_{X}}$ which can be compared to those of $SU(3)$ at both finite lattice spacing $a$ and by extrapolating to the continuum. Using a linear extrapolation, we compute a continuum value of $T_{c}\sqrt{t_{0}}=0.2489(11)$ which is in agreement with $T_{c}\sqrt{t_{0}}=0.2489(14)$ Francis et al. (2015) and $T_{c}\sqrt{t_{0}}=0.2473(7)$ Kitazawa et al. (2016) computed in full $SU(3)$. Similarly, our extrapolated value of $\frac{\sqrt{t_{0}}}{\sqrt{t_{0.2}}}=0.1269(6)$ is in good agreement with the value of $0.1264(4)$ computed for $SU(3)$ Asakawa et al. (2015). Our results for $T_{c}\sqrt{t_{0}}$ are compared to Francis et al. (2015) in Fig. 3. It is interesting to note that the $\mathcal{O}(a^{2})$ corrections appear milder for the modified action used here compared to the Wilson action of the $SU(3)$. This feature of modified actions has been observed previously in $SU(3)$ Blum et al. (1995); Heller (1995, 1996); Hasenbusch and Necco (2005, 2004), which suggest further benefits of using this action for quantum simulations.

Assuming our $S(1080)$ action has the same continuum limit as $SU(3)$, it is possible derive physical values for $t_{0}$ and a lattice spacing $a$ for each of our ensembles. Using the $SU(3)$ relations of Asakawa et al. (2015); Kitazawa et al. (2016) between $w_{0.4}$, $T_{c}$, $\Lambda_{SU(3)}$ and $r_{0}=0.49(4)$ fm, we obtain $\sqrt{t_{0}}\approx 0.16(2)$ fm where our error is dominated by $r_{0}$, agreeing with $\sqrt{t_{0}}=0.1638(10)$ fm determined in $SU(3)$ Sommer (2014). With this, our $24^{4}$ lattice has $a\approx 0.08$ fm or 2.5 GeV*-1*. This suggests it would be possible to extract and compare glueball states Morningstar and Peardon (1999); Chen et al. (2006) or quenched calculations of hadron masses Aoki et al. (2003) to $SU(3)$ values with sub-percent precision.

In this work, we have presented results for the discrete gauge group $S({1080})$ using a modified Wilson action. After mapping the entire phase diagram of the new action, we found a parameter trajectory that avoids the freezing transition, allowing for calculations with lattice spacing $a\approx 0.08$ fm. We have shown that this action is capable of reproducing the physics of $SU(3)$ below 2.5 GeV*-1* by measuring $T_{c}\sqrt{t_{0}}=0.2489(11)$ and $\frac{\sqrt{t_{0}}}{\sqrt{t_{0.2}}}=0.1269(6)$, which agree to remarkable precision with the full group.

The qubit savings from using $S({1080})$ instead of $SU(3)$ are dramatic. In the NISQ era, where small lattice sizes and noisy gates will likely dominate the error, the parameters used in this paper should be a sufficient approximation of $SU(3)$. However, if gluon actions are required at $a<0.08$ fm, the action in Eq. (1) on the trajectory specified by Eq. (3) may be insufficient since $a$ on this trajectory seems to have a minimum. Another trajectory nearer to the freezing transition may provide a smaller lattice spacing. Future work is required to determine whether we can generate arbitrarily small $a$ on a different trajectory with our action or if additional terms are required. Furthermore, questions of how well other observables like hadronic spectra are reproduced and the effect of including of fermions are left for future studies. Given that circuit depth is of concern in the NISQ era, a broad search in modified action space should be undertaken with an eye toward terms that require few quantum gates while still efficiently reaching the continuum limit of $SU(3)$. Quantitative comparisons of these actions will require constructing the quantum gates for $S({1080})$ Lamm et al. (2019).