Gluon Field Digitization for Quantum Computers
Andrei Alexandru Paulo F. Bedaque Siddhartha Harmalkar Henry Lamm Scott Lawrence Neill C. Warrington






aalexan@gwu.edu
Department of Physics, The George Washington University, Washington, D.C. 20052, USA
Department of Physics, University of Maryland, College Park, MD 20742, USA

bedaque@umd.edu
Department of Physics, University of Maryland, College Park, MD 20742, USA

sharmalk@umd.edu
Department of Physics, University of Maryland, College Park, MD 20742, USA

hlamm@umd.edu
Department of Physics, University of Maryland, College Park, MD 20742, USA

srl@umd.edu
Department of Physics, University of Maryland, College Park, MD 20742, USA

ncwarrin@umd.edu
Department of Physics, University of Maryland, College Park, MD 20742, USA


NuQS Collaboration


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 .  
The  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≈ 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.



Gluon Field Digitization for Quantum Computers
    Neill C. Warrington
    December 30, 2023
==============================================


Quantum computers can attack problems in physics which appear intractable on classical computers <cit.>.  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 <cit.>. This is evident from the few existing calculations performed on quantum computers <cit.>. 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 <cit.>. 

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 <cit.>.  
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), , 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 <cit.>. For any discrete group, there is always a finite difference in the action, Δ 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 β_f; that is, all field values except the identity (and gauge-equivalents) are exponentially suppressed.  For values of β beyond β_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 β 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→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 𝒪(0.1 fm). This corresponds β 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) <cit.> and SU(N) <cit.> gauge theories.  While the discrete subgroups all have freezing transitions for the Wilson action, β_f increases with the size of the subgroup.  For U(1), the  theories have a β_f  in the scaling regime (in this case β 1) for  N>4.  SU(2) has only three crystal-like subgroups: the binary tetrahedral, 𝔹𝕋, the binary octahedral, 𝔹𝕆, and the binary icosohedral, 𝔹𝕀.  Using the Metropolis algorithm with 100 measurements separated by 1000 updates, we refined the results of <cit.>, finding that while 𝔹𝕋 has β_f=2.24(8), 𝔹𝕆 and 𝔹𝕀 have β_f=3.26(8) and β_f=5.82(8) respectively, both deep in the scaling regime β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 , designated by their number of elements.  For all these, β_f < 6 , with the largest, , being reported to have β_f=3.58(2) obtained on a 2^4 lattice <cit.>.  Our own calculations on a 2^4 volume with larger statistics for  show a slightly larger value β_f=3.935(5). In any case, it is evident that the  theory with the Wilson action is inadequate to reach the scaling regime. Subsequent work <cit.> showed that extending the elements to include the midpoints between elements of  was sufficient to push β_f≈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 <cit.>, although only in <cit.> were Monte Carlo calculations undertaken. 
There, simulations using S(648) with a Wilson action modified by a | 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  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 <cit.>.  



With this history in mind, we study the viability of  with a different action


    S=-∑_p (β_0/3 U_p +β_1  U_p^2) ,

where U_p ∈ S(1080) indicates a plaquette, and the first term has been normalized such that, for β_1=0, the action matches the SU(3) Wilson action (with β=β_0).  




Simulations can be performed efficiently by employing precomputed multiplication and trace tables.

