Quantum simulation of the universal features of the Polyakov loop

Lattice gauge theories are fundamental to our understanding of high-energy physics. Nevertheless, the search for suitable platforms for their quantum simulation has proven difficult. We show that the Abelian Higgs model in 1+1 dimensions is a prime candidate for an experimental quantum simulation of a lattice gauge theory. To this end, we use a discrete tensor reformulation to smoothly connect the space-time isotropic version used in most numerical lattice simulations to the continuous-time limit corresponding to the Hamiltonian formulation. The eigenstates of the Hamiltonian are neutral for periodic boundary conditions, but we probe the nonzero charge sectors by either introducing a Polyakov loop or an external electric field. In both cases we obtain universal functions relating the mass gap, the gauge coupling, and the spatial size which are invariant under the deformation of the temporal lattice spacing. We propose to use a physical multi-leg ladder of atoms trapped in optical lattices and interacting with Rydberg-dressed interactions to quantum simulate the model and check the universal features. Our results provide a path to the analog quantum simulation of lattice gauge theories with atoms in optical lattices.

Lattice gauge theories (LGT) are fundamental to our understanding of strongly interacting particles in high energy physics. Translating the success of quantum simulations with cold atoms in optical lattices [1] of systems relevant to condensed matter such as the Bose-Hubbard model to the quantum simulation of LGT would open the door to real-time and finite-density calculations which are beyond the realm of classical computations. An important first step is to achieve this goal for models in one space and one time (1+1) dimensions. While the dynamics of the Schwinger model, quantum electrodynamics in 1+1 dimensions, has been explored using a few qubit digital quantum simulation in a system of trapped ions [2] or classical-quantum algorithms on IBM quantum computers [3], the analog quantum simulation of gauge theories with cold atoms requires complex experimental settings. Existing efforts involve mixtures of bosonic and fermionic atoms [4,5] or dipolar interactions of cold molecules [6] and are still in progress.
In this Letter, we propose to quantum simulate the Abelian Higgs model in 1+1 dimensions, the Schwinger model with the electron replaced by a complex scalar field. The validity of the proposal can be checked by measuring the Polyakov loop, an observable that plays an important role in finite-temperature studies [7] and for which we present remarkable finite size scaling (FSS) properties. We invoke a single atomic species in an optical lattice on a multi-leg ladder and recently explored Rydberg-dressed interactions in this platform [8]. The ladder is a physical lattice with a long side ("legs") representing one-dimensional space and a short side ("rungs"), with a shorter lattice spacing, representing rotor angular momentum. We aim at maximal simplicity both on the theoretical and experimental side. In contrast to other approaches [4][5][6][9][10][11][12][13][14][15][16][17][18], we use a manifestly gaugeinvariant formulation [19] where there is no need to enforce Gauss's law. In addition, we consider the limit [19] where the scalar self-coupling becomes large and the Higgs mode decouples from the low energy theory. We are then left with a gauged O(2) spin model with compact field integration. Fourier analysis provides a discrete reformulation, suitable for quantum simulations, in agreement with Pontryagin duality [20].
The remarkable FSS properties are illustrated by the collapse of 24 datasets in Fig. 1. There are four differ-ent spatial sizes N s represented in the figure (4, 8, 16, and 32) and it is possible to probe the critical behavior with systems of modest spatial sizes. We first give some brief explanations about these results and provide the details later in the text. The FSS is related to the energy gap ∆E created by inserting a Polyakov loop (a Wilson loop wrapping around the periodic Euclidean time) or by applying an external electric field. When the gauge coupling g approaches zero, we have an O(2) model and when the hopping parameters exceed their critical value at the Berezinsky-Kosterlitz-Thouless (BKT) transition, we have infinite correlation length at infinite volume and ∆E ∝ 1/N s at finite N s . Fig. 1 indicates that when we turn on g, N s ∆E is a linear function of (gN s ) 2 at small argument and then a linear function of gN s at larger argument. The two parts of Fig. 1 each contain 24 data sets with 12 from a discrete Euclidean time Lagrangian, while the other 12 are from the continuous-time Hamiltonian limit. It is remarkable that the two calculations provide the same universal functions. Ultimately, it is this equivalence between the two formulations which enables the transfer of results from the experimentally accessible Hamiltonian dynamics to the Lagrangian formulation.
The considered Abelian Higgs model is described by the lattice path integral Z = Dφ † DφDU e −S with action The complex (charged) scalar field is φ x = e iθx on spacetime sites x and the Abelian gauge fields U x,µ = e iAµ(x) on the links from x to x +μ. The electromagnetic tensor F µν = ∂ µ A ν − ∂ ν A µ appears when taking products of gauge fields around an elementary square (plaquette) in the µν plane. The notation for this product is U x,µν = e i(Aµ(x)+Aν (x+μ)−Aµ(x+ν)−Aν (x)) and the gauge coupling enters through β pl = 1/g 2 . The parameter κ is the hopping coefficient. The Fourier expansions of the Boltzmann weights lead to expressions of the partition function in terms of discrete sums of products of modified Bessel functions with integer orders on each plaquette and each link of the square lattice. The integer plaquette quantum numbers completely determine the integer link quantum numbers which are the difference of the neighboring plaquette quantum numbers which can be interpreted as dual variables. Explicit formulas and sign conventions are given in Ref. [19], where we also show that the discrete tensor renormalization group (TRG) [21,22] approach and the standard Monte Carlo approach give consistent numerical answers. The link quantum numbers can be interpreted as matter charges and their sum on the time links between two successive time slices stay constant as The Q = 0 sectors can be probed by inserting a product P , called the Polyakov loop, of gauge links in the Euclidean time directionτ with periodic boundary conditions: inside the path integral. This is illustrated in Fig. 2. Conventional Monte Carlo simulations and TRG calculations with a typical bond dimension of D bond ≈ 40 provide consistent evidence [23] for the behavior P can be interpreted in terms of the free energy induced by the inclusion of a static charge. A similar effect can be induced by using asymmetric boundary conditions such as 0 on one side and 1 on the other side for the plaquette quantum numbers. This situation can be interpreted as the introduction of an external electric field and is denoted as 01BC in Fig. 1 and hereafter. A similar situation can be accomplished by using 00BC and moving the Polyakov loop to the boundary as shown in Fig. 2b where an addtional column of zeros is implicit. In order to connect the LGT calculations in the Lagrangian formulation with quantum simulations using the Hamiltonian formulation, we need to take the time continuum limit. This is done [19] by taking κ τ , β pl → ∞ while simultaneously taking κ s , a → 0 (a is the temporal lattice spacing) such that the combinations are kept constant. Note that X here is related toX in Ref. [19] by X = √ 2X. In this limit, using the properties of the Bessel functions, the Hamiltonian can be identified as The sum i , taking 00BC into account, includes (L z 1 ) 2 + (L z Ns ) 2 . The operator U x = (U + +U − )/2 with the special type of ladder operators U + , U − defined by U ± |m = |m ± 1 , where L z |m = m |m with plaquette quantum numbers m = 0, ±1, ±2, . . . If we truncate at |m| max = s, we call it a spin-s truncation even though it is different from the rotation group representation used in Ref. [24].
In the continuous-time limit, the introduction of the Polyakov loop amounts to changing the Hamiltonian intõ where we have assumed that the Polyakov loop is put on the center of the lattice. The 01BC choice provides another way to probe the charge-1 sector. This simply changes (L z Ns ) 2 in the second summation of Eq.
The numerical continuous-time results were obtained using the density matrix renormalization group (DMRG) [25,26] to calculate the ground state energies for both cases. The finite DMRG algorithm with matrix product state (MPS) [27] optimization was performed using the ITensor C++ library [28]. Y = 1 units were used in all DMRG calculations.
We are now in position to provide more details about Fig. 1. As explained in more detail in [23], arguments regarding the behavior at small and large gN s led us to conjecture that N s ∆E is solely a function of the product (gN s ) 2 . Fig. 1 supports this idea and shows a good data collapse across multiple system sizes for both the discrete and continuous-time limits. Note that for the discrete-time (Lagrangian) calculations at various κ, ∆E was rescaled by 2κ while for the continuous time (Hamiltonian) calculations, no such rescaling was introduced for the different values of X. We emphasize that the collapse is by no means automatic. It breaks down for κ not large enough, if we increase g to large values while keeping N s constant (there are hints of this in Fig. 1 for the small N s data), or if the truncation value m max is too small.
Notice that in all cases, N s ∆E 0.5 when g 2 0, which corresponds to the gapless BKT phase in the O(2) limit. We also notice that the energy gap at finite g 2 for 01BC is bigger than the gap for PL-00BC. Because The upper part shows the possible mz-projections. Below, we show the corresponding realization in a ladder within an optical lattice. The atoms (green disks) are allowed to hop within a rung with a strength J, while no hopping is allowed along the legs. The lattice constants along rungs and legs are ar and a l respectively. Coupling between atoms in different rungs is implemented via an isotropic Rydberg-dressed interaction V with a cutoff distance Rc (marked by blue shading).
01BC breaks the inversion symmetry of the system, creating a charge on the left costs more energy than that on the right. We can understand this by doing the transfor- Then the Hamiltonian with 01BC is related to the Hamiltonian with the Polyakov loop by H 01 →H + U Ns i=Ns/2+1 L z i + N s U/4, which is simply adding a linear potential on the right half of the system. But the data collapse is still present, indicating the same universality class. As discussed in Ref. [23], a spin-6 truncation is good enough to probe the finite size scaling of the energy gap in the O(2) limit for 2 ≤ X ≤ 5 up to N s = 32. As a finite gauge coupling g will suppress the contribution of high plaquette quantum numbers, a spin-6 truncation should work well for any finite g and is used for all DMRG calculations in this work.
An important feature of the Hamiltonians considered above is that L z i has positive and negative eigenvalues and cannot be realized as the number operator of a Bose-Hubbard model unless a large chemical potential is introduced [24,29,30]. For this reason, two species of atoms were introduced in Ref. [24]. For a similar reason, a 2-leg ladder with 2s atoms per rung for a spin-s truncation was suggested in Ref. [19], however the hopping along a rung can only emulate the L x operator in the rotation group representation instead of the U x operator in Eq. (5) and (6). Here, we propose a simpler experimental realization to overcome this difficulty, namely an asymmetric ladder of N s rungs of length 2s + 1 each, with lattice constant a l and a r along legs and rungs respectively, see Fig. 3. The tunnel coupling along the legs is vanishing while it has a strength J along the rungs. The number of atoms per rung is held fixed at unity, such that the L z -projection of the spin is encoded in the position m of the atom within a given rung and can be read out with near-unity fidelity in a quantum gas microscope [1]. The initialization of the system can be achieved in such a setup by preparing an atomic Mott insulator and employing site-resolved opti- cal potentials [31].
The Hamiltonian of such a ladder system can be written aŝ Here, we have additionally introduced an interaction V m,m ,i,i between two particles at positions (m, i) and (m , i ) as well as an on-site potential m,i . The term X i U x i in Eq. (5) and (6) directly maps to the tunneling term in Hamiltonian (7)  While the tailored on-site potentials m,i can be gen-erated using optical potentials controlled at the singlesite level [32], realizing the quadratic distance dependence of the interaction between two particles is challenging in cold atomic gases. However, they can still be realized approximately using off-resonant optical coupling of the atoms to Rydberg states. The resulting isotropic Rydberg-dressed interactions [33,34] in cold atoms have recently begun to be explored in a manybody setting [8] and exhibit a characteristic distance dependence V (R) = U 0 /(1 + (R/R c ) 6 ) for two atoms separated by a distance R. The saturation value U 0 can be tuned to be positive or negative, and the interaction range R c is set by the interactions of the coupled Rydberg states and typically reaches up to several sites in an optical lattice [35].
The key idea in implementing quadratic interactions in the ladder model consists in utilizing an asymmetric ladder with different lattice constants along legs and rungs respectively. In the limit of large a l /a r , the interaction potential along the rung approximately acquires the desired quadratic distance dependence for neighboring rungs with |V 0 | = |U 0 |/(1 + (a l /R c ) 6 ) and Y = 6|U 0 |(a l /R c ) 6 (a r /a l ) 2 /(1+(a l /R c ) 6 ) 2 . At the same time, interactions between next-nearest-neighbor rungs can be minimized, see Fig. 4, making them irrelevant for the predicted collapse shown in Fig. 1. This and other imperfections as well as concrete experimental numbers are further discussed in [35].
A strength of the presented ladder implementation is the simple realization of models with different spin. A natural first step would be to check the experimental feasibility of the proposal with just two legs, i.e. s = 1/2 in Eq. (7). The emerging spin model corresponds to the well studied spin-1/2 quantum Ising chain in a transverse field with the Hamiltonian The transverse field is realized by the tunneling of the atoms and has a strength h x = J/2.
is required to realize the other two terms. Expressing all energies in units of the transverse field (h x = 1), this model has a second order phase transition at λ = 1 with known exponents [36]. As quantum simulations are still made on relatively small lattices, it is convenient to study the finite size scaling dictated by the Renormalization Group (RG) analysis of the secondorder phase transition. The zero temperature magnetic susceptibility reads where ... are short notations for Ω| ... |Ω with |Ω the lowest energy state ofĤ. The data collapse obtained with the standard RG rescalings is illustrated in Fig. 5.
The quantum Ising model has for example been quantum simulated in systems of ultracold ions [37] and with atoms in tilted optical lattices [38]. New generations of D-wave machines have more versatile time-dependent capabilities. It seems possible to maintain a transverse field [39] but there are temperature effects that need to be better understood. Multi-mode cavity photon-mediated interactions [40] can also be used to simulate the quantum Ising model. The possibility of extending these setups or related ones to reproduce a multi-leg ladder is being investigated.
In conclusion, we have presented an experimental platform for the quantum simulation of the Abelian Higgs model in 1+1 dimension and outlined a strategy for an initial benchmark of the quantum simulator. An interesting perspective is the experimental simulation of out-ofequilibrium dynamics following a quantum quench, which promises insight into dynamics described by the LGT when inaccessible with classical computing.

