Phase Diagram and Conformal String Excitations of Square Ice using Gauge Invariant Matrix Product States

We investigate the ground state phase diagram of square ice -- a U(1) lattice gauge theory in two spatial dimensions -- using gauge invariant tensor network techniques. By correlation function, Wilson loop, and entanglement diagnostics, we characterize its phases and the transitions between them, finding good agreement with previous studies. We study the entanglement properties of string excitations on top of the ground state, and provide direct evidence of the fact that the latter are described by a conformal field theory. Our results pave the way to the application of tensor network methods to confining, two-dimensional lattice gauge theories, to investigate their phase diagrams and low-lying excitations.


Introduction
Since the formulation of the density-matrix renormalization group algorithm by White in 1992 [1], numerical techniques based on tensor network (TN) ansätze have found widespread application in condensed matter theory [2,3]. While initially restricted to low-lying excitations of one-dimensional (1D) systems, these techniques are nowadays routinely used to investigate real-time dynamics as well as finite temperature. Over the last five years, a flurry of activity has been devoted to show how TN methods can be extended to lattice gauge theories (LGTs) [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] -i.e., to theories including dynamical gauge fields. Numerical simulations based on gauge invariant tensor networks (GITNs) have been applied to several 1D Abelian and non-Abelian LGTs, to investigate phenomena as diverse as, e.g., string-breaking and low-lying spectra (see, e.g., Ref. [19,20]). The vision behind this programme is to ultimately extend these methods to higher-dimensional LGTs displaying confinement, and thus render numerical regimes accessible, which are typically not tractable with Monte Carlo (MC) techniques. These include finite chemical potentials and real-time dynamics -describing, for instance, the time-evolution of many-body systems governed by quantum-chromodynamics [21,22].
In this work, we apply GITNs to a two-dimensional, confining LGT -a U(1) quantum link model also known as square ice [23][24][25][26]. This model represents an ideal testing ground for exploring possibilities and limitations of GITN in 2D. Its phase diagram is rich, presenting three distinct confined phases: a Néel phase, a resonating valence bond solid (RVBS) phase, and a columnar phase. From the simulation perspective, its dynamics is solvable with Monte Carlo simulations, which allows us to qualitatively and quantitatively benchmark 2D GITN approaches. Still, it is worth mentioning that an efficient Monte Carlo algorithm for this model has been devised only very recently [25], indicating highly nontrivial quantum dynamics at low energies.
Our purpose here is twofold. The first objective is to validate and assess the performances of GITNs in investigating 2D lattice gauge theories which display confinement (for deconfined theories, see Ref. [4]). The second one is to exploit the predictive power of these techniques, and in particular, the possibility of gaining novel insights based on entanglement properties -e.g., entanglement spectra and entropies -of the system.
In order to address the aforementioned objectives, we carry out extensive GITN simula-tions using the time-evolving block decimation (TEBD) algorithm [27,28] for various cylinder geometries including up to 600 spins. We start our analysis focusing on expectation values and correlations of local observables, which allow us to clearly distinguish the different phases of the model, and to locate the transition point between Néel and RVBS phases in agreement with alternative studies based on exact diagonalization. Following this, we specifically address aspects related to confinement by calculating Wilson loop expectation values, and string tensions by introducing a pair of static charges. Our results clearly show that the string tension is finite at the Néel-RVBS transition, in agreement with recent MC simulations. This first part of the analysis shows that GITN compare well to other methods in the zero-density regime, even if we are not able to reach sizes as large as the ones accessible to MC.
In the second part, we revisit the phase diagram of square ice from an entanglement perspective, allowed by the fact that GITN provides direct access to the system wave-function of both the ground state and low-lying excited states. We show how the entanglement spectrum -the spectrum of the reduced density matrix obtained by tracing out a part of the degrees of freedom from the ground state wave function -displays distinctive features in different confined phases, reflecting the corresponding symmetry breaking patterns. Then, we carry out a detailed study of string excitations in the system from an entanglement perspective. In particular, we provide entanglement-based evidence for the fact that such string excitations are described by a conformal field theory with central charge c = 1 -effectively behaving as a compactified boson. This approach has never been applied to LGTs, and, compared to other methods based on energy spectroscopy of the Lüscher term [29], it allows us to estimate the central charge from moderate system size simulations.
The paper is organized as follows. In Sec. 2, we introduce the model Hamiltonian, summarize the main features reported in the literature about the phase diagram, and provide details on our simulation strategy, based on imaginary time-evolution by means of TEBD on gauge-invariant matrix-product states (MPS). In Sec. 3, we present our results on order parameters, correlation functions, and Wilson loops, followed by numerical analysis of entanglement properties and string excitations in the RVBS phase in Sec. 4. We conclude in Sec. 5.

