Z(N) lattice gauge theory in a ladder geometry

Under the perspective of realizing analog quantum simulations of lattice gauge theories, ladder geometries offer an intriguing playground, relevant for ultracold atom experiments. Here, we investigate Hamiltonian lattice gauge theories deﬁned in two-leg ladders. We consider a model that includes both gauge boson and Higgs matter degrees of freedom with local Z N gauge symmetries. We study its phase diagram based on both an effective low-energy ﬁeld theory and density matrix renormalization group simulations. For N (cid:2) 5, an extended gapless Coulomb phase emerges, which is separated by a Berezinskii-Kosterlitz-Thouless phase transition from the surrounding gapped phase. Besides the traditional conﬁned and Higgs regimes, we also observe a novel quadrupolar region, originated by the ladder geometry.

Under the perspective of realizing analog quantum simulations of lattice gauge theories, ladder geometries offer an intriguing playground, relevant for ultracold atom experiments.Here, we investigate Hamiltonian lattice gauge theories defined in two-leg ladders.We consider a model that includes both gauge boson and Higgs matter degrees of freedom with local ZN gauge symmetries.We study its phase diagram based on both an effective low-energy field theory and density matrix renormalization group simulations.For N ≥ 5, an extended gapless Coulomb phase emerges, which is separated by a Berezinskii-Kosterlitz-Thouless phase transition from the surrounding gapped phase.Besides the traditional confined and Higgs regimes, we also observe a novel quadrupolar region, originated by the ladder geometry.

I. INTRODUCTION
Gauge theories are both the backbone of the standard model of particle physics and the key to understand a wide variety of condensed matter systems [1].Their pervasive importance, however, is flanked by the extreme difficulty in obtaining exact solutions for such strongly correlated models: many non-perturbative phenomena of quantum chromodynamics and other gauge theories remain open challenges at the core of intense research efforts [2].These difficulties prompted a long-standing endeavor in the simulation of gauge theories, generally based on the framework of lattice gauge theories (LGT) [3] and Monte Carlo techniques, which achieved many accurate results, including, for example, the definition of both the hadrons and light mesons spectra [4].
In addition to these traditional simulations, novel strategies to analyze gauge theories are being explored in the last years, based on the knowledge acquired in the field of classical and quantum simulations of many-body quantum systems (see, for example, the reviews [5][6][7][8]).These efforts develop on several directions and include new effective theoretical approaches to LGTs, the experimental realizations of the building blocks for their quantum simulation (for instance in trapped ion [9] and ultracold atom systems [10][11][12][13][14]), and tensor network calculations based on the Kogut and Susskind Hamiltonian formulation of LGTs [15].
The general development of these novel approaches relies on implementing progressive steps of increasing complexity on several levels.On one side, Abelian LGTs have been the first basic platform to test these techniques, before considering non-Abelian models.On the other, one-dimensional (1D) quantum systems offer the easiest playground to test tensor-network simulations before addressing higher dimensions.
Concerning the simulation of Abelian gauge theories and quantum electrodynamics (QED), for most numerical and experimental quantum simulations, it is useful to restrict the number of degrees of freedom in the considered many-body systems.Two main possibilities have been considered: (i) to maintain a continuous U(1) gauge symmetry by truncating the maximal value of the electric field flux propagating in each link; this is consistent with a quantum link model approach [16][17][18] or with a truncation of the choices of the gauge boson states corresponding to different representations of the U(1) group [20][21][22]; (ii) to consider Z N gauge theories [23][24][25][26][27][28][29][30][31][32] which reproduce the U(1) physics in the large N limit and rapidly converge to the exact compact QED observables in one dimension [18] (see also [19]).
Concerning the dimensionality of the systems, very recently there have been first attempts to investigate two-dimensional models with Abelian gauge symmetries.Truncated U (1) models have been studied based on both theoretical [33,34] and finite-size tensor-network investigations [35]; the phase diagram of the pure Z 3 LGT, instead, has been studied with both infinite-size tensor networks [31] and finite tensor networks combined with variational Monte Carlo procedures [32].
In this work, we address Z N LGTs in the geometry of ladder systems.This geometry offers an interesting compromise between one and two dimensions: on one side, it is the simplest geometry featuring plaquette interactions, thus enabling a full investigation of the Kogut and Susskind Hamiltonian; on the other, its quasi 1D nature allows us to develop an effective quantum field theory based on bosonization [36], which guides us in the exploration of the phase diagram of the model.Furthermore, ladder geometries have been very recently adopted for small-size digital quantum simulations with superconducting qubits of a truncated non-Abelian SU(2) lattice gauge model [37].
Our aim is to investigate the Z N models independently on their continuum and U(1) limit.Models with discrete Abelian gauge symmetries, indeed, have recently been a focus of attention on their own; for instance, several Z N symmetric models have been discussed in the context of 2D topological order [38][39][40][41].Furthermore, the last generation of experimental platforms for quantum simulators based on Rydberg atoms displayed the emergence of phases with discrete Z N symmetries [42].
In the following, we will investigate the phase diagram and main features of the Z N Abelian gauge theories with  (3).Squares and circles refer to Higgs matter and gauge boson degrees of freedom lying, respectively, on the vertices and links of the ladder.The yellow circle depicts the electric field energy, the orange square the matter mass, the green link represents a rung tunneling term and the blue plaquettes displays the magnetic plaquette interaction.Finally, in red, we depicted one of the gauge constraints (dashed ellipse).The dashed purple rectangle depicts the first unit cell on the left (rough boundary); the boundary conditions on the right, instead, are of the smooth kind.
Higgs matter.In Sec.II we introduce the model in the ladder geometry.In Sec.III we investigate the pure lattice gauge theory limit and we show that the electric field term in this quasi-1D system always dominates over the magnetic energy, thus leading the model into a confined phase for every N .In Sec.IV we introduce the Higgs matter and we set up a low-energy field theoretical description of the model based on bosonization.Sec.V is devoted to the renormalization group (RG) analysis and density matrix renormalization group (DMRG) simulation of the model, focusing on its thermodynamic phases and their properties.We show that for N = 2, 3, 4 the phase diagram is, in general, trivial and displays only a single gapped phase interpolating between the confined and Higgs regimes.For N ≥ 5, instead, an extended gapless phase appears which we interpret as a Coulomb phase.In Sec.VI we discuss the possible extension of the model to a larger number of legs and in Sec.VII we present our conclusions.The Appendices are devoted to several details of the analysis of the model and its renormalization group study.

II. THE GAUGE THEORY IN THE LADDER GEOMETRY
The quantum simulation of ladder systems has been broadly investigated in connection with the introduction of artificial gauge fluxes; in ultracold atom setups, ladders can be realized either in real two-dimensional systems [43] or using inner degrees of freedom to implement a synthetic dimension [44,45].On the theoretical side, their study is closely related to the coupled wire construction of interacting topological phases of matter (see [46] and references therein).
In connection with lattice gauge theories, ladder models offer the simplest scenario to study the plaquette interactions which are at the basis for the appearance of deconfined and topological phases in two space dimensions.In the following, we consider an Abelian Z N theory, in which the connection and plaquette operators are unitary [47].Our model is composed by bosonic Z N gauge field degrees of freedom living on the edges of the ladder, which represent N possible values of the electric field, and "frozen" Higgs matter degrees of freedom lying on the ladder vertices, which represent N different charge states (see Fig. 1).Both the edges and vertices, therefore, are characterized by an N -dimensional local Hilbert space and we introduce two pairs of clock operators acting on the gauge and matter degrees of freedom respectively.The first pair is given by τ ∼ e i 2π N E and σ ∼ e iA , which affect the gauge degrees of freedom and are respectively related to the electric field operator E and the magnetic connection A. They obey the algebra of Z N clock operators: N τ σ .
(1) In the following, we will label with indices σ r,↑ and σ r,↓ the clock operators on the links in the upper and lower leg of the ladder at position r.σ r,0 refers instead to the clock operators along the rung r.The same indices apply to the τ operators.
The second pair of clock operators is given by ζ r,y and η r,y which, instead, act on the Higgs matter degrees of freedom lying in the upper (y = ↑) or lower (y = ↓) legs.ζ, in particular, represents the phase of a Higgs field, whose radial mode is frozen to unity (London limit).Therefore, ζ and ζ † respectively annihilate and create a Z N electric charge, whereas η = e i 2πq N is linked to the charge operator q = 0, . . ., N −1 defined modulo N .They obey the same onsite algebraic relations as the previous clock operators: N ηζ .
(2) Based on these definitions, the Kogut-Susskind Hamiltonian [15,48] in the ladder geometry is: The terms appearing in this Hamiltonian are represented in Fig. 1.The first corresponds to the plaquette interaction; it is a function of the σ connection operator which define the mass of the Z N magnetic fluxes in the r th plaquette of the ladder.The second term, with coupling constant g, provides a dynamics to the magnetic fluxes and defines the electric field energy density on each link of the ladder.The term in 1/λ represents the mass of the electric charges of the model, whereas the last line corresponds to their tunneling mediated by the gauge degrees of freedom.In this work, we consider λ and g as free parameters, which are not directly related to the continuum U(1) theory.In this respect, we adopted the notation in Ref. [1] and we label by g the coupling constant for the gauge boson interactions, rather than using the standard particle physics g 2 notation [15].The Hamiltonian (3) corresponds to rough boundary conditions on the left side and smooth boundary conditions on the right side (see Fig. 1).In this situation, the ladder geometry can be described in terms of L unit cells, each including 2 matter sites and 3 gauge-boson links.Such Hamiltonian can be supplemented by boundary operators: This additional term allows for single matter charges to enter or leave the system from the left boundary, thus breaking their global charge conservation if λ b = 0.The total Hamiltonian H + H left bound.is symmetric with respect to the following local gauge transformations for the bulk vertices of the ladder: with 1 < r < L − 1; depending on the boundary conditions, additional boundary gauge constraints may appear.For our choice of rough-smooth boundaries we have the boundary gauge symmetries: The physical Hilbert space (without static charges) is defined by the gauge contraint: G r,s |ψ phys = |ψ phys , for each r, s , which imposes a Z N Gauss law on each vertex.The gauge constraints can be used to rewrite the Hamiltonian in specific gauge choices.In the following we will adopt either the axial gauge, in which the local gauge transformations are used to set all the gauge degrees of freedom along the two legs to the trivial state σ r,↑/↓ |ψ axial = |ψ axial , or the unitary gauge in which all the matter sites are set to the trivial state ζ r,y |ψ uni = |ψ uni .The latter choice is adopted in our tensor network simulations.