QUADRATIC INTERACTIONS IN A LADDER
For the quantum simulation of the Abelian Higgs model in spin-s truncation, we aim to implement a ladder system with quadratic interactions between particles in neighboring rungs of the ladder. The starting point is an isotropic Rydberg-dressed potential V (R) = U 0 /(1+(R/R c ) 6 ) with a cutoff distance R c = (C 6 /2∆) 1/6 and a saturation value U 0 = Ω 4 /8∆ 3 given by the laser coupling Ω, detuning ∆ and van-der Waals coefficient C 6 for the interaction between the coupled Rydberg states [33,34]. Quadratic interactions can be realized in a ladder with different lattice constant a l along the legs and a r along the rungs. To see this, we express the distance between two particles in terms of the rung and leg lattice constants and indices as where we have abbreviated the separations along the legs and rungs with ∆i ≡ |i−i | and ∆m ≡ |m−m |. Inserting this in the interaction potential V (R), we obtain for the nearest-neighbor rung with ∆i = 1 6U 0 (a l /R c ) 6 (a r /a l ) 2 (1 + (a l /R c ) 6 ) 2 ∆m 2 For the expansion, we have assumed (∆m a r )/a l 1 This allows us to express the potential given in the main text in terms of the experimentally relevant parameters U 0 , a l , a r and R c by identifying For simplicity, we have hereby assumed attractive interactions and written U 0 = −|U 0 |. As the interactions between next-nearest-neighbor rungs are assumed to be vanishing in the desired spin model in Eq. (5), it is important to check that they are small enough in this approximation. The next-nearestneighbor-rung interactions (NNNRI) are obtained by just setting a l → 2a l in the above parameters, yielding In order to continue the discussion, it is now helpful to gain some intuition for where the quadratic dependence is a good approximation to the Rydberg-dressed potential. One exemplary Rydberg-dressed potential is depicted in the main text in Fig. 4. As one can see in this example, the points for the different ∆i are grouped. This is generally necessary to achieve a large ratio between ∆i = 1 and ∆i = 2, i.e. a suppression of NNNRI. Furthermore, the ∆i = 1 points should be located close to the most linear part of the potential R ≈ R c . This can be understood from the geometric argument that in the limit a r /a l → 0, where the quadratic approximation of R in ∆m works best, we sample a linear potential with approximately quadratically spaced points. At the same time, however, pushing a r /a l towards zero squeezes the points closer together, effectively reducing the interaction Y . This is also directly obvious from the quadratic dependence of Y on a r /a l shown in Eq. (13). To summarize, one has to compromise between maximizing Y , minimizing Y (2) and optimal quadratic dependence within a rung. Qualitatively, Fig. 4 and our previous discussion indicate to work close to the regime R c ≈ a l . Setting for simplicity R c = a l , we can estimate the ratio between nearest-neighbor and next-nearest-neighbor rung interactions as (approximating 2 6 + 1 ≈ 2 6 ) This shows that NNNRI are suppressed by approximately |Y |/|Y (2) | = 2 6 = 64 in this regime and hence only reach up to 1.5% of the nearest-neighbor-rung interactions (NNRI). Furthermore, it should be noted that already small changes in a l /R c influence the exact numbers due to the strong scaling with (a l /R c ) 6 , which allows for considerable flexibility in fine-tuning the implementation. For example, increasing a l /R c slightly above unity suppresses the NNNRI to values below 1%, such that they can easily be tuned to be irrelevant for the collapse displayed in Fig. 1. Fig. 6 shows N s ∆E obtained from the Hamiltonian of the ladder system with NNNRI 1.5% of NNRI by DMRG. The data collapse is robust to NNNRI for the ladder systems which are mapped to the spin Hamiltonians with and without 01BC, while it is broken dramatically by such a small NNNRI in the ladder systems which are mapped to the spin Hamiltonians with and without a Polyakov loop. In the O(2) limit, NNNRI doesn't change the critical property for the former, while it would open a gap to the gapless phase for the latter, which has been confirmed numerically (not shown here). However, based on 1st order perturbation theory, we can correct the ground state energy by subtracting the expectation value of NNNRI, H N N N RI , from it. The corrected data sets collapse perfectly onto the same lines of collapse in Fig. 1. Experimentally, we can measure the densitydensity correlations between next-nearest-neighbor rungs n m,inm ,i+2 and extract H N N N RI to do the same correction.

