Two-dimensional quantum-link lattice Quantum Electrodynamics at finite density

We present a tree tensor network approach to the study of lattice gauge theories in two spatial dimensions showing how to perform numerical simulations of theories in presence of fermionic matter and four-body magnetic terms, at zero and finite density, with periodic and closed boundary conditions. We exploit the quantum link representation of the gauge fields and demonstrate that a fermionic rishon representation of the quantum links allows us to efficiently handle the fermionic matter while finite densities are naturally enclosed in the tensor network description. We explicit perform calculations for quantum electrodynamics in the spin-one quantum link representation on lattice sizes of up to $16\times16$ sites, detecting and characterizing different quantum phases. In particular, at finite density, we detect a phase separation behaviour as a function of the bare mass values at different filling densities. The presented approach can be extended straightforwardly to three spatial dimensions.

Recent progress in quantum simulations are paving the way to the possibility of studying high-energy physics phenomena with tools developed in low-energy quantum physics [1][2][3][4][5][6][7][8][9][10][11][12][13]. In the Standard Model, forces are mediated through gauge fields, thus gauge invariant field theories -e.g., quantum electrodynamics (QED) for the abelian case or quantum chromodynamics (QCD) for the non-abelian scenario -are fundamental building blocks to our understanding of all microscopic processes ruling the dynamics of elementary particles [14,15]. When discretising the gauge theories, the dynamical gauge variables obey a lattice formulation of the original quantum field theory which is referred to as a Lattice Gauge Theory (LGT) [16,17].
LGTs encode many-body interactions satisfying exact constraints, encoding a latticediscretised version of the local gauge invariance (in QED the Gauss's law ∇·E = 4πρ). Many of the collective phenomena arising from these theories, including the phase diagram, have yet to be fully characterised [18], especially for higher spatial dimensions and at finite charge density.
In this work, we fill this gap and develop a computationally practical Hamiltonian formulation of low-energy QED in two spatial dimensions. We show that TN states allow for an accurate representation of its manybody ground state, thus allowing identifying the different phases of matter, and effectively test the response of the system to a finite density of charge. The study of lattice gauge hamiltonians at finite chemical potential is in general out-of-reach for Montecarlo-based techniques [13,24]: here we show that using an unconstrained Tree Tensor Network (TTN) [72] we can face this highly non-trivial setup. By exploiting the quantum link formalism [2,[51][52][53][54], we put particular emphasis on the natural implementation of the Gauss's law as well as on the emerging additional link symmetry. The quantum-link formalism is of high interest for the quantum science approach to LGTs, as it opens a pathway to develop tensor network and quantum simulations of LGTs [6,10,12,13]. The techniques developed in this paper not only provide the basic ingredients for an efficient calculation of the phase diagram of simple lattice gauge models, but they can be extended to more complex theories and higher dimensions.
We demonstrate the power of this approach by focusing on the low-energy properties, both at zero and finite charge density, of a two-dimensional lattice quantum link theory with U (1) gauge symmetry. Specifically, we investigate a model involving (spinless, flavorless) Kogut-Susskind matter fermions [16,21] and U (1) electromagnetic gauge fields, truncated to a Spin-1 compact representation. We investigate the (zero temperature) phase diagram in the zero global charge scenario: first without magnetic coupling, and afterwards with finite magnetic coupling. We observed that, both magnetic and electric Hamiltonian terms, separately, hinder the creation 2 of a charge-crystal phase, which emerges at large negative bare masses. However, when electric and magnetic terms are mutually frustrated, the charge-crystal phase is restored. Moreover, we study the ground state in the presence of a finite charge density, which we can directly control in the TN ansatz state. Small charge densities impact the zero charge phases in the following way: The vacuum phase allows charges to freely move and reach the system boundaries, manifesting a clear phase separation in space between the bulk and the boundaries; this is reminiscent of the classical electrodynamics properties of a perfect conductor, where ∇ · E = 0 in the bulk and the excess of charge will be redistributed on the outer surface of the conductor, according to a certain surface charge density profile. By contrast, the charge-crystal phase, which is full of matter/antimatter, is characterised by a homogeneous delocalisation of the charge-hole, resulting into a quasi-flat charge distribution in the bulk, therefore reminiscent of a plasma phase [55].
These results further motivate the development and refinement of TN techniques which turn out to be very well suited to study highly nontrivial two-dimensional models, where local gauge constraints have to be properly encoded. Finally, we stress that the quantum link formulation provides the ideal tools to establish a connection between abelian or non-abelian LGTs and atomic lattice experiments [56,57]. In this framework, the dynamical gauge fields are usually represented by spin degrees of freedom, which have a natural mapping to typical condensed-matter models, like Hubbard Hamiltonians or locally constrained Ising-like Hamiltonians. These models can be engineered with cold atoms in optical lattices [58,59], or within the very promising experimental setups involving Rydberg atom chains [60,61].
The paper is structured as follows: In Sec. I, we present the 2D lattice gauge Hamiltonian and its quantum link formulation in terms of the gauge field Spin-1 compact representation. We also give some technical details of the Tensor Network numerical simulations. In Sec. II we focus on the ground-state properties in the zero-charge sector: we explore the phase space of the model by varying the mass and the electric coupling; we then analyse the effect of a finite magnetic coupling. Sec. III is devoted to study the equilibrium properties at finite charge density. We exploit TTN techniques to investigate how the charges redistribute all over the lattice depending on the phase of matter the system belongs to. Finally, we draw our conclusions in Sec. IV and give additional supplementary technical details in the Appendices.

I. MODEL AND METHODS
We consider a field theory on a 2D square lattice with U (1) local gauge symmetry. The sites of a finite L×L square lattice host the matter field, while the quantum gauge field lives on the lattice links, with open boundary conditions. Following the Kogut-Susskind (staggered) Sketch of the 2D lattice gauge theory in the Spin-1 representation. The cartoon of the square lattice shows a specific gauge-invariant configuration of the matter and gauge fields with zero total charge (3 particles and 3 antiparticles). Staggered fermions represent matter and antimatter fields on a lattice bipartition: on the even (odd) bipartition, a full red (blue) site represents a particle (antiparticle). The gauge field points in the positive direction of the color gradient.
formulation [16,21], the discretisation of the matter field is performed by introducing a staggered fermionic field, whose positive energy solutions lie on the even sites and negative ones on the odd sites. The matter field is thus described by spinless, flavorless Dirac fermions, whose operator algebra satisfies the usual canonic anticommutation relations {ψ x ,ψ † x } = δ x,x . In particular, in the even sub-lattice particles represent fermions with electric charge +q ('positrons'), whilst in the odd sublattice holes represent anti-fermions with electric charge −q ('electrons'). Here a lattice site x labels a 2D coordinate x ≡ (i, j), and the parity p x ≡ (−1) x = (−1) i+j of a site, distinguishing the two sub-lattices, is well defined on the square lattice (see Fig. 1).
The gauge field is defined on the lattice links, its algebra constructed by the electric field operatorÊ x,µ = E † x,µ and its associated parallel transporterÛ x,µ , which is unitary U x,µ U † x,µ = 1 and satisfies [Ê x,µ ,Û y,ν ] = δ x,y δ µ,νÛx,µ . Here µ (ν) represents the positive unit lattice vector in one of the two orthogonal directions, namely µ x ≡ (1, 0) and µ y ≡ (0, 1), thus (x, µ) uniquely defines a link. For comfort of notation, we also allow (technically redundant) negative unit lattice vectors −µ x and −µ y , with the conventionÊ x+µ,−µ = −Ê x,µ , and in turnÛ x+µ,−µ =Û † x,µ . With such definitions, apart from a rescaling due to the lattice spacing regularisation [16,21], the two-dimensional lattice QED Hamilto-nian, including a magnetic plaquette term, readŝ where the coordinate µ runs in {µ x , µ y }. The first term in Eq. (1) provides the minimal coupling between gauge and matter fields associated with the coupling strength t. It describes a process of fermion-antifermion pair creation/annihilation, where the parallel transporter operator guarantees that the local gauge symmetries are not violated. The second term in the Hamiltonian represents the energy associated to the fermionic bare mass, and it appears as a staggered chemical potential according to the Kogut-Susskind prescription. For numerical purpose, it has been redefined by adding an overall constant mL 2 /2, thus replacing . This way, a filled local state in the even sub-lattice cost positive energy m and carries charge q; otherwise, when an odd site is empty, the energy cost is still m, but it corresponds to having an antiparticle (a hole) with charge −q. The last two terms contribute to the gauge field dynamics: the electric part with coupling g e , is completely local. The magnetic part, with coupling g m instead, is constructed by considering the smallest Wilson loop -product of parallel transporter U x,µ in a closed loop -the size of a plaquette. Its name is related to the fact that it generates the magnetic contribution to the energy density in the continuum limit.
The LGT HamiltonianĤ commutes with the local Gauss's law generators (in unit of q) where the unit lattice vector µ in the sum runs in {±µ x , ±µ y }, while p x = (−1) x is, again, the lattice site parity. In addition, the model exhibits an U (1) global symmetry, namely the conservation of the total chargeQ = x [ψ † xψx − 1−px 2 ] = − L 2 2 +N , equivalent, apart from a constant, to the number conservation N = xψ † xψx of Kogut-Susskind matter fermions. As a consequence of the convention, usingÊ x,−µ = −Ê x−µ,µ , the sum of all 4 terms of the gauge field around the lattice site x corresponds to the outgoing electric flux, i.e. µÊ x,µ = E x,µx + E x,µy − E x−µx,µx − E x−µy,µy . The gauge invariant Hilbert space is thus given by all states |Φ satisfyingĜ x |Φ = 0 at every site x. As each electric field degree of freedom is shared by two Gauss' generators G x , the generators themselves overlap, and projecting onto the gauge-invariant subspace becomes a nonlocal operation. Only for 1D lattice QED, or lattice Schwinger model [19], it is possible to integrate out the gauge variables and work with the matter field only (albeit with long-range interactions) [62]. However, in two dimensions, a given (integer occupation) realisation of the matter fermions does not fix a unique gauge field configuration, thus requiring explicit treatment of the gauge fields as quantum variables. A numerically-relevant complication, related to the standard Wilson formulation of lattice gauge theories, arises from the gauge field algebra, [Ê,Û ] =Û withÊ =Ê † andÛÛ † =Û †Û = 1, whose representations are always infinite dimensional. Simply put, if a representation contains the gauge field state |α , such thatÊ|α = α|α with α ∈ R, then the states |α±1 =Û ±1 |α belong to the representation as well. By induction, the representation must contain all the states |α + N , which are mutually orthogonal as distinct eigenstates ofÊ, thus the representation space dimension is at least countable infinite.
In order to make the Hamiltonian numerically tractable via Tensor Network methods, we need to truncate the local gauge field space to a finite dimension. For bosonic models, this is typically done by introducing an energy cutoff and eliminating states with singlebody energy density beyond it, while a posteriori checking the introduced approximation. Similarly, for U(1) lattice gauge theories, we truncate the electric field according to the quantum link model formulation. Specifically, the gauge fields are substituted by Spin operators, namelyÊ x,µ = (Ŝ z x,µ + α) andÛ x,µ =Ŝ + x,µ /s, such thatÊ is still hermitian and the commutation relation [Ê x,µ ,Û y,ν ] = δ x,y δ µ,νÛx,µ is preserved [2], howeverÛ is no longer unitary for any finite spin-s representation |Ŝ| 2 = s(s + 1)1. The original algebra is then restored in the large spin limit s → ∞, for any background field α ∈ R. Similar truncation strategies, based on group representations, can be applied to non-Abelian gauge theories as well [35,63]. In the following, we make use of the Spin-1 representation (s = 1), under zero background field α = 0, which captures reasonably well the low-energy physics of the theory, especially in the parameter regions wherein the ground-state is characterised by small fluctuations above the bare vacuum. s = 1 is the smallest spin representation exhibiting a nontrivial electric energy contribution. In fact, for s = 1/2, we have thatÊ 2 x,µ ∝ (σ z x,µ ) 2 = 1 is simply a constant in the Hamiltonian, thus g 2 e plays no role. Finally, in 1D it was observed that truncated gauge representations converge rapidly to the continuum theory, e.g. in the Schwinger model [37,64,65], reinforcing quantitative validity of the results obtained in the simplified model.
Let us mention that, to highlight the connection to the Wilson approach in the Hamiltonian formulation of QED in 2D, the electric and magnetic couplings should be related between each other such that g 2 e ∼ g 2 a −2 and g 2 m ∼ g −2 a −2 , where g is the electrodynamic coupling and a the lattice spacing. Similarly, mass and kinetic coupling should behave respectively as m ∼ 1 and t ∼ a −1 . To recover the continuum theory, matter and gauge fields need to be proportional to the lattice spacing [16,20]. However, in the following we treat each parameter in the LGT Hamiltonian (1) independently: indeed, we consider the model interesting per se, and leaving more detailed convergence analysis along the lines of those already presented in one-dimensional systems to future work [13].
In this work we introduce an algebraic technique, similar to the rishon representations common in quantum link models [2,[51][52][53][54], which has the advantage of automatically accounting for the Gauss' law, while carefully reproducing the anticommutation relations of the matter fermions without resorting to Jordan-Wigner string terms (see next section). This strategy relies on splitting the gauge field space on each link (x, µ) into a pair of 3-hardcore fermionic modes, defined later on. We say that each mode in this pair 'belongs' to either of the sites sharing the link, in this case x and x + µ.
Thus, we writeÛ x,µ =η x,µη † x+µ,−µ , where the 3hardcore fermionic operatorsη satisfyη 3 x,µ = 0 (whilê η 2 x,µ = 0) and anticommute at different positions {η x,µ ,η ( †) y,ν } = 0, for x = y or µ = ν. Moreover, these new modes obey anti-commutation relations with the matter field as well, i.e. {η ( †) x,µ ,ψ ( †) y } = 0. To explicitly build this 3-hardcore fermionic modeη x,µ , we use two sub-species of Dirac fermionsâ x,µ andb x,µ , such that wheren a x,µ =â † x,µâx,µ andn b x,µ =b † x,µbx,µ are the occupation number operators for each sub-species. This construction provides the local algebra and grants only access to the 3-dimensional subspaces for each 3-hardcore fermion mode, spanned by the following three states where |0 x,µ is the Dirac vacuum of both sub-species, i.e. a x,µ |0 x,µ = b x,µ |0 x,µ = 0. The stateb † x,µ |0 x,µ is disconnected from the other three and thus projected away. Such "half-link" local subspace is joined with a similar constructionâ † x+µ,−µ andb † x+µ,−µ on the other half of the link, thus the pair defines the link space, and E x,µ will be diagonal in the occupation basis. While, in principle, a full link space is 3×3=9-dimensional, we can now exploit the following symmetrŷ which counts the total number of fermions in each link and is a conserved quantity since [L x,µ ,Ê x,µ ] = [L x,µ ,Û x,µ ] = 0. By working in the subspace with two fermions per link,L x,µ = 2, we reduce the link space to dimension 3 and we can restore the desired algebra. First of all, we write the occupation basis (see also Fig. 2) so thatÛ x,µ acts correctly, i.e. where the minus sign in the first equation ensuresÛ x,µ |ø = | → and U † x,µ | → = |ø . Then, we express the electric field operator as the unbalance of fermions between the two halves of a link, preciselŷ implementing the correct action of E x,µ . It is worth to mention that this formulation can be extended to higher Spin-j representations, where each link becomes a pair of (2j + 1)-hardcore fermions.

B. The local gauge-invariant dressed sites
One of the common issues while working with a lattice gauge theory, even in the compact representation of the electric field, is to properly identify the gauge-invariant Hilbert space. Due to the overlapping of the Gauss' law generatorsĜ x , the identification of the correct local basis is highly non-trivial, especially for dimensions higher than one [26].
Using the 3-hardcore fermion pairs language gives us a shortcut to this issue. In fact, we are able to recast the Gauss' law generators as non-overlapping operators, at the price of enforcing the link constraint (L x,µ − 2)|Φ = 0. Using this constraint, we can rewrite the electric field operator in Eq. (8) taking only the fermionic operators into account, which act on the "halflink" connected to x, i.e.
valid in the link-symmetry invariant space (L x,µ − 2)|Φ = 0. As a consequence, in the Hilbert

Spin-1 Quantum Link Formulation
Local gauge-invariant dressed site Mapping of the gauge field states in the Spin-1 representation to the fermionic Fock states. Each half link is constructed by employing two species of Dirac fermions. The link symmetry formally reduces the total number of states to only the 3 allowed states with 2 fermions. We therefore construct the local gauge-invariant dressed site by gluing each single matter site together with its neighbouring half links. In the four examples, notice how the quantum number φ = 1 represents the presence (absence) of a charge (anti-charge) in the even (odd) sites.
space with 2 Dirac fermions per link, the Gauss's law generators become strictly local, i.e. containing quantum variables belonging solely to site x, and read Within this picture, it is easy to identify a local gaugeinvariant basis for the dressed site × |φ x |k 1 x,−µx |k 2 x,−µy |k 3 x,µx |k 4 x,µy , 1} describes the matter content, while k j ∈ {0, 1, 2} selects a state, from those in Eq. (5), for each respective half-link. The factor (−1) δ k 1 ,2 +δ k 2 ,2 accomodates the sign in Eq. (7). In this language, the Gauss' law, cast as Eq. (10), simplifies to which fixes the total number of fermions in each dressed site, specifically 4 in the even sites, and 5 in the odd sites. Eq. (12) actually reduces the Hilbert space dimension of each dressed site from 162 to 35, and we use these 35 states as computational basis for tensor network algorithms.
A fundamental feature of this language is that, since the total number of fermions at each dressed site is conserved, their parity is conserved as well, thus the gaugeinvariant model will exhibit no Jordan-Wigner strings (outside the dressed sites) in the computational basis. An operative way to show this property is to consider that the Hamiltonian termψ † xÛx,µψx+µ , decomposes as the product ofψ † xηx,µ andη † x+µ,−µψx+µ : each of these two factors is local (acts on a single dressed site) and commutes with the algebra of other dressed sites. The same applies to the magnetic plaquette term. In conclusion, by working on the dressed-site computational basis, we can employ standard (spin-model-like) tensor network techniques, without the requirement of keeping track of fermionic parity at each site [66][67][68][69][70][71]. Notice that this construction can be exploited also to perform quantum computations of two-dimensional LGTs.

C. Tensor Network for 2D lattice gauge simulations
In order to numerically simulate the quantum system, we use a two-dimensional Tree Tensor Network (TTN) state to represent the many-body wave function [72][73][74]. We work in the computational 35-dimensional local basis for each dressed site, defined in the previous section, which automatically encodes the Gauss' law. Operators appearing in the Hamiltonian (1) can be cast in this basis, either as local operators, or acting on a pair or plaquette of neighbouring dressed site (see Appendix A for the explicit construction). The extra link symmetryL x,µ = 2 must be enforced at every pair of neighboring sites. We do so by introducing an energy penalty for all states violating the link constraints. This penalty term is included in the optimisation by a driven penalty method -similar to an augmented Lagrangian method -which is described in more detail in Appendix B. Under all other aspects, the TTN algorithm employed here for finding the manybody ground state follows the prescriptions of Ref. [75].
In the numerical simulations we fixed the energy scale by setting the coupling strength t = −1. Furthermore, we worked within a sector with a fixed total chargeQ, by using standard techniques for global symmetry conservation in TNs [75][76][77]. We thus characterised the groundstate properties as function of the mass m, the electric coupling g e and the magnetic coupling g m .
In order to exploit the best performances of our TTN algorithm, we ran simulations on square lattices L × L with the linear length L being a binary power; in particular, we considered L = 4, 8 and 16, and varied the TN auxiliary dimension (or bond dimension) up to χ ∼ 300. Depending on L and the physical parameters, we obtain a convergence precision between ∼ 10 −2 −10 −5 , sufficient to characterise the phase properties.

II. ZERO CHARGE DENSITY SECTOR
In this section, we focus on the zero charge density sector ρ ≡ Q /L 2 = 0, where there is a balance between matter and antimatter, and analyze the ground state of Hamiltonian (1) within this subspace. Unless otherwise stated, we consider periodic boundary conditions. We characterise the ground state of the Hamiltonian by looking at the energy density Ĥ /L 2 , and the particle density n = 1 x ) counts how many charges are in the system, both positive and negative, i.e. fermions in even sites plus holes in odd sites. We start our analysis by first focusing on the case in which the magnetic coupling has been set to zero, g m = 0. Before detailing the numerical results, some analytically-solvable limit cases should be considered. For large positive values of the bare mass m t, the fluctuations above the bare vacuum are highly suppressed; the system exhibits a unique phase since there is no competition between the matter term and the electric field term in the Hamiltonian. Indeed, to construct pairs of particle/antiparticle, the matter energy and the electric field energy both contribute to an overall increasing of the ground-state energy. In order to explore more interesting phenomena, we allow the mass coupling to reach negative values. Doing so, we can identify two different regions depending on the competition between the electric coupling g 2 e /2 and the values of the mass m < 0: (i) for g 2 e /2 2|m|, we still have a vacuum-like phase, where we expect a unique non-degenerate ground-state with small particle-density fluctuations. This phase exists, no matter the value of the mass, as far as the energy cost to turn on a non-vanishing electric field on a single link overcomes the gain in creating the associated pairs of particle/antiparticle. Indeed, for any value of the mass and g 2 e /2 → ∞, or for g 2 e /2 = 0 and m → ∞, the pres-ence of a finite electric field, or finite particle density, is strictly forbidden and the ground-state flows toward the only admissible configuration, namely the bare vacuum.
(ii) for −2m g 2 e /2 > 0 the phase of matter is characterised by slightly deformed particle-antiparticle dimers; this phase of course only exists for negative value of the mass and represents the region wherein the energy gain for creating a couple of particle/antiparticle largely overcomes the associated electric field energy cost.
Here the ground-state remains highly degenerate as far as the kinetic energy coupling |t| is much smaller than all the others energy scales (degeneracy being lifted only at the fourth order in t). In particular, for g 2 e /2 = 0 and m → −∞ the ground state reduces to a completely filled state. In order to minimise the electric field energy, particles and antiparticles are arranged in L 2 /2 pairs (where we are assuming L even) sharing a single electric flux in between. All these configurations are energetically equivalent and their degeneracy corresponds to the number of ways in which a finite quadratic lattice (with open or periodic boundary conditions) can be fully covered with given numbers of "horizontal" and "vertical" dimers. This number scales exponentially with the system size as exp(L 2 C/π) for L → ∞, with C 0.915966 the Catalan's constant [78]. For sake of clarity, we stress that such 'dimers' are not entangled clusters of matter and gauge fields; they are roughly product states.
Let us mention that the case g e = 0 with m → ∞ (m → −∞) is more pathological since any gauge-field configuration compatible with the vacuum (dimerised) state is admissible, provided the Gauss's law is fulfilled. In practice, we may draw a generic closed loop with finite electric flux on top of the vacuum state without modifying its energy; similar gauge loops may be realised on top of the dimerised state, provided it is compatible with the occupied links, without changing its energy as well. All these configurations are gauge-invariant by construction, and increase the degeneracy of the ground-state energy sector.
Our numerical results confirm and extend this picture, as it can easily be seen in the phase diagram displayed in Fig. 3, obtained from TTN simulations in a 8 × 8 system. The matter density is roughly zero in the vacuum phase; otherwise, it takes on a finite value whenever the system exhibits "dimerisation", i.e. in the charge-crystal phase. We checked that the numerical data, both the groundstate energy density and the particle density, show an asymptotic tendency toward the perturbative estimates. Interestingly, the particle density experiences an abrupt change mainly in a narrowed region around m −g 2 e /4, where the local slope is becoming steeper as the electric coupling (and the mass) is approaching zero (see left panel in Fig. 4), as roughly predicted by perturbation theory and supported by the exact results in the 2 × 2 case (see Appendices C and D).
As a confirmation of this scenario, we expect particle fluctuations to be enhanced around such region, mainly due to the very low energy cost in creating particle- antiparticle pairs. This feature can be highlighted by tracking the density-density correlation functions, precisely C x,x = Ô xÔx , and its spatial averagesC(v) = L −2 x C x,x+v , using a local order parameterÔ x = (−1) x (2ψ † x ψ x − 1) which does not discriminate between the two phases. By evaluating the corresponding correlation length x + v 2 y , we observe a sharp drop in proximity of the critical point, signaling competition between density-ordered phases and dominance of fluctuations. Such behavior is shown in the right panel of Fig. 4, where various sizes are compared, up to 16 × 16 matter sites. The right panel clearly shows that the actual transition point is slightly below the classical (i.e., t = 0) position g 2 e /2 = −2m since particle/antiparticle fluctuations, induced by a finite value of the hopping amplitude t, naturally discourage the charge-crystal order. This effect emerges already at the second order in perturbation theory treatment (see Appendix C), where the crossing point of the two different ground-state energies, E v (vacuum) and E d (dimer), slightly shifts toward the dimerised phase.
A relevant physical question is whether the system undergoes an actual quantum phase transition across the two regions. Exactly at t = 0, when m crosses the critical value −g 2 e /4, the ground-state exhibits an exact levelcrossing, passing from the bare vacuum to the chargecrystal energy sector. In this limit case, the system experiences a trivial first-order phase transition, since the gauge-field energy term and the matter-field mass term commute between each other. However, if we tune the mass at the classical critical value m = −g 2 e /4, a small hopping amplitude t = 0 is already sufficient to remove such degeneracy: namely, the bare-vacuum energy and the charge-crystal energy get modified in a different way so that a gap opens between the two sectors. At the critical value of the mass, creation/annihilation of pairs of particle/antiparticle has no energetic cost, and the ground state energy sector is characterised by all possible states with any number of dimers; however, creating a pair in the vicinity of the bare vacuum is more favourable than annihilating a pair on top a fully dimerised state; at least, this is true for any finite L (see Appendix C for details).
A crucial insight comes from the features of the overlap between the exact ground-state |GS and the unique bare vacuum |Ω (see Appendix D for exact results in the 2×2 case). Indeed, for the t = 0 trivial case, it experiences a discontinuous transition when passing from the vacuum sector to the full dimerised sector, suddenly jumping from one to zero. Interestingly, for fixed system size L, we may evaluate such overlap in the approximate groundstate |GS (k) at a given order k in perturbation theory. The resulting perturbative expansion of the square of the overlap | Ω|GS (k) | 2 changes continuously in the vacuum phase, while it remains identically vanishing, i.e. | Ω|GS (k) | 2 = 0, when correcting the fully dimerised state up to k < L 2 /2. We thus expect the exact overlap to be continuous at the transition point and identically vanishing in the thermodynamic limit L → ∞. As a consequence, for any finite t, no first-order phase-transition occurs, and we may have a second-order phase transition. Let us stress that, although the perturbative expansion of the fully dimerised state does not produce any change in the overlap with the bare vacuum for k < L 2 /2, local observables do experience perturbative modifications, simply because the state by itself gets modified. In particular, as a consequence of the Hellmann-Feynman theorem, the particle density as a function of the mass coupling m coincide with the derivative of the ground-state energy density, GS(m)|n|GS(m) = ∂ m E GS (m)/L 2 . A second-order phase transition will thus imply continuous profiles of the particle density, with a discontinuous or diverging derivative at the transition point. In fact, we have numerical evidence that the matter density changes continuously when going from one phase to the other (see Fig. 4), however, it remains very hard to infer about its derivative at the transition point.

A. Finite magnetic-coupling effects
We now analyze the case of non-vanishing magnetic coupling g m , especially focusing on how it impacts the many-body quantum phases at zero temperature.
In Figure 5 we show the field-plot representations of the ground-state typical configurations for an 8 × 8 system in the presence of magnetic couplings. For the sake of visibility we only plot a 4 × 4 subsystem out of the complete 8 × 8 lattice simulated with periodic boundary conditions. Both mass m and electric coupling g e have been chosen so that the system is well deep within the two different phases (left panels: vacuum, right panels: charge crystal). As the magnetic coupling g m is increased to commensurate values (bottom panels) we see negligible changes affecting the vacuum phase. By contrast, in the charge-crystal phase, the non-vanishing magnetic coupling introduces a nontrivial reorganisation of the electric fields.
Such an effect can be well understood in terms of perturbation theory: (i) in the vacuum phase, the groundstate is not degenerate and the first nontrivial corrections are given by coupling such state with all the states with a single flux loop over a single plaquette (whose energy is therefore 2g 2 e ). In this regime, the flux loop state has high electric field energy, thus it will impact only slightly the global features of the state. The first order correction to the ground-state energy will be quadratic in the magnetic coupling, i.e. ∼ g 4 m /g 2 e (see Fig. 6 left panel). Let us stress that, even though the state may experience an electric-field reconfiguration due to the "field-loop" superpositions, the fact that in this phase the electric field is almost zero, causes no visible effect in its expectation value; this is pretty clear from the left column of Fig. 5: when passing from a small (g 2 m /2 = 0.2) to a slightly bigger value (g 2 m /2 = 2) of the magnetic coupling, the changes in the expectation value of the gauge field are negligible. (ii) in the charge-crystal phase, the effect of the magnetic interaction is nontrivial. Indeed, the ground-state energy sector in this phase is highly degenerate, and the magnetic field contributes to lift such degeneracy (at much lower order than t). The magnetic coupling introduces first order transitions between different gauge field configurations therefore its first contribution to the ground-state energy is order ∼ g 2 m /g 2 e . Actually, a sufficiently large value of the magnetic coupling g m helps the TTN wave function to restore the square lattice symmetry by introducing a gap between the actual ground-state and all the energetically unfavourable configurations. This property is noticeable in Fig. 5 (bottomright panel) where for g 2 m /2 = 2 the gauge field distribution becomes uniform (on average) in the bulk. In this scenario, the charge crystal does not encourage the formation of dimers, but instead a global entangled state of gauge fields.
The previous considerations are supported by the behaviour of the ground-state energy and of the particle density, as a function of g m . In Fig. 6 we plot the numerical results obtained via TTN simulations in a 8 × 8 system. We fixed the value of the mass to m = −2, and explored the behaviour for g 2 e /2 = 2 and 6, that are slightly below or above the classical transition point g 2 e /2 = −2m. We vary the magnetic coupling in a rather big interval g 2 m /2 ∈ [0, 8]. In the first case, i.e. when the system is initially in the charge crystal phase, we expect linear corrections to the ground state energy as a Numerical results for the ground state energy Ĥ /L 2 density and the particle density n as a function of the magnetic coupling for a 8 × 8 systems. The mass coupling has been fixed to m = −2. The energies have been renormalised to their absolute value at gm = 0; dashed lines are guide for the eyes. Data have been obtained by extrapolating from TTN simulations with different auxiliary dimensions.
function of g 2 m ; this is pretty clear from Fig. 6 where, however, some deviations are visible due to the vicinity of the phase boundary. In the second case, i.e. starting from the vacuum phase, a quadratic deviation of the ground state is clearly visible (see Fig. 6).
Interestingly, in the parameter region we are exploring, the magnetic coupling enhances the production of particles, thus increasing the average matter density, even though the magnetic term does not directly couple to matter. In practice, the magnetic coupling creates resonating configurations of the gauge fields in the crystal charge phase, thus decreasing the electromagnetic energy density of the state itself, which in turn favours the crystal charge phase in proximity of the phase boundary. Hence small g m values effectively enlarge the charge crystal phase. However, in the Spin-1 representation of the gauge field, the crystal charge phase is not stable under an arbitrary large value of the magnetic coupling, and we expect n x = 0 when g m g e . This can be easily understood at the classical level (t = 0) comparing the effect of the a Wilson loop operator in the zeromatter (vacuum) sector and in the full-matter sector: in the former case, each single plaquette is resonating between three different diagonal gauge-field configurations {| , |ø , | }, in the last case, only two configurations are resonating, e.g. {| ↑↓ , | }, since constructing a clockwise(anticlockwise) electric loop ( ) on top of the first(second) state is forbidden by the Spin-1 finite representation. This originates an energy gain which is proportional to − √ 2g 2 m for a plaquette in the vacuum, whilst it is only −g 2 m for a dimerised plaquette. In practice such difference remains at the many-body level as well, and, therefore, although for −2m g 2 e /2 the dimerised configuration represents the lower energy state at g m = 0, it will be getting energetically unfavourite for sufficiently strong magnetic couplings.
We have further analytical confirmation of such behaviour from the exact diagonalisation of 2 × 2 systems (see Appendix D for details): for a single plaquette system, the first visible effect of a non-vanishing magnetic coupling is to mix up the two dimerised states into two different superpositions with different energies. The transition between the vacuum state toward the lower energetic charge-crystal state is therefore sharpened and its position is shifted as well in g 2 e /4 + m (g 2 e + g 2 m /2 − g 4 e + g 4 m /2)/4. Interestingly, depending on the values of g e , this shifting is not monotonous in g m , producing an initial increase in the particle density followed by a definitive decrease toward zero (c.f. Fig. 16), and thus confirming the previous heuristic argument based on perturbation theory.

III. FINITE CHARGE DENSITY SECTOR
One of the most intriguing phenomena we observed in our numerical simulations rely on the possibility to create a charge imbalance into the system. This scenario is challenging for Montecarlo techniques as it produces the sign problem [13,24]. Instead, our gauge invariant tensor network approach is very well suited to overcome such difficulty: the fact that the global U (1) symmetry has been explicitly embedded into the Tensor Network ansatz [75], allows to work exactly within each sector with fixed total charge. In the following we only considered g m = 0. Moreover, in this setup, due to the finite net electric flux coming out from the entire system, we have to work with open boundary conditions. In this geometry, the dressed sites at the boundary are now characterised by one outgoing half-link (two in the corners) In the vacuum phase, the excess of charge prefers to be localised at the boundaries, since such configurations are more energetically favourable. In the dimerised phase, the holes may occupy any position since the system can reconfigure the pairs of dimers in a way to pay always the same amount of energy. However, due to the very high degeneracy of the low-energy sector, the TTN simulations may get stuck into a slightly asymmetric configuration.
which can support electric field to allow the existence of a non-vanishing total outgoing flux.
When a finite density of charge ρ ≡ Q /L 2 ∈ {−1/2, 1/2} is injected into the system, we expect a different behaviour depending on the part of the phase diagram the ground state is belonging within. Indeed, when the ground state is very close to the bare vacuum, any charge created on top of it is forced to reach the boundaries so as to minimise the total energy; this is easily understood already with the classical (t = 0) Hamiltonian, and there are no fluctuations of the gauge fields. In this case, a classical configuration with a single charge located at distance from the boundary costs at least g 2 /2 more than the optimal configuration where the same charge is located at the surface (see Fig. 7). In this phase the diagonal energy term gets modified as E v /L 2 = (g 2 e /4 + m)ρ, as far as Q ≤ 2(L − 1), i.e. whenever the total excess of charge is lower than the number of allowed free sites at the boundaries. When the total charge gets larger, deeper sites start to be filled, e.g. for 2(L − 1) < Q ≤ 4(L − 2) one starts filling the next-neighbouring sites to the surface (e.g. Fig. 8).
In general thus, we have a phase separation between a boundary region attached to the surface, or strip, where all charges are localised, and the bulk region where there is neither charge nor electric field. In practice, defining ρ ≡ 2 (L− )/L 2 being the maximum amount of chargedensity the system can store within a strip of extension from the surface, we have a sharp discontinuity in local charge densities, at the smallest * such that ρ ≤ ρ * , between a finite-charge region (for j < * ) and a zerocharge region (for j > * ). In particular, in the thermodynamic limit we obtain * /L = (1 − 1 − 2|ρ|)/2. We stress that both * and L − * are extensive quantities, thus the phase separation argument remains valid in the thermodynamic limit: in practice, as far as the average charge density is such that |ρ| < 1/2, we always have an extensive region in the bulk of the system, whose linear dimension scales as L − * , which exhibits no charges.
We expect this picture to be slightly modified at finite hopping coupling |t|, but to remain valid as far as the system belongs to the vacuum phase. In practice, a fi- . Notice how the excess of charge is larger than the allowed antimatter sites at the surface: the system has to allocate two extra charges in the next-neighbouring sites to the surface.
nite tunnelling amplitude introduce a small homogeneous particle density, thus slightly rising * and having the effect as well to build up a finite charge-penetration length ∼ |t| such that the transition at ∼ * becomes smooth, with an exponentially small density-charge tail penetrating into the bulk. The overall scenario is confirmed by the field plots in Fig. 7 and Fig. 8, where it is pretty clear that when the couplings are tuned in order for the system to be deep into the vacuum phase, the excess of charges stick to the boundary so as to minimise the length of the attached electric strings. In principle, all possible configurations with all charges at the boundaries are energetically equivalent. Anyhow, the TTN many-body wave function spontaneously breaks such symmetry and picks up a single specific configuration, as it is usually the case with DMRG-like alghorithms.
When the state belongs to the charge crystal phase, a finite positive (negative) charge density is mainly generated by creating holes in the odd (even) sub-lattice; namely, negative (positive) charges are removed from the fully dimerised state. In order to minimise the energy, the holes can be now fully delocalised: a hole in the bulk, or at the boundary, generates a reconfiguration of the charge crystal state in such a way to always pay the same amount of energy, and guarantee the expected total outgoing electric flux. In this phase, the zero-order energy term gets modified as E d /L 2 = (g 2 e /4 + m)(1 − ρ). The entire system is now characterised by a unique spatial phase where we expect a uniform average charge density and finite electric field in the bulk. Let us mention that, for any finite value of the hopping amplitude, we still expect a similar behaviour, where the transition toward the phase-separated phase will be driven by the competition between the mass and the electric coupling.
In order to highlight the different features of the lowenergy state at finite chemical potential, we analysed the behaviour of the surface charge density where D is a square which counts dim D = 4(L + 1 − 2 ) lattice sites as sketched in Fig. 9. Here ∈ {1, 2, . . . , L/2} represents the distance of the domain D from the external surface: namely, as grows, we select domains deeper into the bulk.
In Fig. 9 we plot the surface charge density σ as function of for different point in the coupling-parameter space. As far as the Hamiltonian is tuned into the vacuum phase, the surface charge suddenly drops when getting into the bulk of the system. As expected, for finite value of the couplings, when approaching the critical region, the bulk charge density gets enhanced; finally, once the system reaches the charge-crystal phase, σ acquires a loosely uniform shape.
Finally, we carefully checked the ground-state energy density and the particle density, which are plotted in Fig. 10 as a function of the mass for two different values of the electric coupling. Notice how, for sufficiently large positive (negative) value of the mass the data get closer to the perturbative predictions. The intermediate region, at m ∼ g 2 e /4, is characterised by stronger quantum fluctuations and thus exhibit a smooth transition between uniform and nonuniform charge distribution in space.

IV. CONCLUSIONS
In this work we demonstrated a novel tensor network approach to the study of two-dimensional lattice gauge theories. By exploiting the quantum link formulation of LGT, the fermionic rishon representation of quantum links, and unconstrained Tree Tensor Networks, we investigated the equilibrium properties of a two-dimensional lattice QED within its first compact spin representation. We present results for lattice size up to 16x16, whose Hilbert space dimension is approximatively equivalent to that of a system composed by spins one-half on a square lattices with edges of eighty lattice sites. Whenever possible, we confirmed our results with perturbative analysis and small scale exact simulations.
In particular, we identified different phases at zero chemical potential, a vacuum state and a charge-density one, that reproduce what has been found in the onedimensional case, and investigated the effects of a magnetic term uniquely present in two-dimensions. Finally, we explore the finite density scenario and individuate two distinct behaviours correspondent to the vacuum and charge-density phases: in the former case, the excess charges accumulate on the boundaries to minimize the electric fields to be energetically sustained, as for classical charged conductors. In the latter, the excess charge is distributed uniformly in the bulk and boundaries.
In conclusions, we have shown that unconstrained Tree Tensor Network turns out to be a powerful tool to obtain a non-perturbative description of a lattice gauge theory in two dimensions. We stress that these simulations have been obtained on standard clusters without exploiting heavy parallelization and with simulations lasting a few days. Despite the fact that the presented results are not yet at the level to allow physical predictions in the continuous limit for the system we study, we foresee that upgrading the current software to exploit the full power of High Performance Computing -without mayor changes in the algorithms -larger system sizes, an additional dimension, the continuum and Wilson limit, and more complex theories will be in range of the approach presented here as already shown for the one-dimensional case [31,62].

V. ACKNOWLEDGMENTS
We are very grateful to M. Dalmonte, K. Jansen, L. Salasnich, U.J. Wiese and P. Zoller for valuable comments. Authors kindly acknowledge support from the BMBF and EU-Quantera via QTFLAG, the Quantum Flagship via PASQuanS, the INFN, the MIUR via the Italian PRIN2017, the DFG via the TWITTER project, the US Air Force Office of Scientific Research (AFOSR) via IOE grant FA9550-19-1-7044 LASCEM, and the Austrian Research Promotion Agency (FFG) via QFTE project AutomatiQ.

Appendix A: Constructing the computational Hamiltonian
In this section we sketch the steps needed to obtain the operator matrices, and their elements, which appear in the computational formulation of the quantum link QED model. In particular, we stress how to construct buildingblock operators A (α) j , each acting on a single dressed site j, which are genuinely local, in the sense that they commute, by construction, with every other building-block operator at another site: j =j ] = 0. The electric field term and the bare mass term are diagonal in the 13 occupation basis of fermions and rishons, as Eq. (11), and thus trivially obtained. The non-diagonal terms are decomposed as follows: Matter-Field coupling terms − Matter-field terms decompose naturally as ψ † x U x,x+µx ψ x+µx = A (and its hermitian conjugate) for horizontal 'hopping' terms, and ψ † x U x,x+µy ψ x+µy = A x+µy for vertical hopping. The decomposition into building blocks is based upon Where η x,µ are the 3-hardcore fermionic operators defined in Eq. (3). Both A (1) x and A x are built on an even number of fermionic operators, thus they commute with any operator which does not act on site x, thus genuinely local. The vertical hopping term is similarly decomposed into building-block operators.
Magnetic terms − The magnetic (or plaquette) term decomposes into building block operators, acting on the four dressed sites at the corners of a plaquette. Specifically, we have In this section, we describe the Tensor Network (TN) approach for LGT in more technical details. As mentioned in section I C we already fulfill the Gauss law by choosing the local gauge-invariant states (Sec I B) as the logical basis in the TN simulations. In particular, we use a unconstraint Tree TN (TTN) to represent the manybody wavefunction [72]. We adapt the TTN structure for the 2D system as shown in [73,74]. Following the description in [75], we additionally exploit the abelian U (1)-symmetry which corresponds to the total chargeQ for the TTN representation. In this way we keep the total chargeQ fixed for each simulation by choosing the proper global symmetry sector. As discussed in Sec I A the chosen local basis does naturally not respect the extra link symmetry arising from the division of the Hilbert space for each link into two half-links. Thus, additionally to the LGT-Hamiltonian H (1), we include a term to penalise the states violating the link constraint during the simulation. In conclusion, we simulate the Hamiltonian with µ ∈ {µ x , µ y }, where the penalty term vanishes when the link symmetry is respected and increases the energy for a state breaking the symmetry. Let us mention, that this additional term translates to a nearest-neighbor interaction term in the TN simulations. In theory, the penalty factor ν should be chosen as large as possible, as the link symmetry is strictly enforced for ν → ∞. But choosing ν too large leads to the optimisation focusing on this penalty term only and fails to optimise for the physical quantities. Depending on the physical simulation parameter t, m, g e , and g m the penalty factor ν has to be chosen in a balanced way, such that we are able to optimise for the physical quantities as much as for the link constraint. In fact, when choosing ν too low, we end up with a result where the state does not strictly obey the link symmetry. If ν is too large, artifacts can appear in the proposed ground state, as the penalty term can introduce local minima and thus freeze the state in the optimisation. These artifacts can either be a matter-antimatter pair for the vacuum phase or as shown on the left in Fig. 11 a matter-antimatter hole for the charge-crystal phase. As the total chargeQ is strictly conserved by the chosen symmetry sector during the simulation, there are only two ways to get rid of such an artifact. The optimiser has to either locally violate the link symmetry, or change the state at the neighboring sites together with the artifact -both of which would increase the energy in the simulation given a large value for ν.
In order to improve this approach, we exploit two different methods. First, we start with a random state, which in general violates the link symmetry. One example for this random initialisation is reported on the right side of Fig. 11. Secondly, we drive the penalty term by increasing ν after every optimization sweep. In particular, we start by linearly increasing ν, until we observe an increment in the energy which signalises that the penalty term becomes significant for the optimisation. Consequently, we switch to a quadratic tuning of ν such that in the following few iterations we increase ν slower than in the linear regime. Finally, we as well set a maximum value for ν at which we stay for the rest of the optimisation. The three different regimes of driving the penalty parameter ν are depicted in Fig. 12 showing the energy difference δe to a higher bond dimension together with ν with respect to the iterations for an exemplifying simulation.
With this driving, we optimise the random initial state in the first phase without being too strictly focused on obeying the link symmetry. This flattens the local minima arising from including the penalty in the Hamiltonian and thereby helps to converge to the global minimum. When choosing the linear tuning correctly, most physical observables are qualitatively already captured at the end of the first driving phase without strictly obeying the link symmetry. Thus the second phase enforces the link symmetry while the last phase -with a constant ν -optimises the state for the final quantitative ground state.
Although introducing the driven penalty drastically decreases the number of simulation which are stuck in artificial configurations, this can not be completely avoided. Therefore, we simulate several samples with different random initial state. From these samples, we perform a postselection and check whether the obtained wavefunctions are indeed physically correct ground states. We as well observe the typical convergence for TN with increasing bond dimension when we discard the results with artifacts. From the different samples and the convergence in bond dimension we can estimate the relative error in the energy which depending on the physical parameters typically lays in the range of ∼ 10 −2 − 10 −4 for a 8 × 8 system.

Appendix C: Perturbation theory
Here we describe the corrections to the ground state in both the two phases outlined in Section II. Let us start by considering particles fluctuations due to the presence of a small tunnelling |t|. The system has periodic boundary conditions. a. Perturbation around the vacuum state.-For m |t|, the vacuum state (with zero energy) is cor- Penalty parameter ν (red) and energy (yellow) with respect to the number of iterations for a typical LGT simulation. The energy is plotted here as a deviation δe to the ground state energy obtained with highest bond dimension available. We start with linearly increasing ν. When the energy increases we change to a quadratic driving regime with zero gradient at the transition point. Finally we reach a predefined maximum value for ν. rected by strictly local particle-antiparticle fluctuations. The first nontrivial contribution comes from a local dimer excitation as depicted in Fig. 13, whose average energy is 2m + g 2 e /2. The truncated Hamiltonian reads (a part from the sign of the tunnelling coupling, which however does not affect the results) which is (1 + 2L 2 ) × (1 + 2L 2 ) matrix. The correction to the vacuum energy is therefore b. Perturbation around the dimer state.-Smallorder tunnelling perturbations on top of the fully dimerised states are not sufficient to remove their degeneracy. The ground-state energy sector remains degenerate up to the fourth-order in perturbation theory. Here we focus on the smallest order energy corrections for one specific dimerised configuration, and consider the possible excitations as depicted in Fig. 13. We now have two different excitation sectors, depending where we remove a particle/antiparticle pairs: when the pairs is annihilated on top of a dimer, the energy cost is 2m + g 2 e /2; otherwise, when we remove a pairs in between two dimers, we have to spend 2m − g 2 e /2. The number of possible configurations of the first type coincides with the number of dimers, i.e. L 2 /2; in the other case, we have 3L 2 /2 differ-ent possibilities. The full truncated Hamiltonian is still a (1+2L 2 )×(1+2L 2 ) matrix which now reads (a part from the overall extensive constant E N ≡ (2m + g 2 e /2)L 2 /2) where also in this case the sign of t does not affect the results. The correction to the vacuum energy can be evaluated as well by solving det(H d − ε) = 0; indeed, due to the structure of the matrix, and thanks to the properties of the determinant, we found where ε − is the negative solution of Now it is clear that, when m is approaching the value −g 2 e /4, the biggest corrections, at the lower order in t, solely come from the sector quasi-degenerate with the classical dimerised configuration. The finite-size scaling depends whether the mass is approaching from above (i.e. from the vacuum phase) or form below (i.e. from the dimerised phase): in the first case 2L 2 states contributes to the energy corrections; in the second case, if g e > 0, only L 2 /2 states get involved. An energy gap |E v − E d | ∼ Lt/ √ 2 opens. Notice that, in the pathological situation where g e = 0 as well, there is no gap opening at the second order in t and therefore a sharper transition is expected.
Let us mention that, the correction to the groundstate energy coincides, as it should, with the secondorder degenerate perturbation theory. In practice, ifQ is the projector into the classical charge-crystal sector, andP = 1 −Q the projector into the complementary sector, then we may split the eigenvectors in two contributions: |E k = |φ k + |ϕ k , where |φ k ≡Q|E k and |ϕ k ≡P|E k . The eigenvalue equation (Ĥ 0 − tV )|E k = E k |E k therefore splits in two coupled equations where we used the fact that, in our caseQVQ = 0. Corrections within the degenerate sub-sector are thus given by recursively solving the following equation At the second order in the tunnelling, the dimerised sub- sector degeneracy is not lifted, and the energy changes according to Eq. (C5). Let us stress that, when deep into the charge-crystal phase, these are the dominant corrections. However, close to the classical transition, the creation/annihilation of particle-antiparticle is energetically favourable, and non-trivial corrections to the degeneracy of the ground-state energy sector are induced by fourth-order tunnelling transitions: two different classical dimerised states are coupled whenever they share at least one "resonating" plaquette, which consists in two neighbouring horizontal/vertical dimers (see the 4m mass sector in Fig. 14). This effect partially removes the groundstate degeneracy, making energetically favourable a specific superposition of different dimer states. Incidentally, let us mention that, in the thermodynamic limit, there exist classical dimer configurations, e.g. the state where dimers are all vertically (horizontally) aligned with all local electric fluxes pointing in the same direction, which are not resonating with any other fully dimerised state at any order in perturbation theory.
Appendix D: Exact results of the 2 × 2 system In the zero-charge density sector, the single plaquette system, i.e. 2 × 2, admits only 13 gauge-invariant diagonal configuration, in the Spin-1 compact representation of the electric field. The full Hamiltonian can be easily constructed by considering each mass sector {0, 2m, 4m} independently, and it acquires the following block structure, where D j = D † j , T 20 = T † 02 , T 42 = T † 24 , and all matrix entries are reals. To construct each block, we used the gauge-invariant eigenstates of the electric fieldÊ x,µ and particle numbern x , as depicted in Fig. 14.
The diagonal blocks read where I 4 is a 4 × 4 identity matrix. The out-diagonal blocks are responsible for creation/annihilation of particle/antiparticle pairs and are given by The exact diagonalisation of the HamiltonianĤ 2×2 allows us to explore the behaviour of the ground state in the vicinity of the transition m g 2 e /4. As ex-pected from the enhancement of quantum fluctuations, the gauge-invariant hopping term gets picked at the transition ( Fig. 15 top panel). The overlap of the ground state with the bare vacuum as function of m for different valus of the electric coupling is analysed as well (central panel in Fig. 15). Exact curves are compared with first-order perturbative results.
In order to explore with more care the transition region, we look at the fidelity susceptibility of the ground state [79][80][81][82], χ F (m) ≡ ∂ m GS(m)|∂ m GS(m) − | GS(m)|∂ m GS(m) | 2 , which gives the leading contribution to the ground state fidelity | GS(m)|GS(m + δ) | = 1 − δ 2 χ F (m)/2 + o(δ 2 ), since the linear contribution in δ vanishes due to the normalisation condition GS(m)|GS(m) = 1. This quantity is the perfect indicator of a changing in the geometrical properties of the ground state when varying the couplings. Moreover, from perturbation theory, it can be easily shown that χ F (m) ≤ [ GS(m)|( xn x ) 2 |GS(m) − GS(m)| xn x |GS(m) 2 ]/∆ 2 , where ∆ is the energy gap between the ground state and the lower excitations. In practice, the fidelity susceptibility of the ground state is bounded from above by the number of particle fluctuations (which is an extensive quantity) divided by the gap. Whenever χ F (m) shows a super-extensive behaviour, the ground state of the system should be gapless. From the numerical data we have confirmation that χ F (m) is enhanced in the vicinity of the transition between the two regions, as depicted in the bottom panel of Fig. 15.
Finally, in Fig. 16 we reproduce the behaviour of the matter density as function of the magnetic coupling, for different values of the electric field couplings. As explained in the main text, and confirmed by these exact results in the 2 × 2 plaquette, the local density gets enhanced by applying a small magnetic coupling; however, when g 2 m g 2 e , the particle density starts decreasing and eventually vanishing for g 2 m g 2 e . Let us stress that this phenomenon is strictly due to the finite compact representation of the gauge field.