III. CONFINEMENT OF THE PURE GAUGE THEORY
We begin our analysis from the pure gauge limit (λ → 0), in which the matter degrees of freedom are frozen in η r,y |ψ = |ψ .This limit is conveniently studied in the axial gauge in which the Hamiltonian is expressed as a function of the rung degrees of freedom only (we drop, for convenience the index 0): Based on the axial gauge choice, the plaquette term is mapped on a ferromagnetic coupling between the gauge bosons on neighboring rungs, whereas the electric field energy of the rung is left invariant.These two terms are therefore mapped into a standard one-dimensional Z N quantum clock model (see, for example, [49,50]).The electric field energy of the gauge bosons in the links, instead, corresponds to the last non-local term in H gauge , and it distinguishes this LGT from the model studied in [51] in the context of topologically ordered systems.This term can be easily derived by considering the products of all the gauge constraints L j=r G j,s for the s leg.These string operators are equivalent to the identity on the physical states and, for the pure LGT, relate the electric field operator τ r,s on both legs with the string operators appearing in Eq. (8).
The pure LGT Hamiltonian with the boundary term (4) enjoys a global Z N symmetry given by the t'Hooft string G = r τ r,0 .This corresponds to the Z N symmetry underlying the clock model defined by the first two terms of Eq. ( 8) only.For our choice of boundary conditions, this symmetry directly appears as the first of the non-local terms (r = 1) in the second line of (8), proportional to the coupling g.
The clock model (local terms of H gauge ) is characterized by two gapped phases for N = 2, 3, 4: an ordered phase, in which the Z N G symmetry is broken for small g, and a disordered phase for large g.For N ≥ 5 a critical phase appears at intermediate values of g [50,[52][53][54][55][56].This picture is modified by the addition of the non-local terms, which include all the disorder operators acting on each link.Physically these operators correspond to the introduction of quanta of magnetic flux ±2π/N on each plaquette and these disorder operators always favor the disordered phase.For each finite value of N , H gauge displays indeed only a single gapped (disordered) phase for each value of g > 0.
The only exception is given by the limit g → 0 in which only the plaquette terms survive.In this limit, the ground state is a state without any magnetic vortex and, for rough-smooth boundary conditions, it presents an N -fold degeneracy corresponding to the breaking of the t'Hooft G symmetry.This symmetry can be explicitly broken by different boundary conditions, including, for example, rough boundary conditions with an additional 3-site boundary plaquette terms.
The non-local Hamiltonian H gauge can be mapped through a bond-algebraic duality [50] to the local Hamiltonian of a quantum clock model with both transverse and longitudinal fields, generalizing the analogous Z 2 Ising model.The Ising model in transverse and longitudinal fields displays indeed only a single gapped phase (see, for example, [57,58]), and we show in Appendix A that the same is true for its Z N generalization.
The analysis of the properties of the gapped phase of the pure lattice gauge theory for g > 0 is conveniently performed by considering perturbation theory in the two limits g → ∞ and g → 0. In the first limit, the ground state is just the product state of rungs displaying zero electric field, thus τ r,0 |ψ = |ψ .The plaquette operators with small amplitude 1/g introduce pairs of gapped local excitations which do not qualitatively modify the paramagnetic ground state.For g → 0, instead, it is easy to see that the degeneracy of the ferromagnetic ground states is split by the term gG appearing in H gauge and corresponding to τ 1,s .Furthermore, we observe that the non-local symmetry G is mapped into an holographic symmetry in the dual model [59] and the resulting ground state is symmetric under G for any g = 0 (see Appendix A).
This gapped phase corresponds to a confined phase for static charges; whereas only the limit g → 0 results in a deconfined phase.This can be proved by the introduction of static charges in the system through a violation of the gauge constraint in arbitrary pairs of sites.In particular, we consider a system in which we introduce a pair of opposite static charges in the sites x and y of the lower leg.This is done by imposing that, in these sites G x,↓ |ψ = e i 2π N |ψ and G y,↓ |ψ = e −i 2π N |ψ (see Appendix A for more detail).We find that the interaction energy between the two charges grows linearly with the distance R = |y − x| for any g > 0. We estimated the string tension T as a function of g through perturbation theory close to the limits g → 0 and g → ∞ (see Fig. 2 and Appendix A), and the results of our DMRG simulations show that the string tension interpolates between the two predicted behaviors, as shown in Fig. 2 for N = 5.
We conclude that the pure Z N lattice gauge theory in the ladder is always confined and, in this respect, it behaves as a 1D LGT.
It is interesting to inspect the behavior of the expectation values of the electric fields, E = −i(N/2π) log τ , in the presence of the two static charges.No phase transition characterize the pure lattice gauge theory at finite g; however, we can distinguish two different regimes.For large values of g, when the string tension is strong, the electric field propagates in a straight line from one static charge to the opposite when the two of them are on the same leg [Fig.3 (a)].For larger and larger g, indeed, the ground state progressively becomes a product state of the τ eigenstates in all the links and the energy cost of prolonging the electric field lines from one leg to the other is excessive.For small g, instead, the expectation value of τ decreases in modulo.The average electric field, however, propagates in both legs in the region between the two static charges [Fig.3 (b)], as dictated by the strong plaquette interactions.

IV. AN EFFECTIVE FIELD THEORY DESCRIPTION
A. The quantum clock model limit The next step in our analysis is to examine the role of the Higgs matter degrees of freedom, in the spirit of the seminal work by Fradkin and Shenker [48].To investigate the physics of the Hamiltonian (3), we begin our study from the limit g → 0. In this limit, the gauge bosons are frozen in a state without any magnetic excitation in the plaquettes of the ladder.Therefore, it is easy to rewrite  The system has rough and smooth boundaries at the left and the right sides respectively.For N = 3 and 4, the dashed lines indicate the proposed critical values of λ.For N = 5, the gray region represents the gapless phase determined by the fidelity susceptibility in a system of smooth-smooth boundaries (see Fig. 6).
the Hamiltonian in the axial gauge by imposing all the link states to be aligned in such a way that σ r,s |ψ = |ψ [1].Hence, the Hamiltonian H takes the form of a Z N quantum clock model on the ladder geometry.
The first line corresponds to a ferromagnetic interaction between any pair of neighboring clock operators representing the Higgs matter, both along the rungs and the legs of the ladder.The second line can be interpreted as the sum of the disorder operators in all the ladder vertices.
Quantum clock models of this kind are in general characterized by a symmetry broken ordered "ferromagnetic" phase for large λ and a disordered "paramagnetic" phase for small values of λ.We stress that, in the lattice gauge theory, the order parameter ζ is not a well-defined gauge invariant operator.The corresponding order parameter can be written in a gauge-invariant form only when considering at least an edge with "rough" boundary conditions by including the boundary interaction (4), which explicitly breaks the global conservation of the matter charge associated to the symmetry r,s η r,s .The boundary term (4) allows us to introduce a gauge-invariant order parameter O r,s = r j=1 σ † j,s ζ r,s which matches ζ r,s in the axial gauge.
In Fig. 4 we illustrate the value of this order parameter in ladders with rough-smooth boundary conditions for N = 3, 4 and 5.One can see a phase transition between the ordered and the disordered phases for N = 3 and 4. For N = 5, an intermediate gapless phase appears as indicated in the gray region in Fig. 4(c).
The appearance of a gapless phase in this model for N > 4 is reminiscent of the study of ferromagnetic onedimensional quantum clock models.In these 1D chains, it is well known that the quantum clock model displays an extended gapless phase for N ≥ 5 separating the gapped ordered and disordered phases (see, for example, [50,52,53]).The phase transitions between them are of the Berezinsky-Kosterlitz-Thouless (BKT) kind [55].Similar properties characterize the ladder clock model in Eq. ( 9), such that for N ≥ 5 a gapless phase appears in the system for intermediate values of λ.In the following, we will investigate the properties of this gapless system by describing the low-energy sector of the theory through an effective field theory inspired by bosonization, and we will numerically examine its main features through DMRG simulations.

B. Bosonization of the model
To build an effective low-energy description of the model we construct a representation of the clock operators based on vertex operators of a pair of dual bosonic massless fields, θ and ϕ.The following construction matches the dual sine-Gordon model description of 1D systems with Z N symmetry presented in Ref. [52] and it is inspired by standard bosonization techniques [36].A similar strategy has also been recently applied to the study of the Z 2 lattice gauge theory on the chain [60].
Our first step is to introduce the pairs of dual bosonic massless fields θ s (x) and ϕ s (x), with s = 0, ↑, ↓ that fulfill the following commutation relations: where Θ is the Heaviside step function with Θ(x ≥ 0) = 1 and Θ(x < 0) = 0. Based on this commutation relations, we build the following mapping between the clock operators of the ladder model expressed in the axial gauge and the vertex operators, in such a way that the algebraic properties of the clock operators (1) and (2) are satisfied.For the Higgs matter operators (s = ↑, ↓) the mapping reads: ja) , η j,s → e −iϕs(ja)+iϕs(ja+a) .
For the gauge bosons on the rungs we analogously impose: σ j,0 → e −iθ0(ja) , τ j,0 → e −iϕ0(ja)+iϕ0(ja+a) .( 12) In these relations we introduced the lattice spacing a, which is useful to define a proper ultraviolet cutoff of the theory.In particular, we consider the bosonic fields ϕ and θ to vary slowly in space with respect to the length scale set by a.It is easy to verify that the previous definitions fulfill (1) and (2) based on Eq. ( 10) (see Appendix C).The physical interpretation of the bosonic fields can be deduced as well from the previous equations.The fields ϕ ↑,↓ represent the electric field propagating along the legs, whereas θ 0 is associated to the magnetic field flux along the rungs.Some care is required in dealing with the boundary conditions: smooth boundary conditions on the ladder, for example, impose Dirichlet boundary constraint on ϕ ↑,↓ , since they require that no electric flux is allowed to enter the system from outside.Rough boundary conditions, instead, impose Dirichlet boundary constraints on θ ↑,↓ .Adopting the axial gauge, we obtain the following effective Hamiltonian on the continuum (see Appendix C for further details): In this Hamiltonian, the s = ↑, ↓ contributions of the first line account for the tunneling term along the legs of the ladder and the onsite term for the mass of the charges.The s = 0 contribution of the first line describes instead both the plaquette term and the electric field interaction along the rungs.These terms of the Hamiltonian in the axial gauge can indeed be mapped into a 3-component Luttinger liquid.The mapping between clock and vertex operators suggest that, in proximity to the gapless phase, K ↑ = K ↓ ≈ 1/λ and we expect K 0 to be proportional to g close to g = 1.The velocity is the same for all sectors: v = 4πa/N .The second line in the Hamiltonian (13) describes the rung tunneling term; the third line corresponds instead to the electric field interaction along the leg links, which has a non-local description in the axial gauge but recovers its locality in this description.
For later convenience we labeled their coupling constants as T and G such that: The final terms in the Hamiltonian ( 13) are aimed at restoring the Z N symmetry of the model and we will refer to them as "background interactions".The mapping (11) and ( 12) promote indeed the clock operators from discrete operators to continuous rotors (see, for example, [61]).The background interactions have the role of breaking the system symmetries from U (1) to Z N , consistently with the field theoretical description of clock models [52,56,62].The values of the constants P s and Q s can be roughly estimated by comparing the energy of the kinks of these sine-Gordon interactions with the energy of the domain walls in the corresponding operators (see Appendix C): Let us finally observe that the Hamiltonian ( 13) is invariant under the global Z N transformation θ ↑/↓ → θ ↑/↓ + 2π/N .The fields ϕ ↑,↓ , instead, do not enjoy such discrete global symmetry due to the leg electric field term.The Z N transformation ϕ ↑/↓ → ϕ ↑/↓ + 2π/N corresponds to the addition of a background electric field, thus to a change of the θ vacuum of the theory.This is analogous to similar features in truncated 1D models with U(1) gauge symmetry (see, for example [63][64][65][66]).

