Towards the continuum limit of a $(1+1)$d quantum link Schwinger model

The solution of gauge theories is one of the most promising applications of quantum technologies. Here, we discuss the approach to the continuum limit for $U(1)$ gauge theories regularized via finite-dimensional Hilbert spaces of quantum spin-$S$ operators, known as quantum link models. For quantum electrodynamics (QED) in one spatial dimension, we numerically demonstrate the continuum limit by extrapolating the ground state energy, the scalar, and the vector meson masses to large spin lengths $S$, large volume $N$, and vanishing lattice spacing $a$. By exactly solving Gauss' law for arbitrary $S$, we obtain a generalized PXP spin model and count the physical Hilbert space dimension analytically. This allows us to quantify the required resources for reliable extrapolations to the continuum limit on quantum devices. We use a functional integral approach to relate the model with large values of half-integer spins to the physics at topological angle $\Theta=\pi$. Our findings indicate that quantum devices will in the foreseeable future be able to quantitatively probe the QED regime with quantum link models.

Introduction -The rapid development of quantum technologies culminating in the precise control of large quantum systems [1][2][3][4] has fundamentally altered the scope of physics questions that can be addressed for strongly interacting systems. As a complementary approach to smashing nuclei in colliders to uncover their substructure, one can realize quantum many-body systems in controlled analog quantum simulators or digital quantum computers, and study their ground state properties at finite density or during real-time evolution associated with quenches, which are extremely difficult to tackle using Markov Chain Monte Carlo methods [5,6].Motivated by this possibility, pioneering proposals [7][8][9][10] have been put forward to study properties of lattice gauge theories with the long-term goal of simulating quantum chromodynamics (QCD), the theory of strong interactions. These ideas have triggered an intensive theory effort to devise efficient and feasible implementations (for recent reviews see, e.g., [11][12][13]), resulting in first experimental realizations [14][15][16][17][18][19][20][21] in recent years.
Nonperturbative calculations of quantum field theories (QFTs) require a careful treatment of regularization and renormalization, for which the lattice approach has proven most successful [23]. The lattice Hamiltonian of a gauge theory [24] retains exact gauge invariance while a finite spatial lattice reduces the infinite number of degrees of freedom of the field theory to a finite number of lattice sites and links. The local Hilbert space dimension of the gauge fields, however, remains infinite in the original Wilsonian formulation. Quantum link mod-els (QLMs) [25][26][27] regulate these infinite dimensional Hilbert spaces with qudits while maintaining exact gauge invariance. They are ideal candidates to be studied on quantum devices such as analog quantum simulators or digital quantum computers, which typically work with finite dimensional local Hilbert spaces. Extracting information relevant for the QFT in the continuum limit, especially when realized in the low-dimensional Hilbert spaces available in current quantum devices, requires a sequence of extrapolations which is the main topic of this article.
The continuum limit physics is quite sensitive to the nature of the employed truncation. For example, in the (1 + 1)−d O(3) model the physics of asymptotic freedom could only be recovered with at least a 16-dimensional local Hilbert space in the angular momentum basis truncation [52]. However, using qubit operators [53], it was shown that the same continuum limit only required 2qubits per site [54]. The QLM approach is similar to this qubit-regularization, but uses larger spin-S opera- Figure 1. (a) In a U(1) QLM, matter fields (blue dots) reside on the sites of a lattice while gauge fields, represented by spins (red arrows), live on the links connecting two neighboring lattice sites. Gauss' law ties consecutive gauge fields to the matter in between as indicated by the shaded ellipses. The ground state energy, shown in (c), and the first two excited states (with vector and scalar quantum numbers), shown in (b), in the zero momentum sector of QED obtained from the U (1) QLM using ED and iMPS show excellent agreement with the analytical prediction for m/e = 0. In (b), the grey solid lines indicate the leading order (LO) analytical expansions [22] for small m/e and e/m, respectively. The error bars indicate an estimated systematic uncertainty (see Section IV in SM). tors for the U (1) gauge links. It was analytically shown that large representations for the gauge links in QLMs recover the standard Wilson lattice gauge theory [55]. A fine-tuning free approach to the continuum limit using the QLMs is via the dimensional reduction in the Dtheory formulation [56]. This is, however, only possible if a phase with an exponentially large correlation length is generated [57]. In this article, we explore the former approach and show that careful analysis techniques allow us to reach the continuum limit quantitatively for certain physical observables, even with very small values of S ( 3 − 5), extending the observations made in [42,43]. We solve the Gauss law analytically for a general spin-S representation, and estimate the quantum resources to simulate the continuum limit. Our results may help to improve digital quantum simulations of the massive Schwinger model [58][59][60], through vari-ational quantum eigensolvers and Trotterized real-time dynamics. Finally, we use path integrals to demonstrate how large half-integer spins give rise to the topological angle Θ = π.
Hamiltonian and Gauss' law -We focus on U (1) gauge theories in d spatial dimensions with (staggered) fermionic matter [24]. The gauge fields are described byŜ n,j = (Ŝ + ,Ŝ − ,Ŝ z ) n,j , the raising, lowering and zcomponent quantum spin-S operators on the links (n, j) connecting neighboring sites n and n + e j on a hypercubic lattice of size N d 1 . For any S, the raising (lowering) operators are S ± = (S x ± iS y )/2 with S x/y the x/y component of the spin vector S. For the spin-1/2 representation, for example, the spin operators are related to Pauli matrices, S = σ/2. The spins are coupled to fermionic operatorsψ n on the sites, as described by the Hamiltonian The first two terms are the electric field energy at the bare coupling g, and the staggered fermion mass µ. The gauge-matter interaction is the correlated hopping of fermions along a link with the simultaneous raising or lowering of the corresponding spin. For d > 1, there is the magnetic energy term,Ĥ m = − P4 + h.c. , with P labeling elementary plaquettes consisting of links P 1,2,3,4 forming a square. We identify the spins with gauge fields (Û ,Û † ) ↔ (Ŝ + ,Ŝ − )/ S(S + 1) and electric fieldÊ ↔ S z . This identification preserves the commutation relations, [E, U ( †) ] = (−)U ( †) , as well as an exact gauge symmetry with (2S + 1)-dimensional Hilbert space. To achieve the correct scaling behavior, appropriate factors of S are inserted in the dimensionless couplings in Eq. (1) (see Sec I in SM for more details). The gauge transformations are generated by the Gauss law operator, satisfying Ĥ ,Ĝ n = 0. The Hilbert space thus separates into superselection sectors labelled by eigenvalues ofĜ n . For the physical Hilbert space, H phys , we requirê G n |phys = 0. In the renormalization group (RG) sense, the parameters g, µ, and S can be regarded as directions in the space of couplings, to be adjusted such that the theory flows to a fixed point corresponding to the desired QFT. In the remainder of this letter, we focus on the case of one spatial dimension, d = 1, where Eq. (1) provides a lattice version of the (massive) Schwinger model [61][62][63]. The continuum Schwinger model is parameterized by the (bare) values of the electric charge e and the fermion mass m. In its lattice version at lattice spacing a, these parameters appear through the dimensionless combinations g = ae and µ = am. The QFT limit is reached for large S, large N , and small a, as demonstrated in Fig. 1. Specifically, we first take the infinite spin length limit S → ∞ at fixed a and N ; then the thermodynamic limit N → ∞ at fixed a; and finally the continuum limit a → 0 at fixed µ/g. The different extrapolations have to be performed for appropriately rescaled ("renormalized") quantities that correspond to physical observables (see Sec IV in the SM for details). In general, the final continuum limit involves a rescaling of the dimensionless coupling constants in order to reach the RG fixed point [23]. This complication is absent for the present model (except for a redefinition of the ground state energy) [62,64].
Mass spectrum from the QLM -The massive Schwinger model is considerably simpler than QED in higher dimensions due to the absence of magnetic interactions and the strong Gauss' law constraints. Consequently, both the weak and strong coupling limits, e/m = 0 and e/m = ∞, respectively, are exactly solvable, and analytic expansions around these limits [22,65] can be used for benchmarking the extrapolation. Our principal numerical methods are exact diagonlization (ED) (using the Python package QuSpin [66]) and variational techniques based on infinite Matrix Product States (iMPS) [67,68] directly in the thermodynamic limit. To perform ED, we derive an equivalent spin model constrained by a projector P on neighboring gauge link configurations allowed by Gauss' law. The resulting Hamiltonian has the form (see Sec II in the SM): For S = 1/2, this reduces to a constrained model of hardcore bosons [69], sometimes referred to as PXP model, whose total number of allowed states can be analytically counted [70], and which is believed to explain the anomalous thermalization observed in the 51-Rydberg atom experiment [71]. The relation between a PXP model and the spin-1/2 QLM was first noted in [72], which we generalize here to arbitrary S. Moreover, we extend the analysis of [70], and derive an analytic expression for the dimension of the physical Hilbert subspace, given by (see The remarkably small Hilbert space size, scaling linearly with S at a fixed N (see Sec. II in the SM for an illustration), enables ED calculations for relatively large system sizes, and data up to S = 3 and N = 16 is presented here. Results from iMPS simulations in the thermodynamic limit are also shown for S ≤ 5 for the Hamiltonian in Eq. (1), where Gauss' law is enforced by adding a large energy penalty ∝ nĜ 2 n [73,74], withĜ n defined in Eq. (2). Our results, summarized in Fig. 1, demonstrate how to accurately reach the continuum limit with QLMs for the ground state energy and the energies of the first two excited states, the "vector" and "scalar" particles. We find excellent agreement in the strong coupling limit, due to small fluctuations around Ŝ z = 0, such that our numerics with small S (≤ 3 for ED and ≤ 5 for iMPS) already capture the relevant physics. For higher excited states or towards weak coupling, the fluctuations grow more pronounced and larger spin lengths S become necessary 2 .
The detailed steps of the underlying extrapolation are shown in Fig. 2 for ED, illustrated for the analytically solvable strong-coupling limit (m/e = 0). Despite small spin lengths S = 1, 2, 3, we observe a clear 1/S scaling, enabling a reliable extrapolation to the S → ∞ limit. Similarly, the subsequent N → ∞ extrapolation is performed with the expected leading behavior at large N . The largest systematic error arises from the choice of fit range for the final a → 0 extrapolation, which require increasingly large values of N and S, attributed to increasing electric field fluctuations at the continuum limit. Empirically, we find that systematic errors are minimized by disregarding "far-off" N and S extrapolation where the extrapolated values differ by more than 10% from the one of the largest available system size. We thus select a smallest lattice spacing a for which the underlying data is sufficiently converged with respect to S and N . Note that this procedure naturally depends on the observable. The ground state energy can be extrapolated with lattice spacings down to ae ∼ 0.1, while we only reach ae ∼ 0.3 for the scalar mass. Details of the numerical extrapolations are presented in Sec IV of the SM. Figure 3 shows our final results for the vector and scalar masses with a trivial m/e dependence subtracted. The iMPS results are obtained analogously to the ED simulations, but without the N extrapolation. The  Figure 2. We illustrate the sequence of extrapolations in the panels from the left to the right, required to reach the continuum limit for the ED data with m/e = 0. The energies of the vacuum (bottom row), vector particle (middle row), and scalar particle (top row) are extrapolated to S → ∞ (left column), N → ∞ (middle column), and a → 0 (right column), as discussed in the main text. The circles in the middle row indicate the values obtained from the corresponding S-extrapolations shown in the left column. Similarly, the ticks in the right column indicate the values corresponding to the N -extrapolations in the middle column. For clarity, we only show selected values of ae and N for the first two extrapolations. For comparison, the green crosses indicate the exact analytical results. All fits employed for the extrapolations are polynomials and χ 2 denotes the resulting normalized square error of the fit (see Sec IV of the SM for details, where we also provide a table with our quantitative results). agreement between both approaches demonstrates that the thermodynamic limit is reached, and that the limits S → ∞ and N → ∞ commute for this model. Comparing to previously obtained results [75] at S → ∞, we find good agreement of both the ED and iMPS data for the vector mass, indicating that S = 3 is sufficient to resolve this excitation. As anticipated, the scalar mass requires larger S values, and we observe stronger deviations in ED.
Estimation of required resources on a quantum device -As illustrated above, already very small spin lengths S 3 and system sizes N 16 are sufficient to obtain quantitative estimates for the low-lying mass spectrum. A brute-force implementation would nevertheless still require controlling a Hilbert space of dimension [2(2S + 1)] N ∼ 14 16 ∼ 2 61 . Due to Gauss' law, most of these states are unnecessary.
To estimate the minimal required resources to implement the model on a quantum device naturally working with the corresponding qudits of size 2S + 1 [76][77][78][79][80][81][82][83], consider the equivalent spin model, Eq. (3). Then control over only N ∼ 16 such qudits would be sufficient to reach the continuum limit. According to Eq. (4), a perfect encoding on a digital quantum computer would need only dim H phys (S = 3, N = 16) = 63757 < 2 16 states, enabling our procedure to be carried out on existing quantum computing devices with control over 16 qubits. This fact has been already exploited to carry out the ED calculations on a conventional laptop computer.
To illustrate the applicability of our improved encoding [Eq. (4)], consider two examples: (i) the mass-spectrum discussed above using a perfect qubit-encoding, and (ii) real-time dynamics on qudit hardware. The first example may be tackled with a variational quantum eigensolver (VQE) [20,84] which finds an optimal representation of an input variational ansatz for the ground or a lowlying excited state, by classically minimizing Ĥ . We emphasize that this does not require to actually imple-mentĤ, but only to measure its expectation value. It is then favourable to work in a computational basis with the perfect encoding whereŜ z is diagonal, such that most terms ofĤ can be measured directly. OnlyŜ x has to be treated separately, which we leave for future work, but we note that it will remain local because the projection P acts locally. Having obtained the spectrum via VQE, the masses are extracted by classical post-processing via the extrapolations discussed above.
The first two single-qudit gates can be efficiently executed in parallel. For a fixed n the off-diagonal term witĥ S x n involves a projector acting locally on (n − 1, n, n + 1). These corresponding Trotter evolution can therefore be realized as a three-qudit controlled unitary acting only on the middle qudit n if the triplet (n − 1, n, n + 1) is compatible with Gauss' law. This suggests to parallelize this last term in three layers, and thus scalable to large system sizes. This approach might be used, e.g., to observe the type of quench dynamics discussed in [85], where the required small spin lengths S ≤ 4 are within reach of exisiting technology (see, e.g., [86] for a review of qudit quantum computing and [87] for an experimental realization, as well as [88][89][90][91][92] for recent related developments using qudits.).
Field theoretic description -To address the closely related question of covergence to the Kogut-Susskind 3 Note that here we do not enforce the perfect encoding. The Trotter decomposition is then valid on states |ψ that fulfill P|ψ = |ψ , and since P commutes withŜ z n , the projector only appears with the terms involvingŜ x n . limit, we use coherent state path integrals. Physically, the worldline of a spin-S at spatial site n traces out an arbitrary closed curve on the Bloch sphere under the Hamiltonian evolution in imaginary time, subtending a solid angle Ω, as shown in Fig. 4(a). The path integral Z = TrĜ e −βĤ , (see Sec V of the SM) is dominated by the electric field term, Ŝ z n 2 , which contributes as ∼ exp − g 2 S(S − 1/2) cos 2 (θ n )/2 , where is the Trotter discretization and θ n , φ n are the angular coordinates of the spin. For large g 2 S, θ n → π/2 as indicated in Fig. 4(b). In this limit, the quantum spin is confined to the equator of the Bloch sphere and reducse to the quantum rotor of the Wilson-Kogut-Susskind formulation, consistent with the numerics. An important additional observation is that the area traced out by the closed curve on the Bloch sphere byŜ n admits a topological interpretation as a Berry phase. In the continuum limit, this area is ω[Ω] = β 0 dτφ(τ ) cos(θ(φ(τ ))) = φ0 φ0 dφ cos(θ(φ)) [93], where τ ∈ [0, β] is the Euclidean time. At leading order in the large-S limit, only the number of total windings is relevant, with each winding contributing 2πS. While this term is irrelevant for the integer spins used in this work, for half-integer spins it gives rise to a π-flux in the system. This feature rigorously establishes that QLMs with large values of half-integer spin lead to a topological theta angle Θ = π, often anticipated in the literature [10,19,72,94] already for the spin-1/2 case.
Conclusion and outlook -In this paper, we have numerically demonstrated the continuum limit of QED in (1+1)−d, regularized with quantum spin-S operators for several physical observables. A highlight of our results are the small spin values S 3 − 5 that suffice to reach the continuum limit. The systematic finite size scaling enables us to quantitatively estimate the resources of a quantum device to realize the continuum limit. These results lend hope that near future quantum simulation experiments with small S and limited lattice size can yield valuable data that can be extrapolated to the QFT limit. Using a coherent state path integral, the approach to the Kogut-Susskind limit is derived and the connection of half-integer spins with topological angle π is formalized. For future investigations, the feasibility of our approach may immediately be tested using today's quantum hardware. High-energy aspects of these models, such as the presence of quantum scars [95], and Floquet dynamics [96,97] are worth studying. In higher dimensions, an exciting challenge is to quantify the convergence properties of both Abelian and non-Abelian QLMs with increasing link representations. Finally, we note that the calculation of the mass spectrum of the Schwinger model using quantum devices is merely a stepping stone towards more complex tasks. The situation changes completely when considering, e.g., real-time dynamics where our proposed Trotterization using qudits can provide a substantial advantage of traditional qubit approaches.

I. SPIN LENGTH AND LATTICE UNITS IN THE QLM
Throughout the main text, we work with units where the speed of light and the Planck constant are set to one, c = = 1. In standard lattice QED, all quantities are rescaled with the lattice spacing a to become dimensionless. For example, the dimensionless coupling is given by g 2 = e 2 a 3−d , where e is the bare value of the electric charge in d spatial dimensions, and µ = am with the bare mass m. Similarly, the Hamiltonian is rescaled aŝ H lat = aĤ phys and the fermion operators are related bŷ phys , such that ψ lat n , ψ lat † m = δ n,m .
In the QLM, dimensionless gauge fields are further replaced by spin operators. Explicitly, the electric field E lat = a d−1 /eÊ phys and the compact U (1) link operator U lat = e iaeÂ phys are replaced according tô Here, and in the following the operatorsŜ ±,z fulfill the angular momentum algebra [Ŝ z ,Ŝ ± ] = ±Ŝ ± and [Ŝ + ,Ŝ − ] = 2Ŝ z with fixed spin length S. These identifications are consistent with the canonical commutation relations in the continuum, Â phys (x),Ê phys (y) = iδ(x − y). Here, the fields with subscript "phys" serve as a reminder of the corresponding (dimensional) fields in the continuum theory. As a cross-check, naively sending a → 0 everywhere reproduces the desired continuum expressions. Crucially, the explicit appearance of the spin length S ensures the lattice commutation relations Ê lat ,Û lat =Û lat and Û lat ,Û † lat = 0 in the limit S → ∞. To see this, consider the operatorê iϕ S which fulfills Ŝ z ,ê iϕ S =ê iϕ S exactly [98], given bŷ where the omitted terms are formally suppressed for , the normalization given in Eq. (1) captures the correct behavior for large S.

II. EQUIVALENT SPIN MODEL AND HILBERT SPACE DIMENSION
We derive the effective spin model starting from Eq. (1) in one dimension, i.e.
Here and in the following, we we take N to be even and assume periodic boundary conditions. The Hamiltonian is gauge invariant, which can be expressed as the vanishing commutator Ĥ ,Ĝ j = 0 with the Gauss' law opera-torsĜ Additionally, the total fermion number is conserved, i.e. Ĥ ,N f = 0 with the fermion number operator As a first step, we deal with the fermionic degrees of freedom and perform a Jordan-Wigner transformation, We thus obtain the equivalent Hamiltonian whereα j =N f for j = N and zero otherwise. The full Hilbert space H = H S ⊗ H σ has dimension dim H = (2S(S + 1)) 2N , but the physical Hilbert space H phys = H| Gn=0 = {|phys ∈ H | G n |phys = 0} is much smaller. Within H phys , the configurations ofσ are uniquely determined by those ofŜ and thus we can express all operators in H phys in terms ofŜ. Let P be the operator that projects all possible spin configurations in H S to the allowed ones in H phys . EmployingŜ z j −Ŝ z j−1 = 1 2 σ z j + (−1) j on the physical subspace and shifting indices, we find where α j = N/2 for j = N and zero otherwise. In the main text, we have dropped the phase α , as it becomes irrelevant in the thermodynamic limit. We emphasize the importance of the projector P, which effectively leads to an interaction among neighboring spins and allows us to rewrite (1/2)P σ + jŜ + jσ − j+1 + h.c. P = PŜ x j P (see also [72] for the special case of S = 1/2).
To find the dimension of H phys , we need to count all possible configurations allowed by Gauss' law. In the basis ofŜ z , this is conveniently expressed by a square matrix C S of size 2S + 1, Here, rows and columns correspond to all possible spin configurations of two neighbors and 1 or 0 indicate an allowed or forbidden pair, respectively. Gauss' law alternates between even and odd sites as expressed by the matrices C S and C T S . The total Hilbert space dimension can now be obtained by summing over all allowed configurations, which is (for periodic boundary conditions) equivalent to the following trace The eigenvalues x (S) j of the matrix (S.20) In our main text and our exact diagonalization numerics, we make extensive use of the smallness of this dimension. The scaling of the dimension with S and N is illustrated in figure 5.

III. ANALYTICAL EXPRESSIONS FOR THE VECTOR AND SCALAR MASSES
The exact masses in the strong coupling limit (m/e = µ/g = 0) are known to be For the readers' convenience, we also quote the perturbative predictions which we used for the plots in the main text. We refer to [22] and references therein for details. Explicitly, the strong coupling expansions are given by

IV. EXTRAPOLATIONS
As described in the main text, the ground state energy, vector, and scalar mass of the continuum Schwinger model are obtained by a series of extrapolations. From the ED, we numerically obtain the lowest three (j = 0, 1, 2) eigenvalues E j = E j (S, N, g, m/e) in the zeromomentum sector of the Hamiltonian, grouped by several values of spin length S, lattice size N , coupling g = ae, and mass m/e. In a first step, we rescale all quantities according to their expected behavior in the continuum limit (see e.g. [75]), (S. 26) and the extrapolated function ω j (N, g, m/e) does not have the S-dependence any more. The 1/S dependence in this extrapolation is motivated as the simplest polynomial dependence on S, and also justified a-posteriori from the data. An exponential dependence would need a scale, and we have no reason to postulate the emergence of an additional lengthscale. Similarly, the large-N limit is obtained by the fits ω 0 (N, g, m/e) = ω 0 (g, m/e) + β 0 (g, m/e) N + γ 0 (g, m/e) N 2 , (S.27) ω 1,2 (N, g, m/e) = ω 1,2 (g, m/e) + γ 1,2 (g, m/e) N 2 , (S. 28) where we take into account that the leading finite size corrections of the excited states arise at second order in N [75]). Finally, we extrapolate to the continuum limit via ω 0 (g, m/e) = ω 0 (m/e) + δ 0 (m/e)g + 0 (m/e)g 2 , (S.29) ω 1,2 (g, m/e) = ω 1,2 (m/e) + δ 1,2 (m/e)g , (S. 30) where including the second order for the ground state significantly improves the accuracy of our results, indicating that our lattice data is rather far away from the continuum limit. Our final results are the ground state energy density ω 0 (m/e) and the first two energies ω 1,2 (m/e) above the ground state. The extrapolations of the iMPS data proceed analogously, omitting the fit with respect to N . In general, the data ω j (S, N, g, m/e) diverges as g → 0 for fixed S, N , and m/e in such a way that for every g there are minimal values of S and N required for a reliable extrapolation. In order to avoid additional assumptions about the functional form of ω j (S, N, g, m/e), we have implemented a simple convergence check to discard data with too small values of S and N : To obtain the results presented in the main text, we discarded all bare data corresponding to the same value of g (and fixed j) if the extrapolated values ω j (N, g, m/e) differ by more than 10% from the value corresponding to the largest available spin length S max , i.e. when |ω j (N, g, m/e) − ω j (S max , N, g, m/e)| > 0.1 × |ω j (S max , N, g, m/e)|. For the ED data, we checked for convergence of the N extrapolation analogously. After supplementing the fits by this procedure, and given the relatively small values of S (and N for the ED), we find that the uncertainty of all quantities is dominated by the final extrapolation to the continuum limit. The error bars shown in the plots of the main text indicate systematic errors that arise from the corresponding choice of fit range. We have estimated these errors following the procedure described in the appendix of [75]. Our final results are summarized in figure 6, where we also state the normalized square error  and N data is the number of data points y k used for the particular fit. We emphasize that the χ 2 -value has no statistical interpretation as all our calculations are purely deterministic.