Square Ice Hamiltonian
The Hilbert space of the square ice model is defined -analogously to conventional lattice gauge theories -by a set of spin-1/2 variables residing on edges of a square lattice.
The system Hamiltonian reads [24,25,30,31]: where the sum goes over all square plaquettes and the plaquette operator f = σ + 1 σ + 2 σ − 3 σ − 4 + H.c. inverts the spins s 1 , . . . , s 4 along the boundary of oriented plaquettes, as show in Fig. 1(a). The first term in the Hamiltonian is the equivalent of a magnetic field term, and the second term with coupling parameter λ -also known as Rokhsar-Kivelson (RK) term [32] -imposes a finite energy on oriented plaquettes (i.e. on two out of four configurations where the product of all spin operators in the z basis is +1).  It is convenient to represent spin up (down) by right-and upwards (left-and downwards) pointing arrows on horizontal and vertical links, respectively. In this picture, f flips oriented plaquettes and projects out anything else, while f 2 counts oriented plaquettes. We study the system on finite cylindrical lattices, and label the vertices ν = nx + mŷ = (n, m) with n ∈ {1, . . . , L x } and m ∈ {1, . . . , L y }. We locate the spins s r on centres of horizontal and vertical links, r = (n + 1/2, m) and (n, m + 1/2), respectively. Spins on the open boundary at n = 0 and n = L x are fixed, and periodic boundary conditions apply in thê y-direction with m ≡ m + L y .
The Hamiltonian (1) is invariant under a local U(1) gauge symmetry generated by the vertex operator wherex,ŷ are unit lattice vectors along either dimension (see Fig. 1(a)). G ν counts the difference between in-and outbound arrows at any vertex ν, and the dynamics of the square ice decouples into sectors with different sets of background charges c ν as stated in the Gausslaw We primarily operate in absence of charges G ν ≡ 0, where every vertex joins exactly two in-and two outbound arrows, a condition termed 'ice rule' for its analogy to the hydrogen configuration in hexagonal water ice [33,34]. The ice-rule configuration space is equivalent to that of the six-vertex model (see Fig. 1(b)), and its dimension grows as W L×L with Lieb's square-ice constant W = (4/3) 3/2 ≈ 1.54 on periodic L × L square lattices [35,36]. This can be compared to W = 4 when ice rule is absent. The gauge invariant subspace states can be identified with the low-energy manifold of the spin-1/2 frustrated antiferromagnetic XXZ model on a two-dimensional 'pyrochlore lattice', in the limit of strong anisotropy J xy J z (or, equivalently, from an Ising model with the same term along the z-axis) [23,24]. The plaquette flip operator f arises from perturbative tunnelling elements between ice-states at the order of J 2 xy /J z , on top of which a tunable chemical potential on oriented plaquettes is added to arrive at Eq. (1). In this framework, spin-excitations manifest in the form of nonzero charges and are separated by a gap in order of J z . Analytical and numerical results on the dynamics of various excitations are established in [24,37].
The square ice of Eq. (1) changes into a model of densely packed dimers [32], when instead of enforcing the ice rule, one places alternating charges c (n,m) = (−1) n+m at each vertex [25]. As in the dimer model case [32], one can introduce winding numbers w x , w y that determine the flux (difference from in-and outbound arrows) through a vertical or horizontal cut of the lattice. On the cylinder, we define w y = n s (n,1/2) , and w x is fixed due to the open boundary conditions. As a consequence of gauge-invariance, these winding numbers are constants of motion and each pair of winding numbers constitutes a decoupled sector of the dynamics. Further symmetries of Hamiltonian (1) are global translation-and chargeconjugation invariances, the latter of which inverts all spin and negates any static charge placed on the lattice (see also Appendix A).