One could argue heuristically that the action in Eq. (<ref>) will lead to the same continuum limit as, for instance, the Wilson action of SU(3) (S_SU(3) = -β_SU(3)3∑_p  U_p) by noting that the continuum limit in asymptotically-free theories is obtained by setting the coupling to be small (β 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 β_SU(3) is a certain linear combination of β_0, β_1. The flaw with this argument is that as β_SU(3) is increased towards the continuum limit, so do β_0, β_1 and, at some point, the fluctuations become smaller than the separation between the identity and the nearest element to it in  and fields freeze. Beyond this point it makes no sense to consider small fluctuations in the  theory. This is seen dramatically in Fig. <ref> where the average plaquette of the Wilson action SU(3) theory (black) and the  action S at β_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 β_1=0 to be of order O(10^-6β^5) <cit.>. By setting β_1<0 the gap between the two lowest action values is reduced and freezing occurs only at large β_0 (see blue curve in Fig. <ref>.) 

Our proposal is to find a trajectory in the (β_0, β_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 β_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  theory in the (β_0, β_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 <cit.>.  We perform simulations with a set of {β_0,i} for a fixed β_1.  For every fifth local update, the configuration with β_0,i is randomly selected and swapped with the ensemble with β_0,j where j=i±1,i±2 with probability:

    P_ij(ϕ_i,ϕ_j)=min(1,e^(β_0,i-β_0,j)(S̃[ϕ_i]-S̃[ϕ_j])).


With this algorithm, we were able to map out the full (β_0,β_1) space by searching for peaks in the susceptibilities  χ_0=∂^2⟨ S⟩/∂β_0^2 and/or χ_1=∂^2⟨ S⟩/∂β_1^2 (see Fig. (<ref>)).  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 (β_1 <0) region of the phase diagram.
At β_1=0 we measure β_c=3.935(5), a larger value
than the β_c=3.58(2) found with smaller statistics in <cit.>. 



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

    β_1=-0.1267β_0+0.253,
 shown as a dotted line in Fig. <ref>.
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 β_0, β_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 <cit.>. 

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 <cit.>, 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 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(β_0, β_1)=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 χ_0 which exhibits power law volume corrections <cit.>.  We found that for N_t=4, the difference in (β_0^c,β_1^c) for N_s=12,16 was negligible, confirming the same behavior in . 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 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 (β_0, β_1) along the trajectory Eq. (<ref>) to find (β^c_0, β^c_1) for which s=0. The values obtained are listed on Table <ref>. For each of these the inverse physical temperature 1/T = N_t a(β_0^c, β_1^c) = 1/T_c is the same.




We compute the scale t_0 defined by the Wilson flow for these sets of (β^c_0, β^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 <cit.>: 

    V̇_t(x,μ)=-1/β_0   (∂_x,μS_W[V_t])V_t(x,μ),
    V_t(x,μ)|_t=0   =U(x,μ)

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

    (t^2⟨ E⟩)_t=t_X=X

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  and SU(3) should be found.  Our results are found in Table <ref>. 

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√(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√(t_0)=0.2489(11) which is in agreement with T_c√(t_0)= 0.2489(14) <cit.> and T_c√(t_0)= 0.2473(7) <cit.> computed in full SU(3). Similarly, our extrapolated value of √(t_0)/√(t_0.2)=0.1269(6) is in good agreement with the value of 0.1264(4) computed for SU(3) <cit.>.
 Our results for T_c√(t_0) are compared to <cit.> in Fig. <ref>.  It is interesting to note that the 𝒪(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) <cit.>, 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 <cit.> between w_0.4, T_c, Λ_SU(3) and r_0=0.49(4) fm, we obtain √(t_0)≈0.16(2) fm where our error is dominated by r_0, agreeing with √(t_0)=0.1638(10) fm determined in SU(3) <cit.>. With this, our 24^4 lattice has a≈0.08 fm or 2.5 GeV^-1.  This suggests it would be possible to extract and compare glueball states <cit.> or quenched calculations of hadron masses <cit.> to SU(3) values with sub-percent precision.


In this work, we have presented results for the discrete gauge group  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≈ 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√(t_0)=0.2489(11) and √(t_0)/√(t_0.2)=0.1269(6), which agree to remarkable precision with the full group.

The qubit savings from using  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. (<ref>) on the trajectory specified by Eq. (<ref>) 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  <cit.>.



 The authors thank Y. Yamauchi and Y. Cai for insight into quantum gates.  A. A. is supported in part by the National Science Foundation CAREER grant PHY-1151648 and by U.S. Department of Energy under Contract No. DE-FG02-95ER-40907. A.A. acknowledges the hospitality of the University of Maryland where part of this work was performed.  P. B., S.H., H. L., and S. L. are supported by the U.S. Department of Energy under Contract No. DE-FG02-93ER-40762. H. L. acknowledges the hospitality of U.W. where parts of this work were produced.

apsrev4-1