C. Properties of the clock model limit
We can obtain the main features of the system in the limit g → 0 by considering the Hamiltonian (13).In this limit θ 0 = 0 everywhere (σ r,0 = 1), consistently with the Hamiltonian (9), and the electric field term disappears.We are effectively left with a two-component Luttinger liquid perturbed by the interleg tunneling and the background interactions.As customary in these cases [36], it is convenient to separate symmetric "charge" ρ and antisymmetric "spin" σ combinations of the fields (the terms "charge" and "spin" are taken from the study of one dimensional two-component fermionic systems): The effective bosonized Hamiltonian reads: A first-order renormalization group analysis shows that, for N = 2, the system at g = 0 displays, as expected, two gapped phases separated by a critical point.For N > 2, the competition between the rung tunneling term and the background Q term is non-trivial and yields the possibility of having a phase in which only the spin sector is gapped.This can be understood by comparing the scaling dimensions of the interactions.The scaling dimension of the rung tunneling is D T = K σ /N , whereas the scaling dimensions of the background interactions are In particular, the bare value of the Luttinger parameters matches K σ = K ρ ≈ 1/λ and we will label it with K only.Based on the previous scaling dimension, for K ∈ 4/N, N/ √ 2 , the rung tunneling T term is the dominating interactions, with the P interaction being irrelevant and the Q interaction being suppressed by T .Therefore T gaps the spin sector of the system.
Moreover, a two-step renormalization group analysis (see Appendix E) shows that additional emergent secondorder terms gap also the charge sector for all values of λ in the cases N = 3, 4; only for N > 4 an extended gapless phase appears in the phase diagram, characterized by a gapless charge sector (see the next section and Appendix D).
In general, for small λ (thus large K) the background interaction Q dominates; the fields ϕ are semiclassically pinned to one of their minima and the system is in a disordered phase of the corresponding clock model.For large λ, instead, the tunneling and background P terms dominate and the system is in an ordered state.
Let us discuss next the behavior of the system as a function of N .For N = 2 and 3, we expect a single critical point to separate these phases.This critical point falls in the Ising and Potts universality class respectively.
For N = 4, the limit g = 0 can be examined through the mapping of the Z 4 quantum clock model into two separate copies of the Ising model in the same geometry [50] (see Appendix B for more detail).Therefore, also in this case, the system behaves as its 1D counterpart with only two gapped phases separated by a single critical point corresponding to two copies of the Ising critical point.
We additionally observe from our numerical results that for N = 2, 3, 4, the critical value of λ c is smaller than 1, due to the presence of the rung tunneling inter-action.This is a signature that, indeed, our model interpolates between one and two dimensions.In particular, our rough numerical estimates for N = 3 provide a value λ c ≈ 0.75 based on the expectation value of the order parameter O r,s .In the 1D clock model the critical value is λ (1) c = 1, whereas in the 2D system the paramagnetic and ferromagnetic phases are separated by a first-order phase transition for λ (2) c ≈ 0.498 calculated with tensor network techniques [31,67]; therefore λ c .For N = 5 (or larger), the phase diagram at g = 0 becomes richer and the second-order renormalization group analysis confirms the existence of a gapless phase between the ordered and disordered phases [see Fig. 5(b)].There is indeed a finite intermediate interval of λ, thus of the bare Luttinger parameter, such that the spin sector of the Hamiltonian ( 19) is gapped by the tunneling interaction, but the charge sector remains gapless.In the next section we analyze in detail the phase diagram for N = 5 and we verify through tensor network simulations that such gapless phase exists and extends also to finite g.

V. THE PHASE DIAGRAM AND THE ONSET OF THE COULOMB PHASE
In the following, we analyze the phase diagram for N = 5 based on density matrix renormalization group (DMRG) [73][74][75] and Wilsonian second-order renormalization group (RG).Our results are displayed in Fig. 5. Our bosonization description, however, does not rely on specific assumptions about N and qualitatively similar results hold for N > 5 as well, as shown in Appendix D.
The ladder model presents a single extended gapped phase that surrounds a gapless phase appearing around λ ∼ 0.75 for g 0.05 (depicted in green in Fig. 5).These phases are separated by a Berezinskii-Kosterlitz-Thouless (BKT) phase transition, which can be detected by evaluating the fidelity susceptibility (FS) of the system [55,68].
The gapless phase displays the properties of a Coulomb phase, characterized by an emergent U (1) symmetry, in which electric fields may propagate without mass gaps along the legs of the ladder and they appear to be only weakly screened by the dynamical matter.
The gapped phase displays instead several crossovers connecting different regimes which include a confined and a Higgs regime, as typical for Abelian LGTs [48].These crossovers are signaled by extrema in both the susceptibilities associated with the Hamiltonian terms and fidelity susceptibility.In both cases these susceptibilities do not diverge with the system size, consistently with crossovers rather than phase transitions.For simplicity, we will refer to the different gapped regimes as "phases", although they are not separated phases in the thermodynamic meaning, but rather adiabatically connected regimes, as in the case of Fradkin-Shenker LGTs (in the fundamental representation for the matter degrees of freedom) [48].I. Summary of the regimes appearing in the phase diagram as predicted by the second-order RG flow.The first table displays the phases of the clock limit g = 0 (in which θ0 is always pinned).The gapless Coulomb phase extends for finite g.The gapped phases, instead, cross over into the additional confined phases listed in the second table.Columns with white or colored background distinguish the phases of the pure LGT (λ → 0) and the ones appearing at finite λ respectively [the background colors match the regions in Fig. 5 (a)].The listed dominating interactions refer to the coupling constants in Eqs. ( 13), ( 19) and (D5).
In the following, we will discuss the origin of the different phases based on their bosonized description and we will identify the phase transitions and crossovers based on numerical simulations.Then, we will focus on the properties and identification of the thermodynamic regimes in the phase diagram which can be obtained by the analysis of the behavior of several observables over the ground states of the model, and by the screening properties that characterize the system upon the introduction of static charges.Table I offers a summary of the thermodynamic regimes appearing in the model.

A. The analysis of the Z5 gauge model in the clock limit
We begin our analysis of the Z 5 LGT from the clock limit (g → 0), which is described by the Hamiltonian (19).In the previous section, we emphasized that a simple scaling analysis predicts that the competition of the background interactions P and Q yields the onset of two fully gapped phases, separated by an additional gapless phase for N ≥ 5.In this critical phase, a gap is opened in the spin sector only, as effect of the rung tunneling T ; therefore its central charge is c = 1, and its charge sector is characterized by an emerging U (1) symmetry (θ ρ → θ ρ + α) in a way analogous to the one-dimensional quantum clock models [50].We call this phase Coulomb phase, in analogy with its higher dimensional counterparts.
According to the first-order scaling analysis of the Hamiltonian (19), the three phases alternate in the following way: for K < 4/N , P dominates and the model is fully gapped; for 4/N < K < N/ √ 2, the tunneling term dominates and only the spin sector is gapped; for K > N/ √ 2, the background Q term dominates and the gap in the charge sector is restored.In these estimates The P -dominated phase (large λ) corresponds to both the Higgs phase of the LGT and the ordered phase in the corresponding clock model.In this phase the θ fields are semiclassically pinned.
In the gapless Coulomb phase, only the field θ σ is semiclassically pinned.Finally, the Q-dominated phase would correspond to a symmetric and disordered phase in the clock model, with the ϕ fields pinned, which can be mapped into the deconfined phase of the LGT.
These simple predictions obtained with a first-order renormalization approach, however, are not sufficient to completely describe the behavior of the system.Due to the presence of non-commuting interactions, the renormalization group flow yields in general the appearance of novel effective interactions, which can be more relevant than the original terms in Eq. (19) and must be taken into account for a more rigorous study of the phase diagram.These emergent interactions appear naturally when considering higher orders of the sine-Gordon interactions.In particular, among the set of operators appearing in its second-order analysis, we focus on the following terms: In Appendix D, their derivation is presented in detail.A key feature of the second-order interactions (20) is that the term C ρ commutes with the rung tunneling.Therefore, it allows for the appearance of a new gapped phase in which both T and C ρ flow towards strong coupling, such that both ϕ ρ and θ σ can be qualitatively considered pinned to a semiclassical minimum.Therefore, for values of λ intermediate between the deconfined (small λ) and the Coulomb phase (λ ∼ 0.75), a new phase appears, which corresponds to a clock model ordered along the rungs but disordered along the two legs.We call this regime quadrupolar, since it is characterized by the condensation of pairs of mesons with opposite dipoles along the two legs (see Sec. V D).We additionally observe that the C ρ interaction is responsible for completely gapping the first-order gapless phase for N = 3 and N = 4 (see Appendix E).For N = 5, instead, the gapless phase is reduced at second order but it survives over an extended region (see Fig. 5).The four phases appearing in the Z 5 clock limit (and, in general for N ≥ 5) are summarized in the first part of Table I and their main properties are discussed in Sec.V D.
To understand the phase diagram of the model, it is useful to compare the two panels of Fig. 5, which depict the results of the numerical DMRG simulations [panel (a)] and the numerical solution of the second order RG equations in Appendix D [panel (b)].The DMRG results clearly provide a more rigorous scenario, although our numerical analysis is based on finite system sizes and an accurate scaling analysis of the phase boundaries is beyond the scope of this work.The RG results, instead, give an insight on the thermodynamical behavior of the system and provide a useful reference to compare its different regimes.They can be easily extended to N > 5 (see Appendix D), but they suffer from several approximations adopted in the bosonization procedures and their results are not reliable in the extreme regions λ 1 and λ 1.The numerical solution of the RG flow equations derived from the Hamiltonian (19) determine the g = 0 axes of the phase diagram in Fig. 5.The so-obtained phase diagram displays all four phases (deconfined, quadrupolar, Coulomb, Higgs) in its lowest part.We observe that both the DMRG simulations and the second-order RG equations confirm the appearance of the gapless Coulomb phase between the disordered and Higgs phases.The main difference between the two approaches in the limit g → 0 is about the deconfined phase.The DMRG results present a smooth behavior of the low λ region compatible with the quadrupolar phase extending everywhere but in the limit λ = g = 0.The second-order RG calculation would instead indicate the onset of an extended deconfined phase for small g and λ.
This difference is due to a limitation of our secondorder RG approach, which relies on the assumption that the bare values of the Luttinger parameters are given by 1/λ.This assumption is realistic only for λ in a neighborhood of 1, but it is likely that the bare values of K s do not diverge for λ → 0. Taking into account this limitation, the region on the extreme left of the phase diagram 5 (b) must be considered non-physical and, similarly to what we discussed for the pure LGT case, the deconfined phase shrinks to the single point g = λ = 0.
In order to understand the extension of the phase diagram for g > 0, in the following we present numerical tensor network simulations of the model and we compare them with the bosonization predictions provided by the Hamiltonian (13).