Zero Temperature Phase Diagram
The phase diagram of the square ice has been explored from numerical and analytical arguments [24,25]. At zero temperature, it comprises a RVBS phase, sandwiched by a Néel phase and a columnar phase (see Fig. 2).
For λ → −∞, the number of oriented plaquettes is maximized in the ground state. Long range Néel ordered spin-correlations spontaneously break translational symmetry. Charge conjugation, which leaves the number of oriented plaquettes invariant, is broken as well, and the existence of degenerate symmetric and antisymmetric Néel ground states indicates a Néel-phase that extends for λ ≤ λ c .
At a critical value λ c , the Néel-phase terminates in a weakly first order quantum phase transition, pinpointed by an energy level crossing of the antisymmetric Néel-state with an other excited state. From extrapolating exact diagonalization of periodic lattices including up to 64 spins, the transition has been estimated to occur in the thermodynamic limit at λ c ≈ −0.3727 [24]. Further global fits to spectral gaps, including higher excited zero-momentum states, lead to λ c ≈ −0.359(5) from up to 6 × 6 lattices [25].
For λ c < λ < λ RK , a phase of resonating plaquettes maximizes the kinetic energy in resonating between the two flipped orientations of individual plaquettes. Since resonating plaquettes repel each other onto one of two possible separate sublattices, translational symmetry remains broken, but charge conjugation is restored with two conjugation-invariant ground states. The RVBS phase includes the λ = 0 point, and extends up to the RK-point λ RK = 1, characterized by a vanishing string tension and a ground state wave function which can be written as an equal weight superposition of closed loops [23].
Finally, for λ > λ c , all configurations without any oriented plaquettes become ground states of the square ice. Those columnar configurations are characterized by long-ranged parallel alignment of spins along vertical or horizontal lines of links, with either all horizontal or all vertical lines locked in parallel alignment.
The nature of spinon-excitations (i.e. nonzero charges) is known to be deconfining in the columnar phase, including the RK-point λ RK [24]; and confining for λ < λ RK , including the weakly-first-order transition point λ c , where Monte Carlo simulations on large systems (up to 180 × 180) have revealed a small but finite string tension between static charges [25]. The fact that confinement is the only possible scenario here is related to the continuum-limit behaviour of the theory, where confinement is due to monopole contributions [38].

Mapping to Gauge-Invariant Matrix Product States
In order to approximate the ground state of the square ice with matrix product states (MPSs) [39], we map the cylindrical system onto an open chain of length L x by combining spins with the samex-coordinates into sites n = 1, . . . , L x . Spins on vertical links merge into vectors |s n := m |s (n,m+1/2) .
Spins on adjacent horizontal links to the left and right of vertices at n merge into | n := m |s (n−1/2,m) and |r n := m |s (n+1/2,m) , respectively. We later fix specific boundary conditions 1 and r Lx (see Appendix A).
A local computational basis for the MPS is constituted by states |ϕ n ∈ {| n |s n |r n }, from which we remove all configurations that do not fullfill the Gauss-law Eq. (3) at each vertex (n, m), m ∈ {1, . . . , L y }. This choice of basis immediately restricts our computations to a physical gauge sector, and is a very general recipe for tensor network simulations of lattice gauge symmetries [5,40]. However, it comes at the cost of doubly-counted horizontal spin degrees of freedom (DOFs). This redundancy must be compensated by imposing the constraints on our MPS representation, for instance by means of projection or energy penalties.
Here we follow a more direct approach, which employs the same efficient numerical structures that encode global Abelian symmetries [41], and effectively extends their field of application to lattice gauge symmetries, including non-Abelian gauge symmetries [12]. As a general recipe, it applies to our 2D scenario as well, and it is highly efficient in that it entirely eliminates the redundant DOFs by encoding all MPS matrices in form of U(1) symmetry invariant tensors [41][42][43]. In more detail, this well established technique equips matrix indices with integer quantum numbers: On physical sites, these are the eigenvalues q n of local generators g n in a global symmetry of the form Γ = n g n . Virtual bond-indices may carry any accessible numbers q n,n+1 selected by the MPS ansatz between sites n and n + 1. All non-zero matrix elements must further obey the additive fusion rule q n−1,n = q n + q n,n+1 ; other combinations are removed from the ansatz. As a consequence, the MPS is confined to a fixed global sector Γ |Ψ = q |Ψ of the symmetry, given by the difference q := q Lx,Lx+1 − q 0,1 of quantum numbers assigned to the first and last MPS bond.
We can now impose all constraints expressed by Eq. (5) in form of two global U(1) symmetries, which we generate locally from diagonal generators g The mapping b(·) assigns a unique number to each configuration of spins, for instance by binary counting. The so defined Γ (k) = n g (k) generate all constraints between even-odd (k = 0) and odd-even (k = 1) pairs of adjacent sites. This is a consequence of the fusion rules at sites n and n + 1, which for n ≡ k (mod 2) read Eq. (5) is now equivalent to fixing quantum numbers q n−1,n to the same constant on every other virtual bond. Hence, by starting from an initial MPS under this condition, and because the Γ (k) are conserved quantities, we exactly maintain the constraints throughout the entire simulation, and we have explicitly removed all redundant DOFs by discarding matrix elements that do not obey the fusion rule. We can easily incorporate further global Abelian symmetries into our ansatz (given they commute with each other). In particular, the winding number w y has been fixed to the sector of the ground state in form of an additional U(1) symmetry. Note that w x is already fixed by the open boundary condition 1 and any static charges on the lattice. This allows us to further remove from our local basis |ϕ n all configurations with incompatible net horizontal flux of arrows, and thus further reduce the local MPS dimension d n = dim{|ϕ n }. A brief account of matrix sizes, quantum numbers and their impact on computational cost is given in Appendix B, Further invariances on the cylindrical lattice, within specific boundary conditions, are combinations of charge conjugation, inversion (n → −n, m → −m) and translation (m → m + 1). While possible in principle with a non-Abelian invariant MPS ansatz, we do not explicitly encode those symmetries in our simulations. More details on the choice of boundary conditions and topological sector are provided in Appendix A.