V. PATH INTEGRAL REPRESENTATION
In this section, we outline a derivation of the path integral for the pure gauge U (1) QLM Hamiltonian using coherent states for the quantum spin operators representing the gauge fields, and indicate how the large representations give rise to the Wilson-Kogut-Susskind limit of the gauge theory. Given a Hamiltonian H in d−spatial dimensions which satisfies a local constraint, the corresponding path integral is given by Z = Tr e −βH P = DΩ e −S [Ω] . (S.32) The extent in Euclidean time, β equals inverse temperature T, and P is the projection operator projecting the configurations in a chosen computational basis to the ones allowed by the local constraint. The path integral is constructed by splitting the total Euclidean time β into N t Trotter steps of extent , such that β = N t . At each Trotter step, we construct the transfer matrix, Ω n−1 | T | Ω n = Ω n−1 | exp(− H) | Ω n , where |Ω n denotes the computational basis (which we introduce next) at time-slices n-1 and n. A periodic boundary condition in the imaginary time direction for the gauge field is used. The final expression is an integral over all possible field configurations Ω in (d + 1)−dimensions in the limit → 0, N t → ∞ with β held fixed. Further, for β → ∞ ground state results can be obtained.
The pure gauge Hamiltonian is given by H g = g 2 2 N j=1 (Ŝ z j ) 2 , and the Gauss' law without the matter fields is simplyĜ j =Ŝ z j −Ŝ z j−1 , and in the follow-ing we will denote the projection operator which selects configurations according to this (local) constraint as P G = j PĜ j . We emphasize that the fermions are not included in this Hamiltonian, and consequently the Gauss' law differs from the one in the main text. To derive the path integral, we use the coherent state basis (also known as Bloch states in the literature) as the computational basis, which is denoted as (j denotes a spatial site, and n the timeslice) |Ω n = j |θ j,n , φ j,n = j (cos θ j,n 2 ) 2S exp tan θ j,n 2 e iφj,n S − j,n |0 , where 0 ≤ θ j,n < π and 0 ≤ φ j,n < 2π. In what follows, we further need the resolution of identity in the coherent state basis at each timeslice n: sin(θ j,n )dθ j,n 2π 0 dφ j,n |Ω n Ω n | .
(S.34) Using the above expressions, in the coherent state basis, we write the partition function Z as The overlap Ω n−1 | Ω n constitutes an important piece in the full calculation and gives rise to the Berry phase. After some algebra, it can be expressed as (S.36) where we have assumed that is small such that the fields on adjacent time-slices are sufficiently smooth. The matrix elements of the following operators in the coherent state basis are also needed: (S.37) Combining the above formulae, the transfer matrix be-tween two time-slices is thus , (S. 38) where we have expanded the exponential, evaluated the operator, and re-exponentiated it, valid up to corrections at O( 2 ). The last timeslice contains the projection operator, and can be written in the |Ω -basis as: The ϕ j are the Polyakov loop phases which appear as the Lagrange multiplier enforcing Gauss' law on the last time-slice. To incorporate the Polyakov loop variables with the action, we consider the matrix elements of the Hamiltonian and the Gauss' law together on the last two timeslices: We collect the terms in φ j,N and perform the integral j π −π dφ j,N e iS(− cos θj,N+cos θj,0)φj,N which enforces the condition θ j,0 = θ j,N for all the spatial lattice sites (recall that 0 ≤ θ j,n < π). We can use this to simplify the matrix element × e [iS cos θj,0(φj,N−1−φj,0+ϕj −ϕj+1)] .
In the limit S → ∞, the effective action can be obtained by minimizing the quadratic term in cos θ j,n . This gives θ j,n ≈ π 2 ± α j,n , thus cos θ j,n ≈ cos π 2 − α j,n = sin α j,n ≈ α j,n , and sin θ j,n ≈ 1. The integral then becomes Gaussian, which can be performed explicitly. Further, we note that the term φ j,N−1 − φ j,0 = 2πn j + ∆ t φ j , where ∆ t φ j is the difference in the φ variable across the last timeslice. Note that the the same considerations as discussed here also hold for the other timeslices, as explicitly shown in Eq. (38). In fact, it is even simpler there, since the Gauss' Law does not need to addressed, and again the quadratic term can be extremized to obtain a Gaussian integral. The overlap of the coherent states, contributing to the Berry phase will be addressed separately.
As the analysis reveals, in the S → ∞ limit, the unit vector is forced to point at the equator, with minor fluctuations (denoted by α j,n ). At every time slice overlap, the φ j,n proceeds on the unit sphere, and thus the term φ j,N−1 − φ j,0 contains the total number of windings around the equator. In the large-S limit, we get The term in the exponential is the well-known Villain form [99] of the plaquette action in the Euclidean theory. Identifying ϕ j − ϕ j+1 ∼ −∂ j A 0 (0) and ∆ t φ j,0 ∼ ∂ t A 1 (0), the term in the exponential is the Maxwell action, and can be interpreted as the path integral in the axial gauge. In the axial gauge, all the time-like links are set to unity, except for the ones which cross the temporal boundary. In order to bring the expression to a form which is explicitly space-time symmetric, we can perform gauge transformations to spread the A 0 field away from the boundary time-slice and into the whole space-time lattice. The Boltzmann weight of the effective action contributing to Z is thus e −S eff = e where we have traded n and j for t and x respectively, and A 0,1 is the vector potential for the gauge field in 2-dimensions. We have thus demonstrated that how in the limit of large-S the localization of the spin operator along the equator is responsible for the emergence of the Wilson lattice gauge theory in Euclidean space. Next, we connect the contribution from the Berry phase term of the quantum spin operators to the presence of topological terms in the Wilson action. We note, following Coleman [63], that the Hamiltonian of the continuum massless Schwinger model in the axial gauge (A 0 = 0) can be written as where ψ,ψ are the fermion fields and γ 1,2 are the gamma matrices. In one dimension, due to the Gauss' law, the E is not an independent operator, but can be written in terms of the charge density up to an integration constant. It was shown by Coleman [63] that the integration constant acts as a background field, has the characteristics of a periodic function and is naturally related to the topological angle. Coleman further showed that for certain special values of this topological angle, Θ = π, the physics of the model is very different from the expected confining behavior.
In the context of QLMs, it is quite remarkable that the use of full or half-integer spin representations can induce qualitatively different behavior, at least in certain parameter regimes [10]. In particular, using the toolbox of coherent states, we can demonstrate that for large-S, half-integer spin representations can give rise to the topological angle Θ = π.
Note that this topological term is generated from the overlap of the coherent states Ω 0 |Ω 1 Ω 1 |Ω 2 · · · Ω N −2 |Ω N −1 = exp iS j ω j [Ω j ] (see Eq. (38)), which is the Berry phase in the problem. ω j [Ω j ] is the solid angle corresponding to the worldline (closed curve due to periodic boundary condition) of a spin at site-j as it evolves in imaginary time. While the Maxwell term, S eff is independent of S in the S → ∞ limit, this term is sensitive to the integer vs. half-integer nature of the spin S. To see this, note that Berry phase term can be written as [93] (we drop the spatial index without loss of generality) Here, θ is treated as a function of φ, as the spin S traces out a path S(τ ) on the surface of the 2-sphere, S 2 . Due to the periodic boundary conditions in time, the angle φ 0 must come back to itself. As we argued before, in the large-S limit, the unit vector is restricted to the equator with fluctuations, which makes the phase space effectively a circle. Due to the topological nature of this term, it only matters how many windings the vector makes around the equator. Thus in the leading order, the topological term contributes as 2πS. For half-integer S, this gives an effective phase π, which contributes the same way as the background field Θ = π in the Wilson lattice gauge theory. On the other hand, for integer values of S, the phase of 2π is irrelevant, and only the plaquette action remains. Thus, we have shown that in the large-S limit, the half integer values of S produce the same effective action as that of the Wilson theory with a background field of π, while the integer values produce an overall phase of multiples of 2π, and no extra effects. This derivation can be straightforwardly repeated with the inclusion of fermions.