Imperfections due to decoherence
Experimentally, the achievable Rydberg-dressed interaction strength U 0 is limited via the effective decay rate γ eff of the dressed state due to the coupling to a Rydberg state with a limited lifetime [8,33]. This leads to decoherence in the time evolution of the system not captured by purely Hamiltonian dynamics described by Eq. (7). Generally, both the Rydberg-dressed interaction strength as well as the effective decay rate increase with larger Rydberg-state admixture. One can quantify the coherence in terms of a "figure of merit" Q = |U 0 |/γ eff , which expresses the average number of coherent cycles per decay event.
For simulating the Abelian Higgs model in the spin representation, the additional constraint of realizing quadratic interactions effectively reduces the usable interactions. Considering Fig. 4, this is directly obvious as the strongest interactions (∆ i = 0) are not sampled for a single atom per rung. A useful figure of merit is given by comparing the quadratic part of the interaction Eq. (13) with the effective decoherence, Y /γ eff = 6 (a l /R c ) 6 (1 + (a l /R c ) 6 ) 2 (a r /a l ) 2 |U 0 |/γ eff = 6 (a l /R c ) 6 (1 + (a l /R c ) 6 ) 2 (a r /a l ) 2 Q.
This shows that realizing an optimal quadratic dependence, which requires ideally a r /a l → 0, as well as maximally suppressing the NNNRI, which requires R c /a l → 0, conflicts with the interaction-to-decoherence ratio. Below, we outline a parameter regime showing that the experimental implementation still seems feasible.