B. DMRG phase diagram
In Fig. 5(a) we show the phase diagram obtained from DMRG simulations with bond dimension up to 300.Our results confirm the presence of two extended phases: the gapless Coulomb phase appearing for intermediate values of λ and small g [green area in Fig. 5(a)], and the surrounding gapped phase, which is depicted in blue, red and yellow in Fig. 5(a), to distinguish the quadrupolar, confined rung-dominated and Higgs regimes respectively.The Coulomb phase is separated from the gapped phase by BKT phase transitions, whereas the gapped regimes are adiabatically connected by smooth crossovers.The typical DMRG truncation errors are about 10 −8 in the gapless phase and are smaller than 10 −10 in the gapped phase.
The gapless Coulomb phase is the only characterized by a logarithmic growth of the entanglement entropy as a function of the subsystem size ; therefore, it can be easily distinguished from the others.The entanglement entropy of the system in this phase follows indeed the Calabrese and Cardy formula [76] where c is central charge of the underlying conformal field theory and c α is a non-universal constant.In Fig. 6(b) we show the entanglement entropy as a function of for λ = 0.75 in the gapless phase.The central charge from the fitting is 1.02 which is consistent with the bosonization prediction c = 1, confirming that only one of the Luttinger liquid sectors remains gapless.
To verify that the phase transitions between the gapless and gapped phases are of the BKT kind, we analyze the behavior of fidelity susceptibility (FS).The FS per link is defined by where the fidelity thermodynamic value with a characteristic 1/ log(L) dependence [55,68], which distinguishes BKT from other phase transitions.In Fig. 6 we show χ F as a function of λ for different system sizes.The two peaks correspond to the transitions between the gapless Coulomb and the gapped quadrupolar and Higgs regimes.The BKT transition points are determined by these peaks and we verify in Fig. 6(a) that their finite-size scaling approximately follows the predicted logarithmic behavior.
We observe that the critical FS grows very weakly with the system size when varying the parameter g (see Fig. 7).
To characterize the transition driven by g, we consider a system with rough (smooth) boundary condition at the left (right) boundary, where the critical FS grows faster than the case with smooth-smooth boundaries.Fig. 7 shows the FS as a function of g for λ = 0.75.Based on the system size we can numerically access L ≤ 161, the data display the evolution of the FS towards the formation of a peak.In the definition of the phase diagram in Fig. 5 (a), we set the critical value of g by the position of the maximum for L = 161.
The fidelity susceptibility across the gapped phase displays additional maxima, characterized by curves smoother than the BKT behavior (compare Fig. 8 with Fig. 6).These maxima typically appear in correspondence with the maxima of other susceptibilities of the system and, in particular, we consider the susceptibilities χ τ and χ σ , defined by: where correspond to the second and the fourth terms (proportional to g and λ) in the Hamiltonian (Eq.3).The maxima of the susceptibilities χ F , χ τ (χ σ ) appear at the same value of g (λ) when scanning λ (g).We present an example in Fig. 8, where χ F and χ τ are shown as functions of g for λ = 1.8, thus across the crossover between the Higgs and confined (rung-dominated) regimes.The maxima of χ F and χ τ indicate a crossover at g ≈ 0.67.These susceptibilities, however, do not diverge when increasing the system size.The insets of Fig. 8 show the system size dependence of χ F and χ τ : both susceptibilities clearly converge to a finite value for larger and larger system sizes (see the insets in Fig. 8).
As already mentioned, χ F does not diverge even at a BKT phase transitions.However, by comparing its behavior in the insets of Fig. 6 and Fig. 8, we observe that in the crossovers within the gapped phase, χ F grows even slower with the system size.
Within the gapped phase, we can distinguish three main crossovers that correspond with the boundaries between the Higgs (yellow), quadrupolar (blue) and confined rung-dominated (red) regimes in Fig. 5 (a).In these cases we observe clear maxima in all the relevant susceptibilities and, based on the observables we will introduce in the next subsections, we can distinguish these three phases consistently with the properties listed in Ta-ble I.
In particular, the crossover between the Higgs and confined regime is qualitatively the same with respect to the analogous crossover in the Fradkin and Shenker LGTs in higher dimensions [48].As a function of g we find not only a maximum of χ F and χ g , but also a maximum in the susceptibilities associated with the plaquette energy and other observables.Analogously to other LGTs with Higgs matter [69][70][71], these features suggest that such crossover can be a Kertész line, namely a percolation phase transition in a corresponding two-dimensional classical model at finite temperature, which is not accompanied by any singularity in the thermodynamic properties of the system.
Finally, let us mention that inside the confined and quadrupolar phases, there appear an additional line of maxima of χ F (associated with either very weak local maxima or inflection points of χ λ ).We did not consider this line to be an additional crossover because of the fundamentally equal behavior of the system for λ larger and smaller than these maxima.These maxima of χ F , in particular, seem to be associated with a prolongation within the gapped phase of the BKT phase transition separating the quadrupolar and Coulomb phases [see Fig. 5 (a)].