Imaginary Time Evolution with Time-Evolving Block Decimation
When mapping the 2D lattice problem onto a chain, the Hamiltonian Eq. (1) has only nearestneighbour interactions which are well suited for time evolution with the time-evolving block decimation (TEBD) algorithm [27,28], or time-dependent DMRG [44,45]. To this end, we reshuffle the interaction terms in H into sums over plaquettes at fixedx-coordinates, and split the plaquette operator into local terms where f = σ − 3, σ − 4, σ + 1, and f = σ + 1,r σ + 2,r σ − 3,r apply at MPS sites = n and r = n + 1 separately, acting on the vertical and redundant horizontal spin DOFs (see Fig. 1) encoded in the local basis of either site. With the local Hamiltonian on the L y plaquettes between sites n and n + 1 is of the form r , a translationally invariant sum of 4L y nearest-neighbour Kronecker products.
The TEBD algorithm requires nearest-neighbour evolution exponentials, which we obtain at first order in the time step τ asŨ r . TheseŨ NN (τ ) have the form of local matrix-product operators, where left and right parts can be contracted into the MPS separately, such that their applications remain feasible and the matrices u (j) , u (j) r of manageable size. Note that even with symmetric tensors, a full expansion of exp{−τ H NN } into a tensor with four indices, as is done in standard TEBD, would quickly surpass memory allocations of 100 GB and more, due to the large local basis (see Appendix B).
We approximate ground states by means of imaginary TEBD as follows: First, we pick a maximal bond-dimension χ for the MPS and evolve an initial state with an initial time step τ 0 as long as the energy difference between subsequent time steps exceeds a preselected threshold ∆E. We then continue to evolve from that state at energy E(τ 0 ) with subsequently reduced step sizes τ i = τ 0 r i (we found r ≈ 0.7 to be a practical choice), until the differences between the so obtained final energies E(τ i ) drop below ∆E as well.
We have simulated N s = (2L x + 1)L y ≤ 590 spins, in systems of size L x × L y up to L x = 30 and L y = 10. Our results have been tested for convergence in χ → ∞ and ∆E → 0 by subsequently improving both the bond-dimension up to χ = 2000 and the convergence threshold down to ∆E = 10 −12 , with final time steps as small as τ ≈ 2.3 × 10 −5 . Exemplary results from this process are shown in Appendix C.
We note that TEBD is just one of many well established methods for MPS ground state search, with alternatives such as DMRG [2]. While all of these methods are ultimately limited by the sizes of bond-and local dimensions, and are sensitive to the energy spectrum in some way, they differ in how local updates of MPS matrices are carried out. In the case of TEBD, by dropping the imaginary unit, we have immediate access to the real-time dynamics.

Results: Order Parameters, Correlation Functions, and Wilson Loops
In this section, we discuss the phase diagram of the model, utilizing diagnostics based on order parameters, correlation functions, and Wilson loops. This analysis serves as a benchmark for our approach, by comparing known results to our extrapolation of the transition points, and by evaluating approximate string tensions from the decay of Wilson loops in real space. Furthermore, by investigating different aspect ratios L x /L y as shown below, this comparison allows us to identify the best scaling regime to recover correct results by means of finite-size scaling.

Phase Transition from Néel to Resonating Valence Bond Solid
We first measured the ground state spin correlations, and evaluate the Néel-order parameter where r, s go over spin coordinates and N 1 = N 2 s − N s counts the number of expectation values summed up. Similarly, we detect the RVBS phase from staggered correlations between pairs of adjacent spins: where ν, µ are lattice vertices, andî, ∈ {±x},k,l ∈ {±ŷ} are unit vectors. The sum goes over all N 2 = (4N ) 2 ν − 4N ν contributions, given four diagonal pairs at N ν = L x L y vertices,   [25] (grey with uncertainty) and [24] (black).
with ν = µ orî = ork =l. In order to reduce the impact of the open boundary, we further restricted the analysis to spin between sites n L = b and n R = L x − b + 1 with b = 1 or 2.
With these parameters, we can clearly distinguish both phases as shown in Fig. 3. The underlying two-point correlation functions are shown in Fig. 4, and already provide a qualitative picture of the different phases: In the top two panels, one can distinguish the very different behavior of the σ z r σ z s correlation function, which is long-range ordered in the Néel phase, and short-ranged in the RVBS. The latter is instead characterized by the correlations between spin pairs, which indicate flippable plaquettes on a sublattice, as evidenced by Fig. 4(d).
Order parameters for the broken translational symmetry in the RVBS phase can also be derived from staggered flippability [24,46] where the summations go over all plaquettes (and ) with respective lattice coordinates n, m (and n , m ), normalized by N 3 = N and N 4 = N 2 − N for N = (L x − 1)L y plaquettes. While the first order expectation works in our case as order parameter due to explicit symmetry breaking by boundary conditions, the second order correlations provide a more pronounced result.
We note that the estimates at ratio L x /L y = 1 display severe finite size effects, while geometries with L x /L y = 2 and 3 provide much better agreement between the different order parameters, as they optimally combine large system sizes with modest anisotropies.
From the differences between the two smallest ratios with maximum lengths L x = 10 and L x = 20, we can estimate our uncertainty due to the limiting maximal system width L y = 10. Increasing the latter should however have a much smaller effect due to the periodic boundary conditions.