MEASUREMENT OF THE ENERGY GAP
In order to experimentally measure the energy gap displaying the universal scaling (see Fig. 1), it is necessary to first prepare the ground states of the respective models. Their mean energy E , with . . . denoting the quantum average, can then be reconstructed by two sets of measurements. In a first set of measurements, the atomic distributions are immediately frozen by switching off interactions and tunneling at the same time. The contributions of on-site potentials, Ns i=1 s m=−s m,i n m,i , and interactions, s m,m =−s V m,m ,i,i n m,inm ,i , to the mean energy can then be extracted from the measured atomic distribution by measuring the mean local density n m,i and density-density correlations n m,inm ,i respectively. Both are accessible in a quantum gas microscope with local detection [41]. In a second set of measurements, the contribution of the tunneling term to the mean energy has to be determined. This amounts to extracting first-order correlations of the form â † m,iâ m+1,i , which has recently been demonstrated for atoms in optical lattices employing Talbot interferometry [42]. The reconstruction of the energy relies on knowledge of the on-site potentials, interactions and tunnel coupling along the rungs, each of which can be calibrated in independent measurements.

EXPERIMENTAL NUMBERS
In the following, we give some concrete experimental figures to underline the feasibility of realizing the desired spin models. Hereby, we assume tunneling along the rungs with strength J/h = 100 Hz, which is e.g. readily achieved for rubidium in optical lattices [31].