C. The observables of the system
The results from the RG analysis of the low-energy bosonized description of the system (13) do not allow us to clearly distinguish between crossovers and phase transitions within the gapped phases.However they confirm that the gapless Coulomb phase survives at finite values of g.In the following we summarize the main results from the complete study of second-order RG equations derived from the Hamiltonian (13).The detail of their derivation is presented in Appendix D.
When considering g > 0 the two fundamental interactions which determine the properties of the system are the rung tunneling and the electric field energy.We can express them in the following form: These interactions commute and can be simultaneously minimized.When relevant, they tend to pin the combinations of fields √ 2θ σ − θ 0 , ϕσ √ 2 + ϕ 0 and ϕ ρ .The solution of the second-order RG equations shows the onset of the two additional confined phases listed in the second part of Table I for g 0.4 [red and brown phases in Fig. 5 (b)].The fully confined phase appears only when the electric field and background Q terms dominate, thus pinning all the ϕ fields.This implies, in particular, that the rung clock degrees of freedom are in a disordered state and, based on our DMRG results, this happens only in the λ = 0 limit (for any value of g > 0).As soon as λ > 0, indeed, the numerical results display ordered rung operators (σ 0 in the unitary gauge).
The different regimes appearing in the phase diagram can be characterized and distinguished based on suitable observables calculated over the ground state of the system.Therefore, in the following, we introduce several gauge-invariant correlations and string-order parameters that constitute a diagnostic toolbox to characterize the thermodynamic phases of the model, in order to be able to compare the field-theory prediction with the numerical results.
To investigate the clock limit of the system, we introduced the order parameter O r,s which allows us to distinguish the phases for g → 0. This string-order parameter can be extended to define the creation operators of the mesons in the system: The operator M s introduces a meson lying on the s leg of the system by introducing opposite dynamical charges on the sites x and y linked by an electric flux line.By considering the axial gauge, it is straightforward to derive the bosonized description of these string operators in the right hand side of Eq. ( 28) through the mapping (11).
Our former analysis of the low-energy sector of the model displayed a separation between spin and charge degrees of freedom.To this purpose it is useful to introduce also the following combination of the meson operators: where we specified their explicit form in terms of the bosonic fields.M σ creates a pair of opposite dipoles on the rungs x and y, which are connected by two opposite electric flux lines lying on the two legs of the ladder.Due to this configuration, the resulting doubled meson presents a vanishing total electric dipole, but a non-vanishing electric quadrupole.In the case of M ρ , instead, two parallel mesons are created with two negative charges in the rung x and two positive charges in the rung y, connected by parallel electric fluxes along the legs of the ladder.The meson strings are mostly useful to investigate the properties of the matter in the system: depending on their behavior at large space separation x − y one can distinguish phases in which the mesons condense (for example the Higgs phase), and phases in which the dynamical matter is screened or confined.In this respect, the quadrupolar phase is characterized by a condensation of the spin meson M σ accompanied by an exponential decay of the the charge meson M ρ .
We introduce next another family of observables which describes on one side the charge fluctuations of the system, and, on the other, the behavior of the electric field.We call these string operators t'Hooft operators since they are a generalization of the t'Hooft string G.In particular, we define them for a system displaying smooth boundary conditions on the right edge and rough boundary conditions on the left as in Fig. 1 (the extension to smooth boundaries on both sides is straightforward): The string operator G s corresponds to the exponential of the total charge laying on the leg s between the r th site and the end of the ladder.Through the Gauss law, it also corresponds to the exponential of the total electric field generated generated by these charges.It can also be considered an operator moving a magnetic flux from the right edge, along all the plaquettes of the ladder until the r th rung and then pushing it out of the ladder through the operator τ r,s .Having smooth boundary conditions on the right edge (see the boundary conditions in Appendix C), the t'Hooft operators G s (r) additionally exemplify the physical meaning of the fields ϕ s .
In analogy with the meson strings (which can be considered the dual operators of the t'Hooft strings), it is convenient to define the following t'Hooft operators addressing the spin and charge sectors: We observe, in particular, that the operator G ρ (r) = e −i 2π N Qr is local (due to the smooth right boundary conditions we have chosen) and measures the charge Q r of all the matter sites on the right of the r th links.G σ is instead related to the charge differences between the two legs.
Finally we introduce the two-point correlation function of the rung tunneling operators: Its expectation value gives us information about the gap opened by the rung tunneling term, corresponding to the T interaction in Eq.( 13) and it provides a direct evidence of the fact that both the fully confined and deconfined phases do not extend for finite values of λ.

D. Features of the thermodynamic phases
After defining the gauge-invariant observables of the system, we can proceed with the study of the phases and regimes we observe in the numerical simulation, summarized in the phase diagram Fig. 5 (a).In particular, our main findings are that for small values of g, the Coulomb phase is stable for g 0.05, whereas both the Higgs and quadrupolar phase smoothly cross over towards the confined-rung dominated regime.
From the field theoretical description, we can predict the behavior of the two-point and string correlation functions introduced above.In particular, when a bosonic field, for example θ σ , is semiclassically pinned by a relevant interaction, its fluctuations are suppressed and its correlation are approximately constant across the system; the fluctuations of its dual field, ϕ σ , are instead maximal, yielding an exponential decay of the two-point correlation functions of the related vertex operators, thus of G σ in the example [see Eq. ( 33)].
Our bosonization and RG analysis predicts that the gapless Coulomb phase is characterized by a gap in the spin σ and rung 0 sectors, whereas the charge ρ sector is gapless.This implies that the expectation value of both the meson M ρ and the t'Hooft operator G ρ decay as a power law.In Fig. 9 we display the typical behavior of these two observables in the ground state of the system within the Coulomb phase.In particular, we considered a ladder with smooth boundary conditions on both sides of the system.This implies Dirichlet boundary conditions for ϕ σ and Neumann boundary conditions for θ σ at the two edges.Consequently one obtains the following leading behaviors: where we introduced the chord distance d and a modified chord distance d (see [72] for detailed calculations): These analytical predictions are roughly compatible with the numerical results in Fig. 9: the t'Hooft operator G ρ decays as a power law with its distance from the edges of the system and, for smooth boundary conditions, it is symmetric under space inversion; the meson M ρ approximately decays as a power law of the modified chord distance between its charges, although it presents a bended shape in the logarithmic plot 9 (a), most probably caused by subleading terms in its bosonized description (in analogy with bosonic systems [72]).Strong subleading deviations from the predictions in Eqs.(36,37) are observed approaching the BKT transitions.
Concerning the gapped thermodynamical phase we begin our analysis by observing that the numerical results are consistent with having an approximately constant expectation value of the rung two-point correlation function R(x, y) for any λ > 0 (see Fig. 10).In the limit λ → 0, our findings suggest that R maintains its distanceindependent behavior, although its value considerably decreases and is smaller than 10 −8 when λ 0.01.As previously mentioned, the constant behavior of R indicates that the fully confined and deconfined regimes exist only in the pure lattice gauge theory limit λ = 0, whereas the gapped phase for any finite λ > 0 falls in one of the following three regimes: quadrupolar, confined rungdominated and Higgs [see Table I and Fig. 5(a)].This is consistent with the fact that the Luttinger parameter approximation K σ/ρ ∼ 1/λ breaks for small values of λ and the bare values of the Luttinger parameter must be considered bounded.In particular, our findings suggest that the rung tunneling term T is relevant for all values λ > 0.
To distinguish the quadrupolar, confined and Higgs gapped regimes, we can compare the behavior of the meson strings M ρ and M σ and the t'Hooft parameter G ρ .As summarized in Table I these three regimes are identified, from the RG analysis, by different behaviors of the gauge-invariant observables, which we discuss in the following.
The confined (rung-dominated) regime appears for g 0.4 (red region in Fig. 5).In this phase both the electric field energy and the interleg tunneling flow to strong coupling, thus pinning the field combinations ϕ ρ , √ 2θ σ − θ 0 , ϕ σ + √ 2ϕ 0 , as can be derived by the Hamiltonian (13) [see also the action (D5)].Consequently G ρ is constant, whereas all the mesons display an exponential decay (see the red curves in Fig. 10 for a typical scenario) and, in this region, we have a good agreement between the RG and DMRG predictions.
The picture is more complicated in the gapped phase at low values of g.In the quadrupolar phase the bosonization analysis yields that the interactions T and C ρ pin the fields θ σ and ϕ ρ , whereas θ 0 is pinned by the background P 0 interaction.Therefore, based on RG, the meson string M σ should display a constant behavior in this regime; this implies that pairs of mesonic strings with opposite dipoles at each end condense.This kind of mesons are indeed compatible with having ordered rung operators.
However, what is observed by the DMRG results is that this picture captures the behavior of the system only for g = 0 [blue limit in Fig. 11 (a)].For small but finite values of g, instead, M σ always displays a weak exponential decay.Nevertheless, for g 0.2 the corresponding decay length ξ Mσ is larger than the typical system sizes we can probe (L ∼ 100) and it diverges for g → 0. This smooth evolution from a constant value at g = 0 to a weak exponential decay is the first signature of the crossover between the quadrupolar and confined phase.We additionally observe that ξ Mσ transitions from a behavior proportional to g −4 to a slower decay for g increasing towards the confined phase [see Fig. 11 (b)].
Concerning the meson M ρ , its expectation value clearly decays exponentially with a very short decay length for any value of g in the quadrupolar phase.This can be seen by considering the data for λ 0.7 (blue curves) in Fig. 12, which depicts the observables of the system for g = 0.2 as a function of λ.This is consistent with the fact that this meson string does not commute with the relevant and ordered T interaction.Finally, G ρ is approximately constant, showing that also in this phase the fluctuations of the dynamical charge are suppressed [blue curves in Fig. 12 (a)].
Concerning the Higgs regime, the features of the crossover into the confined phase for growing g are even stronger.The bosonization prediction is that the tunneling and background P terms pin all the θ fields to their semiclassical minima.This implies that all the mesons should condense in this regime and all the M strings should display an approximately constant expectation value.On the contrary, the fluctuations of the charge are maximal and G ρ must decay exponentially from the boundary of the system.What we observe in the DMRG is that, once again, the bosonization predictions are accurate in the clock limit.For any g > 0, instead, all the meson strings acquire a weak exponential decay, with a decay length diverging as g −4 for g → 0 (see Fig. 13), compatibly with results from a quasiadiabatic continuation estimating: See Appendix F for detail.Eq. ( 40) relies on g 1 in general, and we display its result for small values of g in Fig. 13 (b) and (c) (orange lines).For large values of g and λ, instead, the approximate behavior of the mesons can be deduced by approximating the ground state with a product state that minimizes the electric field and tunneling energies only.The results of this approximation are depicted as the green lines in Fig. 13 (b) and (c) and provide a reasonable estimate of the behavior of the meson even for the chosen value λ = 1.4 and g 1.The smooth crossover between these behaviors of the de-cay length is a further signature of the Higgs / confined crossover.
The Higgs / confined crossover has an even more interesting effect over G ρ .In systems with smooth boundary conditions, G ρ decays exponentially away from the edge as predicted for the clock limit.However, for finite g, this exponential decay stops at a certain distance * from the edge, and G ρ stabilizes to a bulk constant (see Fig. 13, showing the results for λ = 1.4).In the limit g → 0, * increases and the bulk region shrinks with smaller and smaller bulk values of G ρ .For g approaching the crossover line, instead, G ρ becomes essentially flat, corresponding to * 1.
We finally comment on the regime at very low λ.As we previously mentioned, the two-point correlation function R displays approximately a constant behavior for any λ > 0, and this is the main reason for which we consider the regions at small λ in the phase diagram 5 (a) to belong to the quadrupolar and confined rungdominated phase, rather than the deconfined and fully confined phases.However, based on the DMRG data, we observe a very weak exponential decay of R(x, y) for the regions at very small λ: the corresponding decay length is typically ξ R > 10 4 for λ = 0.02.This tiny decay may be interpreted as an effective crossover from the fully confined phase of the pure LGT to the quadrupolar and confined rung-dominated phases.At the level of the field theoretical description, this may also be due to the weak mixing between the spin σ and gauge 0 sectors of the theory which appears at second order in the perturbative RG analysis and we neglected in our flow equations (see App. D).

E. Static charges and screening
When introducing dynamical charges in the system by considering λ > 0, the confinement of the pure lattice gauge theory is disrupted in general by screening.In particular, when we insert two opposite static charges through the violation of the Gauss law in two specific sites at a distance R, as done in Sec.III, dynamical charges can accumulate in a screening cloud around them and, in general, they will suppress the electric field propagation in the intermediate region.
To roughly estimate the extent of this screening mechanism, we can compare the electric field string energy T R (see Fig. 2) and the mass of a pair of dynamical charges 2m = 4(1 − cos 2π/N )/λ.We define the length scale R * such that T (g)R * = 2m(λ); R * provides a rough approximation of the size of the screening clouds of the dynamical charges around the static ones.When T R < 2m, thus R < R * , the energy cost of an electric field string connecting the static charges is smaller than the mass cost required to screen them.Therefore ∆E increases linearly with R as in the pure lattice gauge theory, and the corresponding states do not display a complete screening of the static charges such that electric lines clearly propagate between them [see, for example, Fig. 14(b)].When R > R * , instead, screening dominates, and ∆E stabilizes towards an asymptotic value (typically smaller than 2m).In this case, localized clouds of dynamical particles completely screen the static charges and the electric field rapidly decays away from them [see, for example, the behavior in Fig. 14(c)].
In Fig. 15 we display the behavior of ∆E as a function of the static charge distance R within the different regions of the phase diagram.In the quadrupolar phase (blue) the string tension is small whereas the mass of the Higgs charges is large, therefore screening does not occur on the length scales here presented.R * is indeed large, and, as expected, for R < R * the electric field propagates along both legs in the intermediate regime [see Fig. 14(b)].In this regime, in particular, the expectation value of the electric field E decreases only weakly away from the static charges and it is approximately constant along the leg connecting them [see Fig. 16 (b)].The electric field along the rungs is typically much smaller than the one along the legs, consistently with the rung tunneling being a relevant interaction.
By increasing g, the system smoothly evolves into the confined rung-dominated phase [see Fig. 15(d)]; in this regime the string tension is stronger, thus R * decreases.For R < R * the state is again analogous to the confined limit and, depending on g, the electric field either propagates on both legs (intermediate values of g) or only in the leg of the static charges (large g).In the former case, we observe that the electric field predominantly flows from one leg to the other along the same rungs of the static charges.The case R > R * is exemplified instead by Fig. 14(d) and 16(d): the electric field is exponentially suppressed away from the static charges by a screening cloud of dynamical charges.
In the Higgs phase, screening is extremely strong such that R * is typically around one [see Fig. 15(c)].Also in this case the electric field is exponentially suppressed with the distance from the static charges [Figures 14(c) and 16(c)].In this regime the amplitude | τ | is in general very small due to the strong plaquette and tunneling interactions, thus concurring in suppressing further the expectation values E .
Finally, in the Coulomb phase, we observe that R * is very large due to the weak string tension; thus, the static charges display confinement over long distances [see Fig. 15(a)].In this gapless phase, the electric field E decays weakly as a power law away from the static charges [Figures 14(a) and 16(a)] and it propagates on both legs, consistently with the predicted algebraic decay of the correlation functions in the charge ρ sector.

VI. EXTENSION TO MULTIPLE LEGS
The field theoretical approach we adopted for the analysis of the system can be extended to investigate wider ladders with a finite number of legs L y .To this purpose, we can apply a so-called coupled-wire construction (see, for example, the review [46] and references therein): we decompose the system (in the axial gauge) into a set of Luttinger liquids that describe each horizontal stripe of the lattice and interact with each other.In particular, we consider the axial gauge and we generalize the matter bosonic fields ϕ s and θ s in Eq. ( 11) into pairs of dual fields ϕ y and θ y labeled by the coordinate y = 1, . . ., L y which specifies the row they refer to.In a similar way, the fields ϕ 0 and θ 0 in Eq. ( 12) are extended to ϕ 0,y and θ 0,y (with y = 1, . . ., L y − 1), in order to represent the gauge E decays approximately exponentially for r < R/2.The four panels correspond to the same coupling constants chosen in Fig. 14.
bosons on the vertical links between the matter rows y and y + 1.As before, each pair of dual bosonic fields requires the addition of a pair of background interactions of the P and Q kinds.
The major difference between two and multiple legs relies in the form of the electric field and tunneling interactions; Eq. ( 27) takes a more symmetric form: The electric field energy on each horizontal link can indeed be described by a combination of the difference of electric fluxes ingoing and outgoing from the vertical links and the total Higgs charge on the same row.These interactions must be supported by suitable boundary conditions.
This form of the coupling between subsequent rows allows for the existence of a Coulomb phase also for multiple legs (with smooth boundaries at y = 1 and y = L y ).Such phase can be understood by the emergence of a gapless sector characterized by the dual fields: These fields are a linear superposition of all the matter fields and they generalize the charge sector into a bulk mode of the system; the operator e iϕ bulk , for example, is linked to the motion of a magnetic flux along the vertical direction of the lattice.Analogously to its two-leg counterpart, this sector remains gapless for intermediate values of λ (such that the related background terms are irrelevant) and small g (such that the G term in ( 41) is irrelevant and does not open a gap in this sector).
A key aspect of the interactions ( 41) is that the electric field energy involves three bosonic fields in the multileg case (differently from the corresponding equation (C11) for the ladder).This implies that its scaling dimension increases and its growth in the RG flow is slightly suppressed.This has the important implication of increasing the extension of the Coulomb phase to larger values of g.The extension of the gapless phase is further increased due to the second-order terms generated by the G interactions being suppressed as well (the generalization of the C ρ terms in (D5) in particular).
Concerning the further extension to a fully twodimensional system, the analysis becomes qualitatively more complex: an extended deconfined and topological phase appears for small values of g and λ, whereas the gapless phase is supposed to evolve into a U (1) symmetric weakly gapped phase [77].The mechanism gapping the critical phase is non-trivial and can be understood based on the mapping, for g = 0, into the 2D Z N quantum clock model and its classical 3D analog.In three dimensions, it is indeed known that the classical clock model does not display an extended gapless phase due to the P background terms being "dangerously irrelevant" perturbations and the system displays only a single phase transition in the 3D XY model universality class.This is due to the fact that the irrelevant P operators become relevant when approaching the U (1) symmetric gapless Nambu-Goldstone fixed point [78], thus causing a second step in the RG flow towards the ordered ferromagnetic phase of the clock model (see the phase diagrams evaluated in [79,80]).As a consequence, the system is characterized by two different length scales ξ < ξ , both diverging at the critical point, which are associated with the onset of U (1) and Z N symmetric features respectively.This causes indeed the appearance of a crossover between a U(1) symmetric and a Z N symmetric regime in the ordered phase.This analysis has been numerically well-verified for the classical 3D model (see, for example, [77,79]), and it has been recently confirmed also for the quantum 2D case [80], thus leading to the conclusion that no gapless phase exists for the clock limit g → 0 and the dual pure LGT λ → 0 [77].

VII. CONCLUSIONS
Ladder setups offer the simplest realization of a lattice gauge theory whose dynamics crucially relies on the plaquette interactions.In the path towards the experimental analog quantum simulations of gauge theories, therefore, the realization of gauge models in the ladder geometry would constitute an important milestone bridging one-dimensional chains and higher dimensional setups.The first steps in this direction have already been ac-complished in ultracold atom systems trapped in optical lattices: a recent experiment [11] has proved that a tunneling term mediated by an effective Z 2 gauge degrees of freedom can be realized based on density-dependent laser-assisted tunneling techniques.This is indeed the rung tunneling interaction required for the realization of a gauge theory in the ladder based on the axial gauge [81].With respect to the double-well systems proposed in Ref. [81], our model additionally includes the plaquette tinteraction, which, in the axial gauge, must be engineered as an operator acting on two neighboring rung double-wells.
Previous works focused on the analysis of several aspects of dynamical gauge theories in ladders [33,37,51,81].In this article, we explored the general features of the full Kogut-Susskind model with a Z N LGT and Higgs matter degrees of freedom.We analyzed its phase diagram based on both a low-energy field theoretical description inspired by bosonization, and DMRG numerical simulations.The model displays different features for N ≤ 4, where a single critical point is observed in the limit of strong plaquette interactions (g → 0), and N ≥ 5, where instead an extended critical Coulomb phase, with emergent U(1) symmetry appears.
Our numerical analysis focused mostly on the N = 5 case, but, based on the renormalization group study of the field theory description of the model, we conclude that the existence of this extended gapless phase is stable for larger (but finite) values of N , in which the gapless phase covers broader and broader domains in the λ coupling constant.
On the technical side, our bosonization description of this gauge theory can be applied, in general, to describe a broad family of quantum clock models and our description can be extended to the study of more complicated quasi-one-dimensional ribbon geometries.In this appendix we provide further information about the pure lattice gauge theory and the calculation of its electric string tension.
Eq. (8) shows that the pure lattice gauge theory in the axial gauge is equivalent to a Z N quantum clock model with the addition of non-local disorder operators.This Hamiltonian can be reduced in a completely local form by use of a unitary bond-algebraic duality transformation, as presented in [1,50]: This mapping preserves indeed all the commutation relations and maps the disorder operators into a local term.The dual Hamiltonian reads: This dual Hamiltonian describes the p-clock model in a longitudinal field, with the important feature that τ1 does not appear in the dual model.As a consequence, the global symmetry G is mapped into the holographic symmetry defined by G = σ1 , such that [σ 1 , H] = 0 (see [59] for more detail on holographic symmetries in Z N models with Higgs matter).From Eq. (A4) we see that the ground state is degenerate in the limit g → 0, because the first clock degree of freedom may assume any of its N states.This corresponds to the ordered ferromagnetic limit for the Hamiltonian (8).For any value g > 0, however, the degeneracy is split by the term proportional to G entering the Hamiltonian.The ground state is then symmetric and non-degenerate for any g > 0 and the gap separating it from the first excited states grows linearly with g.For N = 2 this corresponds indeed to the Ising model in a longitudinal and transverse field which presents only one gapped phase [57,58].
This feature distinguishes our ladder model from its two-dimensional counterpart, since, for pure Z N lattice gauge theories in 2D square lattices, there exists a phase transition between a confined (g > g c ) and a deconfined (g < g c ) phase for a finite g c > 0 [1,48].
To verify that the extended gapped phase characterizing the pure lattice gauge model (8) at g > 0 corresponds to a confined phase, we consider the behavior of the system in the presence of external static charges and we estimate their string tension.
In order to simplify our analysis, we modify the rough boundary on the left of the ladder by including an additional 3-site boundary plaquette term which does not violate any of the gauge constraints: This magnetic field term must be added to the boundary term in Eq. ( 4) and its inclusion explicitly breaks the G symmetry, favouring a symmetry-broken ferromagnetic ground state with all the operators σ 0 aligned along their eigenvalue 1 in the axial gauge choice of Eq. (8).
A confined phase is characterized by the presence of a linear string tension for the electric flux lines.When introducing two opposite static charges, the energy of the ground state displays a linear dependence on the distance R between them.The static charges are introduced by modifying the gauge constraints (7) associated to two arbitrary ladder sites.In particular, we impose that the physical Hilbert space fulfills G x,↓ |ψ = e i 2π N |ψ and G y,↓ |ψ = e −i 2π N |ψ .For simplicity, we introduce the charges in the same leg and we consider y > x.In this physical subspace, the Hamiltonian in the axial gauge takes the form

(A6)
To verify that this Hamiltonian only supports a confined phase for any g > 0, we consider the behavior in the two limits g → 0, ∞.For g → ∞ the ground state is a (paramagnetic) product state with all sites obeying τ r |ψ axial = |ψ axial .The electric energy of this product state with the two static charges is where ∆E describes the energy difference between the ground states with and without the static charges separated by the distance R = |x − y| and T denotes the string tension (see Fig. 2).The phase in this limit is therefore confined.Including also the plaquette interaction, for g 1 the system still supports a confined phase, and from perturbation theory one finds (A8) In the other limit, g → 0, the ground state is instead a product state with all sites obeying σ r |ψ axial = |ψ axial .Exactly at the g = 0 limit, the ground state is deconfined since ∆E(g → 0) = 0.For small values of g > 0 we apply a standard non-degenerate perturbation theory and the lowest order correction to the energy yields ∆E(g 1) = 2g 3 R + O(g 7 ).
(A9) Also in this case, the static charge energy presents a linear dependence with respect to their distance R, thus showing that that the deconfined phase is unstable under any small g perturbation (see Fig. 2).In this respect, the ladder model behaves like a fully 1D system and the ground state of the pure lattice gauge theory displays confinement of the static charges for any g > 0.
Appendix B: The case N = 4 in the limit g = 0 In the limit g = 0, the lattice gauge theory model is equivalent to the quantum clock model ( 9) in the ladder geometry.In the specific case N = 4, we can apply a unitary mapping from the Z 4 clock operators into two pairs of Z 2 Ising operators [50].We introduce an additional index j = 1, 2 to label these two sets of Pauli matrices: where σ a label the Pauli matrices.Based on this unitary mapping, the model of Eq. ( 9) for N = 4 becomes: The resulting Hamiltonian corresponds to two copies j = 1, 2 of the Hamiltonian (9) for N = 2 (up to an overall rescaling of the energy by a factor 1/2). Therefore the critical value of λ at g = 0 is the same for N = 2 and N = 4.

Appendix C: Details about bosonization
In this appendix, we present the detail about the construction of the effective Hamiltonian (13).
The fundamental criterion to built a low-energy description in continuum space of the Z N LGT on the ladder is to create a mapping from the clock operators to the bosonic fields θ s and ϕ s which preserves their algebraic properties.To this purpose, we verify first that Eqs.(10,11,12) fulfill the commutation relation in Eq. ( 2): ζ r,s η r ,s → e −iθs(r) e −iϕ s (r )+iϕ s (r +a) = e −iϕ s (r )+iϕ s (r +a) e −iθs(r) e −[θs(r),ϕ s (r )−ϕ s (r +a)] = e −iϕ s (r )+iϕ s (r +a) e −iθs(r) e i 2π N [Θ(r−r )−Θ(r−r −a)]δ s,s = e −iϕ s (r )+iϕ s (r +a) e −iθs(r) e i 2π N δ r,r δ s,s → η r ,s ζ r,s e i 2π N δ r,r δ s,s , (C1) which verifies Eq. ( 2).In the previous calculation we used the Campbell-Baker-Haussdorf formula and we adopted the convention that Θ(x) = 1 for x ≥ 0. The analogous property in Eq. ( 1) can be verified in the same way.
Concerning finite systems, we observe that for the right edge, characterized by smooth boundary conditions (see Fig. 1), the definition of η L,s and τ L,0 must be taken with Dirichlet boundary conditions ϕ s (L + a) = 0 such that: By introducing different boundary conditions for ϕ ↑/↓ at the two edges it is possible to add a background electric field thus modifying the θ vacuum of the theory.Based on the mapping (11,12), we are now ready to derive the effective Hamiltonian (13).We list in the following all the Hamiltonian terms in the axial gauge and their continuum approximations.The intraleg tunneling terms read: where we considered that the bosonic fields vary slowly over the length scale a, and we neglected constant contributions.This term clearly contributes to the first line of Eq. ( 13) for s = ↑, ↓.An analogous contribution, for s = 0, is obtained from the plaquette term: The additional quadratic terms of the fields ϕ s are derived from the mass and rung electric field contributions.
The former reads: for s = ↑, ↓.The rung electric energy has an analogous form: The sum of these four quadratic terms determines the Luttinger liquid part of the Hamiltonian (13) with the parameters: ) These values of the Luttinger parameters provide the bare values entering the renormalization group flow, whereas the velocities are equal in all the Luttinger sectors and are invariant through RG due to the emergent Lorentz symmetry of ( 13).
The interacting terms in (13) are determined by the interleg tunneling and the leg electric field term.The former is straightforwardly obtained by the mapping (11,12): (C9) The latter must be estimated by considering its string operator formulation in the axial gauge, which is derived from the form of the gauge constraints (in the case of a smooth right edge).In particular, on the physical states, we have: This expression is derived from Eqs. (5,6,7) similarly to the pure LGT case in Eq. ( 8).Based on Eqs.(11,12,C2) we obtain: (C11) The analogous expression for the lower leg completes the electric energy terms in (13).
Finally, we must consider the background terms.They are meant to restore the original Z N symmetry of the system, and, indeed, when translated back to the clock operators, the background terms become proportional to the identity because they correspond to the N th power of the clock operators.Despite this, they play a crucial role in determining the phase diagram of the system and their interplay is crucial to understand the transitions from gapped to gapless phases for N > 4, in a way similar to the 1D quantum clock model.
The coupling constants P s and Q s must be determined by comparing the energy spectra of the field theory (13) with the lattice model.This task is non-trivial but a rough estimate of their value can be obtained by considering the limiting cases λ → 0, ∞ and g → 0, ∞ and neglecting all the irrelevant terms.Let us consider, in particular, the clock model limit with frozen gauge boson degrees of freedom.For simplicity, we also neglect the coupling between the two legs of the ladder, and we focus in the following on a single 1D clock chain.
In the case λ → ∞, thus K ↑ , K ↓ → 0, only the P ↑,↓ terms are relevant and we neglect all the other interactions.In this case we expect that the corresponding clock model is deep in its ferromagnetic phase, where the elementary excitation is provided by the domain walls with mass M cl = 2λ (1 − cos 2π/N ).We compare this mass with the mass of a kink in the classical and static sine-Gordon model with quadratic part corresponding to the one in the Hamiltonian (13).This classical mass can be obtained by following standard techniques (see, for example, Chap.16 in [82]): By comparing the two masses we derive: which provides Eq. ( 15) through K = 1/λ for s = ↑, ↓.This expression sets the bare coupling constant of the θ ↑ and θ ↓ background terms and it must be considered as an approximation valid for small K s .The result for the small λ limit, corresponding to the paramagnetic phase of the clock model, can be easily retrieved through the duality θ ↔ ϕ and K ↔ K −1 .It results into Eq.( 16).The calculation of P 0 and Q 0 follows the same procedure.In this case, though, it is convenient to consider first the pure lattice gauge theory limit for g → 0. In this way, the plaquette term plays the role of the ferromagnetic coupling of the clock model ( 8) and one obtains Eq. ( 15) for K 0 = 1/g.Finally, the duality K 0 ↔ K −1 0 allows for the estimate of Q 0 .
Appendix D: The second-order renormalization group equations Our renormalization group analysis is based on Wilson's approach and, in particular, on a second-order perturbative calculation in momentum space of the Euclidean action corresponding to the Hamiltonian (13).
In this Appendix we summarize the main steps for the derivation of the second-order RG equations we adopted in the study of the phase diagram, and, in particular, we focus on the onset of the main effective interactions that are generated through the flow of the sine-Gordon terms in (13).
Our perturbative approach relies on considering all the interaction terms in (13) as perturbations of the Gaussian action S 0 corresponding to the quadratic terms in the bosonic fields.The quadratic action can be written as a function of either the ϕ or the θ fields and, with the former choice, it reads: (D1) where we consider a two-dimensional Euclidean spacetime.In the following, we will use both the charge and spin degrees of freedom, and the spin ↑ and ↓ fields, depending on the most convenient notation.
Based on Wilson's prescription, we distinguish fast and slow oscillating modes for each of the bosonic fields.In particular, we introduce an effective cutoff Λ in momentum space such that the slow modes are characterized by k < Λ, whereas the fast modes are defined by the choice Λ < k < Λ, with Λ = 2π/a being the ultraviolet cutoff of the system.To perform the Wilsonian RG, we will integrate out the fast modes of each field, and, in particular, we are interested in the limit Λ/ Λ = 1 + dl, with dl infinitesimal.The decomposition of the bosonic fields in fast and slow modes reads: ϕ s (x, τ ) = ϕ s,s (x, τ ) + ϕ f,s (x, τ ) , (D2) θ s (x, τ ) = θ s,s (x, τ ) + θ f,s (x, τ ) . (D3) To derive the RG equations, we will define an effective action for the slow modes in the form: where the average values are taken over the fast oscillating modes.The interacting part S I of the action matches the interacting part of the Hamiltonian (13).S I , however, collects also many effective interactions whose bare coupling constants vanish, but acquire non-trivial values during the RG flow.As discussed in the main text, we keep track only of the most relevant of these terms appearing at second order of perturbation, and we neglect simple powers of the terms appearing in the bare Hamiltonian which do not qualitatively affect the RG flow.We list here, for reference, the main interaction terms we consider: (D5) The last two integrals, in particular, refer to interactions that appear only at second order in perturbation theory, but bear important implications for identifying the physical regimes of the system.The C interactions, in particular, appear only when g > 0, differently from the C interactions which influence the system in the clock limit as well.We emphasize that these interactions are only a small subset of all the terms appearing at second order.Let us consider, for example, the operators appearing in the clock limit.Beyond the interactions (20), one should consider similar interactions in the θ fields.We point out, however, that such interactions would be in general less relevant than the P and T terms and commute with them, in such a way that their effect in determining the phase diagram is marginal, as we verified through a numerical solution of the second order RG equations including also these additional terms.
Before proceeding in the evaluation of the main second order terms, we summarize here some of the properties of the bosonic fields we will utilize.Concerning their duality relations in Euclidean time τ = it, we have: Concerning the correlation functions of the fast modes, we will adopt the following approximations: ϕ f,j (x 1 , τ 1 )ϕ f,j (x 2 , τ 2 ) f ≈ C(r) N K j ln Λ Λ , (D8) θ f,j (x 1 , τ 1 )θ f,j (x 2 , τ 2 ) f ≈ K j C(r) N ln Λ Λ .

