Gluon Field Digitization via Group Space Decimation for Quantum Computers

Efficient digitization is required for quantum simulations of gauge theories. Schemes based on discrete subgroups use fewer qubits at the cost of systematic errors. We systematize this approach by deriving a single plaquette action for approximating general continuous gauge groups through integrating out field fluctuations. This provides insight into the effectiveness of these approximations, and how they could be improved.


I. INTRODUCTION
Large-scale quantum computers can simulate nonperturbative quantum field theories which are intractable classically [1]. Alas, Noisy Intermediate-Scale Quantum (NISQ) era systems will be limited both in qubits and circuit depths. Whether any gauge theory simulations in this period are possible depends upon efficient formulations. The situation is similar to the early days of lattice field theory when computer memory was limited and the cost of storing SU (3) elements was prohibitive.
In this work, we systematize the proposal of replacing continuous gauge groups G by their discrete subgroups H [11,28] by deriving the lattice actions using the group space decimation procedure of [36,37]. After deriving the general third order action, we will investigate the behaviour of discretizing three distinct gauge groups U (1), SU (2), and SU (3). We begin with a review of the discrete subgroup approximation in Sec. II. In Sec. III we discuss the general aspects of the group space decimation procedure. Following that, in Sec. IV we derive the decimated action up to 3rd order. Sec. V the continuous group limit of the procedure is studied, and we conclude in Sec. VI.
Here, we move toward improving the systematic errors of one digitization scheme-approximating gauge theories by replacing G → H [11,28]. The best candidate to replace SU (3) is its largest crystal-like discrete subgroup, the Valentiner group V with 1080 elements 1 . As a rough comparison, a SU (3) gauge link represented by 9 complex numbers using double-precision format requires 9 × 2 × 64 = 1152 qubits. In contrast, V might need as few as 11 qubits per link. This digitization scheme was explored in the early days of Euclidean lattice field theory. The viability of the Z n subgroups replacing U (1) were studied in [42,43]. Further studies of the crystal-like discrete subgroups of SU (N ) were performed [40,44,45], including with fermions [46,47]. These studies met with mixed success depending on the group and action tested.
The fundamental issue of group discretization can be understood by considering the Wilson gauge action where U p indicates a plaquette of continuous group gauge links U (for discrete groups, we will denote plaquettes by u p and links by u). As β → ∞, gauge links near the group identity 1 dominate, i.e. U ≈ 1 + ε, where ε can be arbitrarily small. Therefore the gap ∆S = S[1 + ε] − S[1] goes to zero smoothly. For discrete groups, ε has a minimum given by the nearest elements N to 1, and thus ∆S = S[N ] − S[1] > 0. This strongly suggests a phase transition at some critical β f = c/∆S, where c ≈ O(1) depends on spacetime dimensionality, gauge group, and entropy. For U (1) → Z n in 4d, β f = 0.78 1−cos(2π/n) [45]. Above β f , all field configurations but u = 1 are exponentially suppressed. Thus, H fails as an approximation for G for β > β f . Another way to understand this behavior follows [48], where it was shown that the discrete theories are equivalent to continuous groups coupled to a Higgs field. The Higgs mechanism introduces a new phase missing from the continuous gauge theory when β → ∞.
Both arguments suggest the discrete group can be thought of as an effective field theory for the continuous group with a UV-cutoff at Λ f . Provided the typical separation of scales of physics m IR Λ f , the approximation could be reliable up O(m IR /Λ f ) effects.
In lattice calculations, one replaces the UV cutoff by a fixed lattice spacing a = a(β) which shrinks with increasing β for theories with asymptotic freedom. To control errors when extrapolating to the continuum a → 0, one should simulate in the scaling regime of a m −1 IR . We denote the onset of the scaling regime by a s , and β s . For In the important case of SU (3) with β s = 6, there are five crystal-like subgroups with V as the largest. For all subgroups, β f < β s , with V having β f = 3.935(5) [28] and thus inadequate to reach the scaling regime. Other work [49] has shown that extending to a subset including the midpoints between elements of V raises β f ≈ 7. However this require more qubits and -potentially more worrisome-sacrifices gauge symmetry completely which is potentially dangerous on quantum computers [22,50,51].
To increase β f , attempts to supplementing the Wilson action with additional terms were made [28,36,37,40,43,[52][53][54][55], although only in [28,40] were Monte Carlo calculations undertaken for SU (3). Two reasons suggest this would help. First, additional terms which have a continuum limit ∝ Re Tr F µν F µν , but take different values on the element of H (e.g. Tr(u p ) and | Tr(u p )| 2 − 1), can change the gap ∆S and thus a f . Second, new terms can reduce finite-a errors as in Symanzik improvement.
The term usually added was the adjoint trace, giving (2) where u p ∈ V, and the first term is normalized so for β {1,−1} = 0, the S[u] matches the Wilson action (with β {1} = β). In these works, no relationship was assumed between β {1} and β {1,−1} . That Eq. (2) improves the viability of V over the Wilson action will be shown in [56].
For a different action, smaller values of a f were demonstrated in [28]. Comparing to Eq. (2) which, up to a constant, adds an additional group character to the action, the second term in Eq. (3) is only a character for Abelian groups. Later, we will rewrite this action in terms of characters of non-Abelian groups.
With these actions, the dimensionless product T c √ t 0 of the pseudocritical temperature and the Wilson flow parameter were found to agree in the continuum with SU (3), allowing one to set the scale of those calculations. a > 0.08 fm was achieved without the effects of a f being seen. This suggest that V can reproduce SU (3) in the scaling region with a modified action, such that practical quantum computations of SU (3) could be performed. While promising, the choice of new terms was ad-hoc and left unclear how to systematically improve it or analyze relative effectiveness. In the next section, we discuss how one can systematically derive lattice actions, discovering that the terms added in these two actions are in fact the first terms generated in a cumulant expansion.

III. GROUP SPACE DECIMATION
Our ultimate goal is to approximate the path integral of group G faithfully by a discrete subgroup H by replacing integration over G by a summation over H. Group space decimation can be understood in analogy to Wilsonian renormalization, where we integrate out continuous field fluctuations instead of UV modes. The typical method used with discrete subgroup approximations is to replace the gauge links U ∈ G by u ∈ H such that the action S[U ] → S [u]. This corresponds to simply regularizing a field theory. For strong coupling, this appears sufficient. As β → ∞, correlations between gauge links increase and the average field fluctuation becomes smaller. When the average field fluctuations decrease below the distance between 1 and N of the discrete group, freeze-out occurs and the approximation breaks down-similar to probing a regulated theory too close to the cutoff. Therefore, improving this approximation and understand the systematics can be done by considering these discarded continuous field fluctuations. To do this, instead performing the replacement U → u, we will integrate out the continuous fluctuations, following the decimation formalism developed by Flyvbjerg [36,37]. He derived the second order decimated action for U (1), SU (2), and SU (3). An important general feature of the decimated action though is missing from this second order action -while new terms are generated at every order, it isn't until third order that the coefficient of existing terms are modified. One would expect such terms are critical to understanding deviations from the continuous group and therefore we compute them for the first time in Sec. IV. . N for H are given by red points. We have applied the S 2 metric to obtain the Ω. In groups representable in two dimensions, this region resembles a polygon while in higher dimensions, it becomes a polytope.
It is natural to associate every subgroup element u ∈ H with an unique set, or region, Ω u containing all closest continuous group elements U ∈ G: where the distance is defined as d 2 (U, u) = Tr (U − u) † (U − u) . By such a definition. the continuous group is fully covered, i.e., G = ∪ u∈H Ω u and a graphical demonstration of Ω ≡ Ω 1 can be found in Fig. 1. Note that for any U ∈ G, there exist a unique u ∈ H and ∈ Ω such that U = u , where we may treat as the error of u approximating U . In this way, without approximation, the Euclidean path integral integrating over G can be written as a summation over H and integration over ∈ Ω: where Z is a functional integral over all gauge links U on the lattice, or equivalently a functional integral over and a functional sum over u. In this expression, is defined by replacing each gauge link U by u . We then expand the exponential in the path integral and integrate over producing a moment expansion where we have introduced the notation f = Ω D f with normalization Ω D = 1. What we are really after is an expansion for the action S[U ], writing Z in terms of a cumulant expansion allows us to match Eq. (6) with (7) to obtain an effective action. In this way, after integrating over , the contributions to the action depend only on the discrete group gauge links u i.e.
One may worry about poor convergence in the region of interest β ≥ β s ≥ 1 . As will be discussed more thoroughly in Sec. V, β n terms are suppressed by powers of the average field fluctuation. Thus, the size of the discrete group, which determines the size of field fluctuations integrated out, also determines the series convergence. Starting with the second order terms computed in Refs. [36,37], the decimated action generates multiplaquette contributions. Their inclusion in quantum simulations brings substantial non-locality which requires high qubit connectivity and increases circuit depth. Luckily these contributions will be shown to be small in Sec. IV.
In the following, we calculate Eq. (8) to (10) in terms of linear combination of the group characters starting from the Wilson action of Eq. (1): Here we introduced χ r , the character of the group representation 2 r. This is the natural basis for the decimated action. All characters required for our β 3 calculation are in Table I. In the interest of deriving a decimated action for general gauge groups, we have chosen a nonstandard basis for U (1) and SU (2). This allows for one general scheme for U (N ) and SU (N ) groups. This basis is not linearly independent and relations between representations exist. This dependence is typically used to write U (1) and SU (2) in reduced sets of representations. For U (1), the resulting identities are For SU (2), one finds that and for SU (3), the set of dependent representations needed up to third order in the cumulant expansion are In deriving the decimated action, integrating out the field fluctuations require us to reduce expressions of the form i1j1 · · · injn . To simplify these, we use an identity derived in [57] for SU (N ) and U (N ) groups that for any where ε is Levi-Civita symbol, A j i are the contracted dummy indices, and B n is the Bell number accounting for the number of ways that one can put the open indices i k , j l on ε such that no i k and j l appear in the same ε. In [57], Eq. (15) was derived for integrating over the entire group G. Hence in our case, we need to determine the constants c i 's when integrating only over Ω for . This is done by contracting the tensor structure on each side of Eq. (15) with products of Kronecker delta's and solving the resulting linear equations. What we are left with are expectation values of χ r over Ω where d r is the dimension of representation r. At first order, only one integral is needed: At second order, there are two relations At third order, there are four structures, but by complex conjugation one can reduce this to two unique ones: For U (1) → Z n , there is only one representation at each order of the cumulant expansion, V {h} = h . These terms can be computed analytically by a change of variables = e iφ [37]: with h = 1, 2, . . . being integers and the normalization Extending this to non-abelian groups, e.g. SU (N ), Ω becomes a high-dimensional polytope. In [36], the V r for BI and V were computed up to second order by approximating these polytopes with hyperspheres to two significant figures. It is crucial to remove these approximations for our purpose because the uncertainty δV r ∼ O(1%) is magnified in the coupling constants of the decimated action. These couplings are combinations of powers of V r with extreme cancellations making the fraction errors grow rapidly. Hence we avoid making the hypersphere approximation and numerically compute all the V r necessary for the 3rd order actions to 4 significant figures. (Results found in Table I.)

IV. ORDER-BY-ORDER DECIMATION
In this section, we will undertake the task of computing the decimated action order-by-order. The first order is relatively straight-forward, and only contains a single plaquette term. Working from Eq. (8) After applying Eq. (17), the remaining S 1 [u] is found to depend only on u: where β (n) r is the n-th order contribution to the coefficient of 1 dr Re χ r . It is comforting that at first order, no new terms are generated in the decimated action. This allows for rescaling β (1) {1} → β, recovering the regularization procedure of simply replacing U → u in the Wilson action. Although this rescaling is permitted, V {1} contains content about the approximation of G by a given H. As the number of elements of H increases, Ω should shrink and V {1} → 1 from below. This means V {1} quantifies how densely H covers G and thus how small field fluctuations can be before freeze-out occurs. Since β (1) {1} β, these decreases in V {1} signals the poorness of approximating using Eq. (24) alone. This is discussed further in Sec. V.
We now proceed to calculate the second order decimated action while fixing a few typos in [36] along the way. The second order decimated action c can be made into three terms based on how the two plaquettes p and q are related: p = q (oneplaquette contribution), p ∩ q = 1-link (two-plaquette contribution), and p ∩ q = 0-links. To all orders, the p ∩ q = 0 contributions to the decimated action vanish. The first term in Eq. (9) for case p = q reads: where we have utilized Eqs. (18) and (19) to contract the u's after integration. The second term of Eq. (9) is obtained from first order action of Eq. (24) which reads, where we have used the identity Hence we conclude that the second order effective action term for two identical plaquettes is where the β (2) r can be found in Table II. Next, we calculate the case of p ∩ q = 1-link for the second order decimation. Contracting the δ's in Eqs. (18) and (19), where unlike Eq. (25), we only identify one link as the same between the two plaquettes. This leads to the following expression, where we have used the fact that all the V r 's are real due to our choice of the integration region. The explicit expressions for the couplings are found in Table II. Note that this expression is also applicable to U (1).
contributes to the action. It would be desirable if these terms could be neglected, because they require substantial quantum resources to implement. By inspecting Table III, one observes that the two-plaquette β r are O(0.1) or smaller than the single-plaquette terms. The largest coupling, β 2i , multiples a term Im χ 1 Im χ 1 ≈ 0. Strong cancellations are expected from correlations between the remaining terms (shown in Fig. (2)) as evident by the observation β 2t ≈ −β 2u .
It is reasonable to expect these individual reasons to persist at higher orders, suggesting that at a fixed order all multi-plaquette terms can be neglected compared to their 1−plaquette counterpart. But can we argue that the multi-plaquette terms generated at order O(β n−1 ) are still negligible when the O(β n ) contribution is introduced? To do this, we look at the continuum limit of each term being introduced. In this way, we recognize that the two-plaquette terms are related to the Lüscher-Weisz action [59]. k-plaquette terms corresponds to applying 2k − 2 derivatives to a 4 Tr F F and are thus O(a 2k+2 ).
Here F is the field strength tensor projected onto the lattice directions. Combining this with the observation that for a coupling β j generated at O(β n ) has the scaling β j ≈ 10 −n β n , we estimate that where D is a covariant derivative projected onto the lattice directions. The combination of higher powers of a and their associated expectation values of higher-dimensional operators of slow-varying fields should be sufficient to suppress the mild enhancement from the couplings, at least for third-order decimated actions. For these reasons, we will neglect higher order multi-plaquette terms in this work.
For the third-order terms of Eq. (10), we, therefore, only focus on the case where three plaquettes are identical. This will be done term by term, where the first term is: Re Tr(u 5 5 u 6 6 (u 7 7 ) † (u 8 8 ) † ) Re Tr(u 9 9 u 10 10 (u 11 11 ) † (u 12 12 For the mixed-order term in Eq. (10): where the second line was simplified with the identities: The final term in Eq. (10) follows from another identity: Combining Eqs. (31), (32), and (36) we arrive at the third order contribution to the single-plaquette decimated where the overall factor of 1/3! has been absorbed into the definition of β (3) r . Note that, unlike the second order results where only certain decimation programs generate renormalization for existing terms, the third order S 3 [u] introduces corrections to Re χ 1 for all G → H. Additionally, a number of the specific group identities in Eqs. (12)-(14) also lead to renormalization.
Putting together Eqs. (24), (28), and (37), the decimated action of Eq. (7) to the third order for a general gauge group including only the single plaquette terms read, where the β r are collected in Table II. Note that this decimated action is correct for any G → H with its associated V r . Referring back to Eqs. (12), (13), and (14), for a given gauge group simplifications occur. We write S[u] for 3 groups of prime interest, U (1), SU (2), and For SU (2): For SU (3):

V. FINITE GROUP EFFECTS
With Eq. (38), it is possible for us to investigate systematically the effect of replacing the continuous group by its finite subgroup. In order to proceed, it is useful to introduce a new parameter which approximately represents the field fluctuations. To do this, consider the representation of a continuous group lattice gauge link in terms of the corresponding generators λ a in the adjoint representation, U = e iλaAa , where a summation over color indices a is implied. In this form, we see that the gauge fields correspond to amplitudes in each of the generators. For ∈ Ω, inserting its small parameter expansion where DA is a measure over all A a which respects gauge symmetry and c (n) r are representation and groupdependent constants. From this, we see that as the subgroup H incorporates more elements, Ω → 0 and V r → 1 from below. This means that for finite Ω the domain size of A a in Ω is an indicator for deviations from G of H. Flyvbjerg defines a parameter R as the radius of a hypersphere with equal volume to Ω to get a handle on domain of A a . This allows him to approximate V r analytically [36,37]. Here, we can use this idea to roughly understand the scaling of V r .
For U (1) → Z n , the hypersphere is exactly Ω and R cleanly defines ≤ R = π/n. Beyond U (1), the connection between Ω and a single value of R is complicated because the Ω of H form polytopes in the hypervolume of their continuous partner. In this case, while Ω is contained by a hypersphere centered at 1 whose boundary incorporates elements in Ω, some element of the hypersphere are not included in Ω. On the other hand, there exists a largest hypersphere centered at 1 that only contains elements in Ω. In this way, we define an upper and lower bound for R. Note, this is different from [36,37] where the polytopes of H were always approximated by hyperspheres with definite radii. For SU(2) with BI, we find 0.09 ≤ R 2 ≤ 0.15 which can be compared to R 2 sphere = 0.12 of [36,37]. In the case of SU(3) with V, 0.42 ≤ R 2 ≤ 0.93 compared to R 2 sphere = 0.62. While superficially the cumulant expansion has appeared as a strong-coupling expansion in β, the actual behavior is controlled by both β and R with R controlling V r . As pointed out in [36,37], the leading order behavior for small R for a given power of β α is actually O([βR 2 ] α R −2 ). Therefore one would predict that the relative smallness of R 2 for BI compare to V signals that β f should be larger for BI which is indeed the case.
For the subgroups of SU (3), this scaling behavior be-  [36,37] at second order.
0.13314β 2 0.0046477β 2 0.0031836β 2 comes unsatisfactory because R 2 ∼ 1. It is possible to study this breakdown in U (1) → Z n where the systematic effect of decimation can be studied in detail both because errors can be made arbitrarily small for large n and because V r and β r are known analytically. In terms of R, one can expand the β r for the U (1) action of Eq. (39) to find: The first thing to note is that the O([βR 2 ] α R −2 ) scaling found in [36,37] continues to the third order. One might be tempted to use this leading behavior to estimate the β f or the radius of convergence of this series, but this would be incorrect. Instead, it would behoove one to notice that for both 2nd and 3rd order contributions, the expansion coefficients to the subleading terms [R 2 ] k R −2 with k > α initially grow until a 1/k! factor dominates over all the other factors.
But what is the origin of this behavior? For simplicity, we can understand this behavior by considering the expansion of V 4m j which form β r . The specific combination of V 4m j dictated by the cumulant expansion ensures that orders lower than O([βR 2 ] α R −2 ) cancel in β r . The j representation contributes to β {r1,··· ,r k } in the form of V 4m1 n } under the constraint |r| = m i j i where |r| = |r 1 | + · · · + |r k | and j i = |j i 1 | + · · · + |j i l |. One might worry that studying the expansion of V 4m j isn't representative, but one can verify that the scaling behavior observed below persists in β r , although the numerical factors become cumbersome. For these V 4m j , we have from which, we see that the coefficients of the [R 2 ] k R −2 contributions to β r are accompanied by a factor ∝ 1 k! m k−1 j 2k−2 . While the factorials ensure the series converges, β r for higher representations r have larger m, j, or both leading to higher order terms in the expansion being large for moderate R. This helps explaining why Z 4 with the Wilson action fails to replicate U (1) substantially above β = 1 -while the naive scaling would suggest R √ β k 1−k would be enough to suppress higher representations, in reality a stronger bound of max{ 1 k! m k−1 ( jR) 2k−2 } 1≤ j≤|r| 1 for ∀β r is required for all subleading terms to be small. Considering the range of m with fixed |r|, j, the bound is strictest when |r| = j yielding R 1/ j in order for the lowest order contribution to dominate such that R √ β α 1−α provides a reasonable estimate for the range of β where the decimated action provides a reasonable approximation for its continuous partner. While these conditions are satisfied for BI, they Another feature observed in the R expansion of the Z n group is that because V r ∝ sin rR, the sign of the O([rR] k ) terms oscillate, and therefore the sign of β (n) r can depend sensitively on R. Since Re χ r (1) > Re χ r (N ), the overall sign of β (n) r determines whether or not the r-th term in the action enters the frozen phase in the limit of β → ∞. This behavior is observed in Fig. 5 where β (1,2) From the behavior observed in U (1), we can improve the quantitative understanding of how well H can approximate G, even when β r are not known analytically. Clearly, V r → 1 indicates that the R → 0, and in that limit the two actions would agree. Therefore, the difference between the two actions serve as an indicator of β f . This proxy can be compared to others in the literature. The simplest estimate is β −1 f ∝ ∆S = Re Tr(1) − Re Tr(N ) [45]. While this estimate finds monotonic behavior for discrete groups of U (1) and SU (2), different O(1) factors are needed. It also fails completely for SU (3), as seen in the left panel of Fig. 3.
Observing the differing O(1) factors, [45] suggested a different estimate. For discrete Non-abelian subgroups near β f , S[u] is dominated by contributions from u p = N . From duality arguments , the action near β f could be approximately rewritten as a Z C action where C is the minimal cycle such that u C = 1 for all u ⊂ N . Since C = n for Z n , these arguments predict a single curve β f ≈ 0.78/(1 − cos(2π/C)) directly from the study of β f in Z n for all discrete subgroups. The discrepancy between SU (2) and U (1) was reduced from ∼ 300% to ∼ 50%. The authors of [45] warned that this approximation could be poor for SU (3) albeit without numerical evidence.
Since then β f for the subgroups of SU (3) have been found and as anticipated, this estimator proves to be poor as presented in the center of Fig. 3. In the plot on the right of Fig. 3, β f is plotted as a function of (1 − V 4 {1} ) −1 . We find that monotonic, linear behavior is observed within the uncertainties for each continuous group. Best fit lines have been included for each group to guide the eye. This suggests that our estimator captures some of the non-perturbative behavior near the freezing transition better than ∆S −1 or C. Physics of the different groups differ, as signaled by their different scaling regimes.
If we divide β f by a rough estimate of β s = [1, 2.2, 6] for U (1), SU (2), SU (3) respectively, we might expect to further remove some of this group dependence. Doing so in Fig. 4, we find that SU (2) and SU (3) collapse onto a single line and U (1) within 25%.
Using our higher order results, one can then gain insight into the effectiveness of the ad-hoc actions of V. Each of these actions corresponds to terms that are generated at 2nd order in the decimated action. The first ad-hoc action used in [28] can be rewritten as The trajectory parameters were chosen to be parallel to the freezing point at large β 0 by eye. From Fig. 5, we see that in both ad-hoc actions, reasonably agreement is found for intermediate β for the 3rd order action. Here β is the coefficient in front of Re χ {1} for the ad-hoc actions. The ad-hoc trajectories are known to poorly reflect G at low β, because they lack curvature to fix the known requirements at β = 0. At large β, we expect higher order terms in the cumulant expansion to become relevant and thus disagreement is expected to occur. This surprising agreement in the intermediate region of β suggests that actions formed by neglecting terms in the cumulant expansion are optimized in their character basis by setting the couplings to results given by the resummed cumulant expansion.

VI. CONCLUSION
In this work, we used the cumulant expansion to develop a systematic method for studying and improving lattice actions that replace continuous gauge groups by their discrete subgroups. This is a step in the ongoing trek toward developing accurate and efficient digitization on quantum computers. These decimated actions, through the factor V {1} , have superior predictive power for finding the freezing transition compared to prior estimators. We further computed the third-order, single-plaquette contribution for the general group. These higher-order terms are necessary for systematizing the decimation procedure of SU (3) → V where it has been observed that the inclusion of terms generated in the second-order cumulant expansion with ad-hoc couplings improve the approximation of SU (3). The most immediate work in these directions would be a classical simulation of the full decimated action of Eq. 38 and compare it's effectiveness to the results of [28,56]. Given the large corrections from second to third order for V, additional work should be devoted to computing the fourth-order contributions.
Another important step in studying the feasibility of this procedure is to explicitly construct the quantum registers and primitive gatesà la [60] where smaller discrete groups were investigated. Together with classical lattice results, this would allows for resource counts.