Confining Properties: Wilson Loops and String Tension
We now turn to a diagnostic which is directly related to confinement -i.e., Wilson loops [49]. In typical lattice gauge theory formulations in Euclidean space-time, Wilson loops are defined as the product of parallel transporters over a spatio-temporal closed path. Differently from traditional Monte Carlo simulations of lattice actions, wave-function based methods such as GITN can only access real-space Wilson loops. As such, extracting the string tension from the latter comes with a caveat, as these two types of Wilson loops coincide only in the presence of Lorentz invariance. We comment on this aspect at the end of next section.
In our simulations, we measured the expectation values of the operator that flips oriented loops of spins on the boundary of a rectangular area A with n ∈ [n 1 , The connection to the more traditional definition of the Wilson loop is immediately established considering that the present model is a quantum link formulation of a compact U(1) gauge theory [25,40]. For a single plaquette A = , one recovers the plaquette operator f . In order to smoothen out corner contributions (which will lead to additional corrections due to ultra-violet effects), we tried to employ Creutz's ratios [50]. However, we found out that the latter did not facilitate our analysis, probably due to the fact that we could only access modest sizes of the area enclosed by the Wilson loop itself.
Instead, in systems of dimension L x = 3L y , with L y ≤ 10, we fit the expectations W (A) = Ψ GS |f A |Ψ GS to an exponential decay with the enclosed area where A = ∆ x ∆ y is of width ∆ x = n 2 − n 1 and height ∆ y = m 2 − m 1 , and both the area contribution α ∆y as well as β ∆y are free parameters that may depend on L y and ∆ y (besides λ) in order to account for finite size corrections and the sharpness of the boundaries. Typical results for two points (λ = −0.4 and 0.6) are shown in Fig. 6(a).
Overall, we find that the numerical results are very well captured by Eq. (15) even for modest areas. For the accessible system sizes, the dependence of α ∆y and β ∆y /∆ y on ∆ y appears to be of the type c 0 + c 1 atan(c 3 /∆ y ) with some coefficients c k . Remarkably, the exponential decay of the Wilson loop seems to be reproduced extremely well even when the value of the correlation itself is smaller than our truncation error. This indicates that, given a certain truncation error, the corresponding error in the Wilson loop correlation can be considerably smaller. Nevertheless, the slopes α ∆y remain affected by simulation precision and we extrapolated our results in both χ and ∆E (see Appendix C). The reason is that even modest variations of the Wilson loop correlator affect the estimate of α ∆y .
In order to characterize the area contribution α ∞ for large ∆ y and in the thermodynamic limit (TL) L y → ∞, we fixed the rectangle width at half the system size, ∆ y = L y /2, and plotted α Ly/2 over λ ∈ [−0.6, 1] for L y = 4, 6, 8 and 10 (see Fig. 6(b)). Apparently, the value of α Ly/2 generally decreases with L y -most prominently around λ c and λ RK due to finite size effects, while the decline occurs at a slower rate in the bulk of the RVBS where we expect a finite α ∞ > 0. Even though we cannot conclusively determine the asymptotic behavior from the available system sizes, we estimated α ∞ from linear extrapolation to 1/L y (black crosses in Fig. 6(b)). As expected, the value of the area contribution approaches zero when moving towards the deconfining RK-point λ RK = 1. Moreover, we observe a visible drop of the decay coefficient for L y > 6 in the almost-critical region around λ c , while the TL estimates suggest a finite value well above zero -testifying that, at a critical point, there is a finite mass gap.

String Tension from the Potential Energy between Static Charges
A direct estimate of the string tension is provided by introducing two static charges in the system, and monitor the amount of excess energy, compared to the ground state in absence of charges, as a function of their distance from each other. In our simulations, we considered a lattice with two static charges c + = +1 and c − = −1, placed at the sites n ± = (L x +1)/2± /2 in horizontal distance of = n + − n − lattice spacings. This configuration has been chosen in order to alleviate potential boundary effects (both charges are at the same distance from the boundaries). At the MPS level, these static charges can be readily implemented by changing the value of the Gauss-law at a given site.
In Fig. 7(a), we show a typical result for the local energy distribution in the cylinder along thex-axis. The plot already shows how the insertion of static charges significantly affects the expectation value of the Hamiltonian density, not only close to the charges themselves, but in the entire region between the two.
The string tension is evaluated by considering the energy difference between the system ground states with (E c ) and without (E 0 ) static charges, and by plotting it as a function of the intercharge distance : Two typical sets of results are presented in Fig. 7(b): Already for modest values of , the linear dependency captures the numerical data well.
The fitted values of σ are shown in Fig. 7(c) for the parameter regime λ > 0.2. As expected, the string tension vanishes at the RK point, and increases almost linearly within the RVBS phase. However for λ 1, finite size effects, that confine the string and lead to self-interaction, become visible with nonlinearities in Eq. (16) and values σ start to drift with increasing system size. Eventually our system sizes become too small for a direct measurement of the string tension for λ 0.2, including λ c . We note that the string tension results are in reasonably good qualitative agreement with the values of α extracted from the expectation values of the Wilson loops in the previous section. This indicates that, despite the RVBS order breaks translational symmetry at the lattice spacing level, the breaking of Lorentz invariance affects the expectation value of the Wilson loop, but not its parametric dependence on λ. At the quantitative level, it is instead impossible to draw conclusions on this parallelism: The reason is that, as discussed in the last subsection, the system sizes we simulated are too small to properly extrapolate α Ly to the thermodynamic limit.

Entanglement Properties of the Ground State and String Excitations in the Resonating Valence Bond Solid Phase
Tensor network simulations provide direct access to the entanglement properties of the ground state and of the low-lying states of the theory. In the context of the quantum ice model, these properties are essentially unexplored, yet intriguing: Given the fact that this model describes a gauge theory, its Hilbert space cannot be written as a tensor product for bipartitions of any size while upholding the Gauss-law.
In this section, we first discuss the entanglement spectrum (ES) properties of the RVBS phase, and then turn to the entanglement entropies of string excitations.

Entanglement Spectrum in the Resonating Valence Bond Solid Phase
We characterize the entanglement properties of the system by investigating the reduced density matrix of a bipartition A: obtained by considering the ground state wave function, and taking the partial trace with respect to the complement of A, which we denote as B.
In our MPS representation, we have direct access to the ES in form of the singular values {s k } associated with a cut through the cylinder between any sites n and n + 1. A sample of our results from the largest system size we simulated is shown in Fig. 8, where the cut is made at an off-center MPS bond n < L x /2 in order to reveal degeneracies in the spectrum 0 I IIa IIb B A Figure 9: Spin configurations s (m,n+1/2) , m = 1, . . . , L y at a vertical cut, associated with the lowest levels 0-IIb of the entanglement spectrum in the RVBS (see Fig. 8). Oriented plaquettes on the resonant sublattice (green) can be destroyed (red) by virtual plaquette-flips in the opposite sublattice (red dots).
that are not merely caused by lattice(-inversion) symmetries. We can characterize each level by the spin configuration on the horizontal links joining the bipartition. According to this criterion, we have grouped states 0-IIb (see Fig. 9): Level 0 comprises all d 0 = 2 Ly/2 configurations available under the condition that spins are anti-aligned on the sublattice hosting the resonant, oriented plaquettes. A subsequent flip of a single pair of antiparallel spin on any of the L y /2 plaquettes of the opposite sublattice leads to the boundary configurations of level I, which is d I = 2 Ly/2−2 L y degenerate. This subspace describes 'excitations' on top of the RVBS state.
In levels IIa and IIb, two such defects occur, reducing the number of antiparallel alignments on the resonant sublattice to L y /2 − 2 or L y /2 − 4, depending on whether the defects are adjacent or not. Corresponding degeneracies are d IIa = 2 Ly/2−3 L y and d IIb = 2 Ly/2−4 L y (L y /2−3), respectively. The drop in singular values with each level is roughly proportional to the number of defects.
This spectrum gives a direct signature of the RVBS state. In particular, all counting above just reflects the structure of the symmetry broken state, and excitations on top of that. In principle, one can try to apply finite-size scaling analysis to extract the critical point λ c from the ES: However, our results suggest that the latter suffers from more severe finite size effects than the correlation functions, rendering this procedure rather inaccurate.

Entanglement Entropy of String Excitations
String excitations on top of the vacuum play a paradigmatic importance in lattice field theory. Historically, they have been vastly employed to provide direct information about confining properties of a theory (without the need of adding dynamical charges) [51]. In the context of magnetism, string excitations have been widely discussed as a distinctive feature of classical spin ices [52,53], and some of their properties have been recently observed [54,55] and discussed at the quantum level [56,57].
While string excitation are typically addressed employing diagnostics based on energetics (such as, for instance, the analysis of Lüscher's terms in the string potential), we employ here a direct investigation of the string entanglement properties to show how, in the RVBS phase, string excitations are well captured by a conformal field theory with central charge c = 1.
In order to estimate the central charge, we have evaluated the entanglement properties in form of the bipartite von Neumann entropies S ±1 (n, ) in presence of charges ±1 in distance , from the MPS bond between sites n, n + 1 (same configuration as the one employed to investigate the string tension). From these results, we have subtracted the background values from simulations without charges S BG (n). The difference S diff (n) = S ±1 (n, ) − S BG (n) is then interpreted as the string entanglement entropy. Finally, we fit our numerical results with [58] S fit (n) = where the free fit parameters are c (central charge) and offset S 0 . Two typical sets of result are shown in Fig. 10(a). Despite the small number of points available even for our largest system size, we observed an accurate fit of our data (in one-dimensional systems, Eq. (18) is sometimes not accurate at distances of order of the lattice spacing due to strong UV corrections). When approaching the critical region λ c with λ ≤ −0.2, these results tend to become unreliable and deviate with increasing system size, as the string excitations start spreading over a region of space in theŷ-direction which is of the order of L y . Instead, deep in the RVBS phase, we observe that, for our largest system size, the string width is sufficiently small such that self-interaction effects are negligible. In this regime, the value of the central charge (blue dots in Fig. 10) is close to 1 within the fit error (which is of order 10%, including the bare fit error and an estimate due to the few points available). This observation is consistent with a conformal field theory description of the string excitations in terms of a compactified boson theory. While this behavior, in the context of non-Abelian gauge theories, has already been widely discussed [29,59], this was so far only diagnosed based on long-distance behavior of particle-antiparticle potentials. Here, we have instead provided an entanglement-based proof of this behaviour, that has the key advantage of requiring more modest system sizes to be precisely diagnosed.

Conclusion
We have reported a gauge invariant tensor network investigation of square ice -a U(1) quantum link model in two spatial dimensions. Our results on the phase diagram are quantitatively consistent with previous results in the literature, and represent an important systematic benchmark for the accuracy and reliability of tensor network methods applied to confining lattice gauge theories in more than one dimension.
In addition to diagnostics based on correlation functions and the excitation spectrum, we have investigated the entanglement properties of the system. First, we have studied the entanglement spectrum of the ground state wave function: While being too finite-size sensitive to be employed for estimating phase transition points, it reflects the characteristic symmetry breaking pattern of the RVBS phase. Secondly, we have provided evidence that the string excitations in such models behave as bosonic strings, using diagnostics based on the entanglement entropy. This prediction can be in principle tested in quantum simulators for quenched lattice gauge theories [46,[60][61][62], and testifies the fact that universal properties of gauge theories (such as the conformal description of string excitations) can be observed in such systems despite achieving modest system sizes. Our work indicates that, whilst challenging from a computational physics perspective, tensor network methods for the simulation of confining lattice gauge theories can reach system sizes which allow a distinct, non-ambiguous characterization of phases of matter. While our study here was on a quenched gauge theory for which accurate Monte Carlo simulations can be carried out -so that we could have a direct benchmark -, the tensor network methods we employ are sign-problem free, and can be applied to real-time dynamics and to theories in the presence of fermionic matter. As a direct extension of this work, it would be interesting to study the effect of dynamical matter in U(1) quantum link models, which is also connected to recently formulated problems in the context of strongly correlated electrons [63,64]. From the methodological perspective, testing the performance with more specialized and powerful tensor network classes such as projected entangled pair states or two-dimensional tree tensor networks appears promising in view of the entanglement results we presented [65][66][67][68][69].
. . . and lattice-inversions (I x : n → L x − n + 1 and I y : m → L y − m + 1). With 'aligned' BCs, permissible operations are for instance P 1 = C h ×I y , P 2 = C v ×I x and P 3 = C ×T y , while 'antialigned' conditions instead support Q 1 = P 1 , Q 2 = C × P 2 and Q 3 = P 2 3 . The latter explicitly breaks the (otherwise spontaneously broken) translational symmetry that transfers resonant plaquettes of the RVBS from one sublattice onto the other: Namely, oriented (and resonating) plaquettes may fully cover only one of the two sublattices, while the other sublattice can not be oriented fully (with at least ∝ L y defects at the boundary, as depicted in Fig. 11(a). As a consequence, resonances on the former sublattice are preferred without frustrating the resonant structure, and the finite size gap is increased. Conversely, the two symmetry breaking Néel configurations remain coupled by Q 2 . Instead, in the Néel phase λ λ c the gap can only be increased by selecting 'aligned' BCs, under which no such transition is permitted and in similar fashion, only one of the two conjugate Néel-configurations can maximize the number of oriented plaquettes.

B Performance: Local-and Bond-Dimensions
Our implementation of the TEBD algorithm, as outlined in Sec. 2.4, performs local updates in O(χ 3 dK 2 ) + O(χ 2 d 2 K)-time and consumes primary memory of order O(d 2 K), where the factor K = L y + 1 originates from the first-order truncated Taylor-series of the imaginary evolution exponential in Eq (9). For our linearized two-dimensional systems featuring large local dimensions d ≥ χ, this is an essential advantage over standard TEBD, which evolves with full nearest-neighbour exponentials at O(χ 3 d 3 ) + O(χ 2 d 4 )-time and in our case requires impermissible large amounts of memory, scaling as O(d 4 ) [27,28]. The lower order evolution errors (at worst O(τ 2 )), on the other hand, remain fully under control by extending the total simulation run-time and by actively monitoring convergence against ∆E while dynamically reducing the time step τ when necessary. As a consequence, for system sizes where we could still employ both TEBD implementations (up to L y = 6), we obtained ground state approximations with comparable accuracy. Further performance improvements are enabled by the constraints inherent to the square ice model. Most importantly, the local dimension d of the MPS grows exponentially in the system width L y , and eventually limits our simulations to about L y ∼ 10. However, in comparison withd = 4 Ly of an ordinary spin 1/2 system without constraints, we achieve a reduction to d ≈ 0.53 × 2.84 Ly due to the locally encoded Gauss-law and the selected topological sector w x = L y /2 (see Fig. 12).
Note that our local growth rate 2.84 is still larger than Lieb's result W ≈ 1.54 because the ice rule between adjacent sites n, n + 1 is not implemented in the local Hilbertspace, but instead imposed through additional constraints mediated by global U(1) symmetries over MPS bonds. Together with the winding number w y , these quantum numbers are carried through the MPS bond-indices and the overall bond-dimension χ can be decomposed into blocks of indices carrying equal tuples of quantum numbers q (which in our case are q (0) and q (1) , as detailed in Sec. 2.3, along with a third number keeping track of w y ).
Sizes of these blocks χ q are shown in Fig. 13 in final simulation states for various λ. While in general, block sizes depend on the current simulation state, we observed that they very quickly freeze out in the final state configuration, in which they remain until the end of the simulation. We additionally found that, within the simulated parameter regime, even the largest block was more than one order of magnitude smaller than the overall dimension χ, and block sizes χ q attain a rather flat distribution when approaching the RK-state at λ → 1. This feature immediately reflects the nature of the state at the transition point, which is made of an equal weight superposition of all possible closed loop coverings.
We also point out that because every local basis state carries a different set of quantum numbers (owed to the encoded spin configuration), matrices decompose into blocks with physical dimension d q ≡ 1. From a performance point of view, we benefit greatly from these small block-sizes, since the dimension χ q determine the actual matrix dimensions in the TEBD linear algebra manipulations. Roughly, memory and runtime complexity reduce by a factor comparable to the relative size of the largest block encountered [42,43].

C Convergence
The precision of our simulation results depends both on the MPS bond-dimension, χ and the simulation convergence threshold ∆E which controls the size and number of imaginary evolution time steps (see Sec. 2.4). We have repeated simulations over a wide range in both parameters in order to check the validity of our simulations and to obtain reliable error estimates. Exemplary results are presented in Fig. 14.
While all results converge asymptotically in both χ → ∞ and ∆E → 0, certain order parameters and Wilson loops with small enclosed area are predominantly susceptible to the energy threshold (see, for instance, Figs. 14[(b),(d)]). In those cases, contributions to excited states have not sufficiently been suppressed, and very small timesteps are eventually required to cope with the second-(and higher) order errors in the time step. These errors then manifest, for instance, in the form of anisotropies in Wilson-loop measurements due to an unresolved finite-size gap (as discussed in Appendix A). As a counter-measure, we average our measurements over both even and odd sub-lattice positions. We also found that certain simulation results (such as energy and order parameters) can be well approximated by the following two-dimensional function O(χ, ∆E) =Ō + c 1 e −|c 2 |χ + c 3 ∆E −|c 4 | In situations where our data O(χ, ∆E) (typically around 10 datapoints) was neither converged satisfactorily in χ nor ∆E, we thus found the unknown extrapolated valueŌ among constants c 1 , . . . , c 4 (of which c 1 and c 3 can have different sign, as exemplified by Figs. 14[(b)-(f)]).