(D10)
Here the logarithm captures the dominant scaling behavior, whereas C(r) is a function of r = v 2 (τ 1 − τ 2 ) 2 + (x 1 − x 2 ) 2 .In the following we will consider C(r) to be suitably short-ranged; for a sharp momentum cutoff, C(r) ≈ J 0 (Λr) and this assumption is not satisfactorily fulfilled; however, C(r) can be made sufficiently short-ranged with more refined cutoffs [83,84].
The first-order contribution A of the interacting action (D4) provides the standard dependence from the scaling dimensions of the RG equations.We focus in the following in the second-order contributions and, in particular, on the non-trivial terms appearing in B.
We begin our analysis by studying a part of B which appears already in the clock limit (g = 0) and we consider, in particular, the following terms: B ⊃ d 2 r 1 d 2 r 2 Q 2 q,q =↑,↓ cos N (ϕ q,s (r 1 ) + ϕ q,f (r 1 )) cos N (ϕ q ,s (r 2 ) + ϕ q ,f (r 2 )) f q,q =↑,↓ µ,µ =±1 e iN µ(ϕq,s(r1)+ϕ q,f (r1)) e iN µ (ϕ q ,s (r2)+ϕ q ,f (r2)) q,q =↑,↓ µ,µ =±1 e iN (µϕq,s(r1)+µ ϕ q ,s (r2)) e iN (µϕq,f(r1)+µ ϕ q ,f (r2)) In this expression: (i) the terms with q = q and µ = µ return the C ρ interaction; (ii) the terms with q = q and µ = −µ return the C σ term; (iii) the terms with q = q and µ = −µ provide a correction to the quadratic part of the action which we must evaluate to obtain the RG equations for the Luttinger parameters.The last terms, q = q and µ = µ , result instead in highly irrelevant interactions which we neglect.We consider the contributions (i) first.We obtain: e iN (ϕ ↑,s (r1)+ϕ ↓,s (r2)) e iN (ϕ ↑,f (r1)+ϕ ↓,f (r2)) where we used Eqs.(D7,D8), we decomposed the fields in the charge and spin sectors, we imposed Λ = Λ(1−dl), and we applied a general rescaling of the coordinates d 2 r = (1 + 2dl)d 2 r .The coordinate dependence of the correlation C has been suppressed for ease of notation.Among the terms in the previous expression, only the ones proportional to C contribute to the second-order correction of the action.The others are simplified by the analogous terms in A 2 in Eq. (D4).In this expression, we approximate C(r) ≈ a 2 /v δ(r).This results in the following contribution to the renormalized action: (D13) By following the same approach, in the case (ii), we ob-tain instead: The terms (iii) in Eq. (D11) provide a paradigmatic example of the feedback of the interactions in the definition of the Luttinger parameters.For q = q and µ = −µ adopt a two-step RG approach, where the flow is divided into separate parts, where each part is terminated when a coupling constant reaches a suitable upper threshold, which indicates when a given interaction semiclassically pins the related fields.After each separate flow, an effective Hamiltonian for the remaining unpinned sectors is considered.For K < N/ √ 2, where the rung tunneling is the most relevant term, θ σ is the first field being pinned to an energy minimum; hence, the effective Hamiltonian of the charge sector for the second RG step results where P is a suitable renormalized parameter deriving from the background P term.We observe that, after the initial flow pins θ σ , the scaling dimension of this background term decreases by a factor of 2, D P = KN/4.Based on a scaling analysis of the second step Hamilto-nian H step 2 , the extended gapless phase disappears for N = 3.
In the above analysis the background Q term was neglected, since we considered first-order contributions only.By following the two-step approach adopted in [85], we refine H step 2 by including the most relevant second-order term, which matches the C ρ interaction in Eq. (D5):