Spin-s model with quadratic interactions
First, we focus on the feasibility of reaching the parameter regime for the ladder implementation of the spin-s truncation of the Abelian Higgs model to study the collapse displayed in Fig. 1. Aiming at J = Y (X = 1 in Fig. 1), the required interaction strength for R c = a l ≈ 7 a r is |U 0 |/ ≈ 2π × 3.3 kHz (using Eq. (20)). Assuming a coupling strength of Ω/ = 2π × 100 MHz to the Rydberg state, which is within reach in future, specialized experiments, the desired interaction |U 0 | and hence Y can be achieved with a detuning of ∆/ ≈ 2π ×1560 MHz. The interaction strength between two Rydberg atoms with 80P 3/2 , measured by the van der Waals dispersion coefficient, is C 6 ≈ 5500 GHz µm 6 . As a result, for the chosen detuning, the cutoff distance is R c = (C 6 /2∆) 1/6 ≈ 3.5 µm and hence an experimentally realistic lattice spacing of a r ≈ 500 nm would lead to R c /a r ≈ 7. For the lifetime τ = 250 µs of the admixed Rydberg state 80P 3/2 , the dressed state acquires an effective decay rate of γ eff ≈ 4 s −1 , resulting in J/γ eff ≈ 25. Tuning to larger values of J/Y by reducing Y to study the collapse shown in Fig. 1 improves this ratio by a factor J/Y for constant ∆ due to smaller Rydberg-state admixture. A direct increase in the figure of merit should be possible by using lighter elements like potassium or lithium, for which larger tunneling rates are feasible. This allows for working at a larger Rydberg-state admixture, which is favorable to increase the quality factor. Furthermore, improvements are possible by increasing Rabi couplings to Rydberg states, by implementing more advanced dressing schemes [43] or in future experiments in a cryogenic environment, where the Rydberg lifetime is expected to increase (to 1.2 ms for 80P ) due to suppression of black-body-radiation-induced decay. Combining these steps, we estimate a possible increase in the figure of merit by well above one order of magnitude.

Spin-1/2 Ising model
The spin-1/2 Ising model (Eq. (8)) discussed in the main text constitutes a first, easier step for implementing the described ladder systems. Contrary to the more complex spin-models, in this case no asymmetric ladder is required, i.e. a r = a l . Keeping the condition R c = a l , in this case λ ≈ U 0 /4 using the same definition as in Eq. (8). For a Rabi-coupling of Ω/ = 2π × 100 MHz and a detuning of ∆/ ≈ 2π × 1560 MHz, the accessible spin interaction strength is λ/ = 2π × 825 Hz. A tunneling strength of J/ = 2π × 100 Hz along the rungs translates to a transverse field of h x / = 2π ×50 Hz, such that λ/h x ≈ 16.5. For a coupling to the Rydberg state 80P 3/2 with an effective decay rate of γ eff ≈ 4 s −1 for the above admixture parameters, the figure of merit in this regime becomes h x /γ eff ≈ 12.5. Towards stronger transverse fields, the figure of merit increases proportional to h x /λ, such that it reaches up to 50 in the critical region