Figure 1 .
Figure 1.Schematic representation of the terms in the gaugeinvariant Hamiltonian(3).Squares and circles refer to Higgs matter and gauge boson degrees of freedom lying, respectively, on the vertices and links of the ladder.The yellow circle depicts the electric field energy, the orange square the matter mass, the green link represents a rung tunneling term and the blue plaquettes displays the magnetic plaquette interaction.Finally, in red, we depicted one of the gauge constraints (dashed ellipse).The dashed purple rectangle depicts the first unit cell on the left (rough boundary); the boundary conditions on the right, instead, are of the smooth kind.

Figure 2 .
Figure 2. String tension T of two static charges on the ladder model in the pure Z5 lattice gauge theory.T is estimated by evaluating the ground state energy difference ∆E of the system with and without static charges.∆E behaves linearly with the separation R of the static charges for each g > 0, as shown in the insets for g = 0.02 (a) and g = 10 (b) (see Eqs. (A7) and (A8)).

Figure 3 .
Figure 3. Behavior of the electric field in the pure Z5 lattice gauge theory (λ → 0) for (a) g = 1.5 and (b) g = 0.3 in the presence of two opposite static charges.The static charges are indicated by the circles and are introduced in the central region of a ladder of size L = 41 with smooth boundaries.The thickness of the ladder edges indicate the expectation value of the electric field E in arbitrary units.Some example of their value is reported in the blue labels.

Figure 4 .
Figure 4. Averaged order parameters r O r,↑ /L as functions of λ for systems in the clock limit (g = 0) for N = 3, 4 and 5.The system has rough and smooth boundaries at the left and the right sides respectively.For N = 3 and 4, the dashed lines indicate the proposed critical values of λ.For N = 5, the gray region represents the gapless phase determined by the fidelity susceptibility in a system of smooth-smooth boundaries (see Fig.6).

Figure 5 .
Figure 5. Phase diagram from DRMG (a) and second-order RG (b).In panel (a), a gapless (green) and a gapped (other colors combined) phases are separated by a BKT phase transition (black dots).The gapped phase is further distinguished into three "phases", the quadrupolar phase (blue), the Higgs phase (yellow), and the confined rung-dominated phase (red), separated by crossovers.The cyan dots indicate finite peaks in the fidelity susceptibility of the system.The inset zooms in the Coulomb phase.In panel (b), different phases are shown by different colors, including the deconfined phase (purple), the quadrupolar phase (blue), the Coulomb phase (green), the Higgs phase (yellow), the fully confined phase (brown), and the confined rung-dominated phase (red).See Table.I for their properties.

Figure 6 .
Figure 6.Fidelity susceptibility per link χF as a function of λ for systems in the clock limit (g = 0) with smooth-smooth boundary conditions.(a) Finite size scaling of the peak values of χ for the left (blue) peak and the right (red) peak.(b) Entanglement entropy as a function of for λ = 0.75 (in the gapless phase), for L = 101.The blue dots are from DMRG and the red curve is the fitting by the Calabrese and Cardy formula.

Figure 7 .
Figure 7. Fidelity susceptibility per link χF as a function of g for λ = 0.75 for systems with rough-smooth boundary conditions.The peak in χF indicates a transition between the gapless and the gapped phases.

Figure 8 .
Figure 8. Fidelity susceptibility (dots), χF , and the susceptibility of τ (crosses), χτ , as functions of g for L = 21, 41, 81, 161 and 321, for λ = 1.8.χτ of different L are indistinguishable in the scale of the figure.The insets shows the finite-size analysis of the peak values of χF and χτ .

Figure 9 .
Figure 9. Mρ(x1, x2) and Gρ(r) in the gapless Coulomb phase as functions of d(x1 − x2|2L) (with x2 = 2) and d(r|L) in log-log scale.The system is a ladder of length L = 101 with smooth boundary conditions for g = 0.001 and λ = 0.75.The black lines are the results of fits with a power-law decay.

Figure 10 .
Figure 10.Examples of R(x1, x2), Gρ(r), Mρ(x1, x2) and Mσ(x1, x2) in semi-log scale for each regime of the gapped phase.The related values of the coupling constants are indicated by the stars with the corresponding colors in Fig. 5.The system is a L = 81 ladder with smooth boundaries.x1 is chosen to be 2.

Figure 11 .
Figure 11.(a) Meson Mσ for λ = as a function of its length.Different colors represent different values of g.Mσ is a constant in the limit g → 0 and decays exponentially for g > 0. (b) Decay length ξM σ of Mσ as a function of g.For small g, ξM σ is approximately proportional to g −4 .In the crossover between the quadrupolar and confined phases, the behavior of ξM σ changes.The vertical dashed line in panel (b) indicates the value of g corresponding to the peak of the fidelity susceptibility at g ≈ 0.375.

Figure 14 .
Figure 14.Examples of the dynamical charge and electric field distributions in the intermediate region between two opposite static charges.The size of the dots represent the expectation value of the dynamical charge q; the thickness of the links represents the expectation value of the electric field E in arbitrary units.The electric field and charge smaller than 10 −6 are not shown in the figures.Black and blue numbers label examples of the q and E expectation values.Panel (a) represents a typical example in the Coulomb phase (g = 0.01, λ = 0.75); panels (b), (c) and (d) correspond to the values of the coupling constants in Fig. 10 (the stars in Fig. 5) and they represent typical examples taken within the quadrupolar, confined rung-dominated and Higgs regimes respectively.

Figure 15 .
Figure 15.Energy cost ∆E for the introduction of two static charges in the same leg of the ladder as a function of their distance R. The four panels correspond to the same coupling constants chosen in Fig. 14.The grey dashed lines indicate the static charge distance R = 8 used in Fig. 14.

Figure 16 .
Figure 16.Decay of the electric field E between two static charges as a function of the distance r from the negative charge.The distance between the static charges is R = 34.(a) Coulomb phase (log-log scale): E decays approximately as a power law for r < R/2.(b) Quadrupolar phase: E is approximately constant (no screening).(c) and (d) Screened states in the Higgs and confined regimes (log-normal scale).E decays approximately exponentially for r < R/2.The four panels correspond to the same coupling constants chosen in Fig.14.
ACKNOWLEDGMENTS M.B. warmly thanks A. Negretti, E. Rico Ortega and D. Rossini for discussing in depth the physics of LGTs in the ladder geometry and sharing their numerical result for the U(1) truncated gauge symmetry.M.B. is also indebted with E. Cobanera, A. Milsted and G. Ortiz for a useful collaboration about 1D systems with Z N symmetry, and with P. Silvi and H.-H. Tu for insightful correspondence.This work was supported by the Villum Foundation (Research Grant No. 25310) and by the EU's Horizon 2020 programme under the Marie Sk lodowska-Curie grant agreement No. 847523.The DMRG calculations were performed using the ITensor C++ library, https://itensor.org/.

Figure 17 .
Figure 17.Phase diagram from second-order RG for N = 8 (a) and N = 15 (b).The different phases are shown by different colors, including the deconfined phase (purple), the quadrupolar phase (blue), the Coulomb phase (green), the Higgs phase (yellow), the fully confined phase (brown), and the confined rung-dominated phase (red).

2 − 2
dx P cos N θ ρ √ 2 + Cρ cos √ 2N ϕ ρ .(E2)Here Cρ is a suitable renormalized coupling constant determined by the first step in the flow.Its scaling dimension is D Cρ = N/K, such that this term reduces the extension of the gapless phase for N > 4, and it completely removes it for N = 4.In conclusion, the two-step RG procedure indicates how the P and Cρ terms are responsible for gapping the gapless phase for N = 3 and 4. In particular, the gapped phases for N = 2 and 4 coincide, consistently with Appendix B.
Constant Gρ Exp.decay of R and mesons Constant Gρ, Mσ and R Exp. decay of Mρ Or,s = 0 Alg.decay of Os Alg.dec. of Gρ and Mρ Constant Mσ and R Or,s = 0 Exp.decay of Gρ Constant Mρ, Mσ and R