Lattice gauge tensor networks

We present a unified framework to describe lattice gauge theories by means of tensor networks: this framework is efficient as it exploits the high local symmetry content native to these systems by describing only the gauge invariant subspace. Compared to a standard tensor network description, the gauge invariant model allows one to increase real and imaginary time evolution up to a factor that is square of the dimension of the link variable. The gauge invariant tensor network description is based on the quantum link formulation, a compact and intuitive formulation for gauge theories on the lattice, which is alternative to and can be combined with the global symmetric tensor network description. We present some paradigmatic examples that show how this architecture might be used to describe the physics of condensed matter and high-energy physics systems. Finally, we present a cellular automata analysis which estimates the gauge invariant Hilbert space dimension as a function of the number of lattice sites that might guide the search for effective simplified models of complex theories.


Introduction
In modern science gauge theories play a fundamental role, with examples ranging from quantum electrodynamics to the standard model of elementary particle physics [1,2]. They represent a cornerstone in our understanding of the physical world and lie at the heart of theories dealing with such diverse systems as quantum spin liquids and the quark-gluon plasma. Starting from the breakthrough contribution by Wilson in 1974 [3,4], lattice gauge theories (LGTs) have attracted significant attention across several branches of theoretical physics. While the lattice formulation of gauge theories has intrinsic fundamental interest in high-energy physics due to the prominent role played by gauge fields, emergent gauge models have also been introduced in different condensed matter setups in relation to exotic many-body phenomena such as quantum spin liquids and topological states of matter [5][6][7]. Furthermore, there is good reason to believe that certain types of gauge structures could open up new possibilities for quantum computation [8,9].
Quantum simulation of gauge theories is receiving an increasing degree of interest [10], due to the fact that these types of platforms could provide the tools to simulate dynamical properties and/or out of equilibrium physics in lattice gauge models including fermionic matter fields, as they are sign-problem free simulators by construction. In the context of atomic, molecular and optical physics, several proposals for the quantum simulation of LGTs have been made recently [11][12][13][14][15][16][17][18][19]. These quantum analog simulations can be seen as a complementary tool to the existing classical ones-the latter could be used to benchmark the outcomes of the former.
The prediction power of quantum gauge theories is often limited due to the fact that the computational resources needed to perform simulations involving large numbers of particles are typically forbidding. The present manuscript presents a partial solution to this problem by showing that the gauge invariant subspace of these theories can be represented exactly by a local set of tensor networks, thus it is possible to apply the well developed and successful architecture of tensor networks to LGTs. In other words, we present an exact and efficient tensor network representation of abelian and non-abelian LGTs [20].
Tensor networks are one of the mainstream paradigms for simulating quantum many-body lattice systems, both in and out of equilibrium, via a representation of the quantum state with tailored variational ansatz wavefunctions. They originated from the understanding that the density matrix renormalization group (DMRG) technique [21] could be recast in a variational formulation by means of matrix product states [22][23][24][25][26][27][28][29]. This stimulated the further development of such a framework in the last decade, extending the tensor network paradigm to encompass higher dimensionality [30], peculiar geometries [31,32], and the limit to the continuum [33].
One of the most appealing features portrayed by tensor networks is the possibility to encode and control global symmetries for the local degrees of freedom [34,35] that characterize several condensed matter models. In fact, a general, robust and numerically efficient formulation of any such symmetries in the tensor network framework is known [36,37]; it is commonly used in simulation to achieve an enhancement of the algorithm performance, as well as precise targeting of irreducible representation sectors [38][39][40].
Lattice gauge symmetries differ from global ones, since they have quasi-local supports and are typically homogeneous, yielding a combined Lie algebra of generators which grows extensively with the system size. Nevertheless, several physical contexts have been found where tensor networks are an exact description of the ground states of gauge-invariant Hamiltonians (e.g., 2D toric code that is an Ising gauge theory [8,41,42]). More recently, this framework has been successfully applied to LGT related problems [20,[43][44][45][46][47][48][49]. In fact tensor networks represent microscopically the local Hilbert spaces and at the same time are tailored on a real-space wave-function representation, so they can be used to describe real-space locality and local symmetries altogether.
Here we show how tensor networks can exactly encode lattice gauge symmetries providing an architecture that is completely general and computationally efficient: our approach outperforms straightforward approaches that do not explicitly exploit gauge symmetries. To achieve this goal, the use of alternative formulations of gauge theories are highly desirable, with the principal motivation being the identification of models with a finite dimensional Hilbert space at each link or site which can be simulated by tensor networks algorithms. Thus, we develop this architecture in the quantum link model (QLM) formulation [50][51][52]

of Hamiltonian
LGTs. Wilsonʼs formulation of LGT has an infinite-dimensional Hilbert space at each link due to the use of continuously varying fields [3]. QLMs provide a complementary formulation of lattice gauge theories introducing generalized quantum spins associated with the links of a lattice. In fact, under some physically motivated assumptions, Wilsonʼs lattice gauge theories can be obtained from QLM [53,54].
One possible way to approach the continuum limit for QLMs is by dimensional reduction from a higher dimension where continuous gluon fields arise as collective excitations of discrete quantum link variables, in the same way as magnons arise as collective excitations of quantum spins [16]. Such an extra dimension is expected to be exponentially smaller than the actual spatial dimension, thus not posing a threat to numerical methods. A second strategy is to increase the number of rishons per link, which also have been seen as a way to achieve the continuum limit without requiring the extra dimension [17] (we will investigate numerical feasibility of this strategy later on).
In addition there are several examples of condensed matter models, characterized by lattice gauge symmetries, where the gauge degrees of freedom are inherently finite-dimensional. This is the case, for instance, for spin-ice or quantum dimer models [55] or in discrete gauge models like the Ising gauge theory [41].
We present the formulation of the LGT network in detail, allowing one to represent efficiently and exactly the gauge constraints of this class of systems, with a performance that improves up to quadratically with the quantum link dimension, thus increasesing its efficiency at the Wilson limit.
The paper is organized as follows. In section 2 we review the framework to describe lattice gauge theories into quantum link formulation. In section 3 we provide a constructive scheme to embed the QLM within the tensor network framework, which relies on matrix product formalism in 1D and projected entangled pair formalism in higher dimensions. The algorithm to exploit such a representation in a numerical context is described in section 4, mainly focusing on time evolution (both in real and imaginary time). In section 5 we investigate theoretical scaling of effective Hilbert spaces growth, under the QLM constraints, made easily available with the tensor network model. Finally, in section 6 we draw our conclusions.

Quantum link models
From now on, as we focus on numerical simulations, we assume that the space of the gauge degrees of freedom is finite dimensional. Starting from this assumption, the formulation in quantum link model language of lattice gauge theories follows without additional loss of generality [50][51][52]. We define the gauge invariant model of interests by defining three elements: x x ]-we describe as quantum degrees of freedom both the lattice sites, which we will refer to as 'matter field', and the 'gauge field' and its canonical conjugate variable or 'electric field' located on the links (the lattice bonds between neighboring sites, every link being shared by a different pair of sites).
• The gauge symmetry generators [ ν G x ]-unlike global symmetries, which operate upon the whole lattice, gauge symmetry generators have a localized support, each one involving a single matter field site and all the gauge fields connected to it.
• The gauge invariant dynamics [H]-the dynamics are defined via a Hamiltonian which commutes with the whole algebra of gauge generators, which guarantees that gauge invariance is conserved throughout the time evolution (as in figure 1, panel a).
In this section, we analyze in detail these elements in a quantum link formulation, while stressing the connection to typical lattice gauge theory models.

Local degrees of freedom
As mentioned, there are two types of degrees of freedom in lattice gauge models, which we describe as finite-dimension quantum variables: • Matter fields ψ x are located on the vertices of the lattice x. They are usually fermionic fields that describe the 'quarks' of the model, ψ ψ δ = { , } x y x y † , . They can also be bosonic fields describing, for instance the Higgs field. In non-abelian models, fermions ψ a x carry color degrees of freedom a. For example, in U (2) x exist on the links of the lattice μ 〈 + 〉 x x , x . They are bosonic fields that describe the gauge bosons of the model. We use the quantum link formulation to recast these fields as bilinear operators: = x , as sketched in figure 1, panel d.
In the QLM formulation, the gauge boson is split into a pair of rishons, linked together by a U (1) symmetry constraint.
As discussed in [16,53], such bilinear decomposition is well-defined once a representation of the symmetry group has been selected for the gauge boson U ab . As a result of this formulation, every lattice link now hosts two field modes, typically called 'rishons' in the usual terminology of QLMs, and respectively labeled as μ 〈 + 〉 x, x . The meaning of these auxiliary modes will become clear when we elaborate on some particular cases and models, nonetheless we advance that they can be seen as a generalization of the Schwinger representation for the gauge field U ab .
Such bilinear representation of the gauge fields can be made either fermionic or bosonic, by setting the appropriate commutation relations for these operators .
on every link is a conserved quantity. This is due to the fact that the rishon degrees of freedom μ c x a , x appear both in the gauge symmetry operators ν G x and in the Hamiltonian H only via x and by construction = . In other words, in the QLM formulation of lattice gauge theories an additional, artificial local symmetry arises: the conservation law of the total number of rishons on a given link, which is always U (1) symmetry generated by x (regardless of the symmetry group generated by ν G x which may as well be non-abelian). There are different representations of the same symmetry depending on the number of rishons per link N one selects. In any case, we restrict the Hilbert space to the 'physical' states φ 〉 | phys which satisfy x . For simplicity, we will refer to this symmetry selection rule as link constraint, as opposed to the gauge constraint which is generated by ν G x instead (see next paragraph). With a little abuse of notation, in cases where the total number of rishons on a link x is independent of the link itself (i.e. homogeneous and isotropic QLM) we will sometimes omit the link label subscript.

Local generators of the gauge symmetry, and gauge constraint (Gauss' law)
The gauge symmetry is defined via the set of its generators ν G x : they all commute with the Hamiltonian = ν H G [ , ] 0 x , and have localized support. To properly characterize the generators ν G x , it is convenient to define the elementary transformation on the gauge fields beforehand. We will separately consider the abelian and non-abelian parts U (1) and SU(N) respectively as follows.
• Abelian U (1): here the elementary transformation is generated by the difference of the rishon occupation numbers on the same link, i.e. = − x , which plays the equivalent role of the electric field in quantum electrodynamics. Its action on the gauge field changes the field with a phase, or infinitesimally = , , , , , , , with μνω f the structure constants of the SU(N) algebra and λ λ δ = μ ν μν Tr ( ) 2 .
Having properly defined the elementary transformations of the gauge fields, we can now easily introduce the complete gauge symmetry generators.
• The local generator of the U (1) part of the gauge model is defined by , , where the first line expresses the generator in terms of the electric field component and the second line in terms of matter and rishon fields. Thanks to the the QLM formulation, it is possible to write G x as an operator acting only on the QLM degrees of freedom on vertex x. Around every vertex, the gauge-invariance (or gauge-covariance) constraint is given by ⎡ which is equivalent to Gauss' law, and the physical states φ 〉 | phys are those which satisfy it. • The generators for the non-abelian SU(N) part of the gauge transformations fulfill the usual algebra Again, in QLM formulation ν G x acts on matter and rishon fields belonging to lattice site x only. The gauge invariant subspace corresponds to the trivial irreducible representation subspace of the symmetry group generated by ν G x , i.e. the singlet subspace: . This provides an extension to the non-abelian gauge symmetries of Gauss' law.
It is important to stress that every single element of the algebra ν G x is local as it acts nontrivially only on the degrees of freedom sharing the vertex x, as this will be the key ingredient for the computational improvement. Such a gauge symmetry locality is sketched in figure 1.b: every gauge generator acts only on one matter field site ψ x a and z gauge field sites x (or rishon sites μ c x a , x in QLM formulation), with z being the lattice coordination number, i.e. all the quantum degrees of freedom sharing vertex x.
As a result, the total number of generators ν G x in the whole lattice gauge algebra scale extensively with the system size, a property which dramatically reduces the manifold dimension the system lives in, as we will see later on.
At the same time, the combined gauge invariance constraints acting on a given vertex x, can be assembled into a single linear mapping, which reads x r s s s s x j x , x , have been chosen. This mapping defines a 'reduced' local basis 〉 j | x r which spans exactly and solely the local gauge-invariant subspace. In the next sections, the set of states 〉 j | x r will be adopted as the logical or computational basis for all numerical purposes, and will be the starting platform for building a tensor network formulation.

Gauge invariant dynamics
The last element that has to be to defined is a gauge invariant model, its dynamics formulated via the Hamiltonian H. By construction, a gauge invariant Hamiltonian must commute with the local generators of the gauge symmetry and those of the link symmetry in the QLM formulation, i.e.
Clearly the class of Hamiltonians satisfying these requirements is still extremely wide. Here we will focus on short-range Hamiltonians that encompass the physics of typical lattice gauge models. A pure gauge model, which embeds non-abelian gauge symmetry content, is given by the following Hamiltonian: The electric terms quantify the energy of the flux for the abelian or non-abelian part of the gauge group. While the first term infers zero abelian electric flux on the link, the second favors singlets of rishons in the non-abelian color variables. The magnetic term associates a positive energy density to every non-zero magnetic flux on every plaquette. The terms g abel 2 , − g non ab 2 and g magn 2 are the coupling constants for the abelian part of the electric field, non-abelian part and the magnetic term, respectively. A physically meaningful choice of constants is the one that recovers the Kogut-Susskind (KS) Hamiltonian [4]: that is , where a is the lattice spacing. Indeed, with this special choice of the couplings, one expects to recover the physics of the usual U(N) gauge theories in the continuum limit. Alternatively, one expects to approach the strong coupling limit by setting non ab magn . It is important to remark that the previous quantum-link Hamiltonian satisfies a U(N) gauge invariance by construction: it is however possible to reduce the gauge symmetry into a pure SU(N) by adding artificial Hamiltonian terms which explicitly break the U (1) part of the gauge symmetry, as proposed in [16,53].
The coupling of the gauge fields with the matter fields is done with the lattice version of the 'minimal' coupling, i.e. a hopping term of fermions mediated by the gauge field. Also, the mass term of the fermions is a gauge invariant term, hence, where we have defined site dependence hopping constants μ J x, x and mass term m x , in case a specific distributions of signs, depending on the sites, is needed for a particular type of fermion introduced on the lattice. This type of minimal coupling is also sketched in figure 1, panel c.

Examples
We have presented all the ingredients that are necessary to define a quantum link version of a lattice gauge theory, however for the sake of clarity and concreteness, we now present four particular examples: the simplest + (1 1) dimensional QLM with the abelian U (1) symmetry, the simplest + (1 1) dimensional QLM with non-abelian U (2) symmetry, and applications for two relevant models in condensed matter physics: quantum dimer [6,59] and spin-ice [56,57] models on the square lattice [60].
U(1) quantum link model-the gauge invariant quantum Hamiltonian is given as where the last term is a staggered chemical potential profile for the matter field which is a spinless fermion field ψ ψ Here J is the strength of the matter-gauge field coupling, g 2 the electric-field energy density and m the staggered mass. The gauge fields can be written in terms of rishons = , † , , . The two independent local symmetries in this U (1) QLM are: x appears because we introduce ψ x spinless fermionic operators (matter fields with a staggered mass term m) usually denoted as staggered fermions [4,61]. The vacuum of the staggered fermions is given by a quantum state at half-filling describing the Fermi-Dirac sea.
In what follows, we aim to understand in more detail two limits dependant on the occupation N . Thus, we characterize the action of the gauge operators and electric field operators on a Hilbert space defined by the occupation of rishons + n x, and + − n x 1, or equivalent by the total number of rishons on the link and the electric flux = |¯, , where we have omitted the labels of the link 〈 + 〉 x x , 1 : [17]: Wilson formulation of compact U (1) gauge theories start with an infinite local dimensional Hilbert space defined with two conjugate variables: the electric field E and an angle ϑ that fulfil the usual commutation relation of position and . Then by defining the link operator = ϑ or, in an eigenstate basis of the electric field operator, In U (1) QLM for general occupation N , the link operator and the electric field fulfil with a unitary operator or parallel transporter of a U (1) gauge model. • The other extreme limit is = N 1: In this case there is only one rishon per link and the dimension of the gauge invariant Hilbert space around every vertex is three, having one empty and two occupied modes on the odd vertices and two empty and one occupied modes on the even vertices. U(2) quantum link model-the generators of the SU (2) gauge transformations fulfill the usual algebra The gauge invariant subspace corresponds to a singlet of this operator, i.e.
The g a 2 and g na 2 terms describe the abelian and non-abelian electric field energy contributions respectively, m represents the staggered mass and t the interaction between matter and gauge fields. The non-abelian part of the gauge selection rule requires that the matter and rishon particles (both spin 1 2 ) on a vertex form a color singlet, therefore they must be an even number. Still, the possible combinations of total particle number on a vertex and on a link + N x x , 1 are various. A possibility, discussed here, is the configuration that includes the uniform half-filling matter state (one matter fermion per vertex). In 1D, a simple way to achieve this is by setting = N 1, and . The local gauge invariant basis is four dimensional: , and ϕ〉 | is the doubly-occupied site, with the two spin- 1 2 particles forming a spin singlet. Later we will consider the former scenario as an example, and also discuss a slightly more complex configuration (with fermionic rishons, = N 2 rishons per link, and − − 3 ( 1) x particles on vertex x). Quantum dimer and spin-ice models-in these models the matter field is fixed, and constitutes no quantum degree of freedom. The dynamics involves only gauge degrees of freedom, which are encoded in spins (hereafter we use spins- . The spin-ice and dimer model share the same gauge symmetry generator, which reads , , x y x y however, in the two cases a different symmetry sector (irreducible subspace) is selected. The QLM prescription splits the spin-1 2 in a pair of rishons, which are spinless fermions in both cases: . This model has been introduced to describe the presence of a Cooper pair or valence bond formed by a pair of electrons on the nearest neighbor vertices (the dimers). The gauge constraint arises from the fact that every electron can only pair with one of the neighbor electrons, which results in the local conservation σ + μ + . This gauge constraint reduces the Hilbert space from 2 4 to just 4 valid configurations around a vertex.
The quantum spin-ice model is similar but not identical. In this case the local gauge symmetry conservation originates from a strong antiferromagnetic Ising-type interaction between every pair of spins around a vertex: x y x y Effectively this interaction projects the Hilbert space to the zero magnetization subspace , phys x y x y . The local gauge invariant space is reduced to configurations with two spins ↑ 〉 | and two spins ↓ 〉 | around a vertex, resulting in a local gauge vertex space dimension of 6 instead of 2 4 .

Matrix product formulation of the QLM constraints
In this section we embed the previous lattice gauge picture within the tensor network framework. We first sketch a general technique, based on projected entangled pairs on the links, which allows one to operatively take into account the QLM constraints defined previously, while reducing the computational space dimension and thus the complexity of the related algorithms. The idea is to exploit the Gauge constraints to reduce the local space dimension, and at the same time combine all the link constraints into simple projectors which act directly upon the reduced space and, in 1D, are conveniently written in the matrix product operator (MPO) formalism.
As we have seen in the previous examples, the gauge constraint and the link constraint in the QLM formulation result in a description of the system, composed by logical sites, that groups a vertex of the original model and the the nearest-neighbor interacting rishon sites. Therefore, we can introduce a computational vertex site that is formed by the tensor product of a matter site and the rishon sites at that vertex, of compound dimension = ψ D d d ( ) c z , where ψ d is the matter local Hilbert space dimension, z is the coordination number of the lattice, and d c is the local rishon space dimension (equal to = + d N 1 c in the abelian gauge case, larger otherwise). We show in the following that the gauge constraint can be solved by reducing the local site Hilbert space and that the remaining link constraints can be exactly written in a simple tensor structure that we can exploit to develop efficient implementations of numerical algorithms.
To be precise, we restrict the local physical space to the trivial irreducible representation subspace φ 〉 |  G . Now let P x be the projector upon the physical space related to vertex x, and 〉 j | x r an orthonormal basis for its range (which coincides with its support, since = = P P P x x x 2 † ). The subscript r indicates that we reduced the effective dimension to = d P rnk( ), since the rank of P x is always smaller than the original dimension D of the combined degrees of freedom of vertex x, so that < d D. Then we have, for a one-dimensional QLM, the generalization to any lattice and dimensionality is given by equation (6). The linear transformation of equation (15) implements the map from the original D-dimensional basis to the d-dimensional basis of gauge constrained states and has a rectangular matrix representation A with j as the row index and the combination of ψ s , − s and + s as the column index (we dropped the vertex index x for comfort of notation). Since we chose an orthonormal reduced basis it follows that  In a QLM formulation the link constraint has to be satisfied simultaneously. As previously stated, the link symmetry group is always U (1) and thus generated by a single operator per lattice link which reads = +  . Now, since all the P x act on mutually disjointed degrees of freedom for different x (and so do the A x and the x which represent the constraints combined over the whole lattice. Now we enforce the link constraint, and then we contract the space onto the gauge-reduced basis. Basically, if we start from a generic, unconstrained many-body state Ψ 〉 | we get | r is now a generic many-body state in the gauge-reduced space, and ≡ Q AQĀ¯¯r † is the link constraint projector expressed in the reduced space. Notice that Q r is again a projector, . Moreover it is possible to write Q r as follows: Equation (21), diagrammatically represented in figure 3, is the MPO formulation of the projector Q r , with the common index q ℓ shared by two neighboring tensors, F ℓ [ ] and + F ℓ [ 1] and assuming = + distinct values (all integers from 0 to . This integer m is often referred to as bondlink dimension, and it has a physical relevance in tensor networks, since it relates to the entanglement properties of the state or operator described via tensor network ansatz [62]. For instance, in the DMRG, the entanglement entropy under a left-right partition of the variational many-body state is bound by m log . By construction the effective Hamiltonian expressed within the reduced space will preserve the link symmetry as it did in the original formulation.  In conclusion, to simulate the dynamics of a QLM, one can work completely in the reduced space and start the evolution in a quantum state of the form Ψ 〉 Q | r r 0 , where Q r enforces the link constraint. The gauge-symmetric reduced Hamiltonian H r will then preserve the link constraint since Moreover, it is possible to apply the projector Q r at any time during state evolution, for instance to prevent the state from violating the link constraint due to uncontrolled numerical errors. As previously mentioned, the MPO formulation for the reduced link projector Q r can be generalized for any lattice and dimensionality in a straightforward manner: what one obtains is a projected entangled pair operator (PEPO), again with the bondlink dimension bounded by

Canonical link-gauge basis
As an additional remark, we will demonstrate that by introducing a particular basis 〉 j | r for the reduced space, the picture simplifies further: in the new basis Q r reads as a diagonal operator without increasing the previous MPO bond link dimension. We start by recalling that in the original QLM picture, the gauge generators ν G x conserve the number of rishons on their related links − n x, and + n x, separately, i.e.
This means that there exists a basis ψ 〉 | j x in the space defined by 〉 ψ − + s s s | , , x which diagonalizes simultaneously all the operators appearing in equation (25). Within this set, we identify those that satisfy the gauge constraint, and select them as the reduced basis ψ 〉 . Such simplified decomposition is sketched in figure 4 (left panel). Notice that + is exactly the Schmidt rank of the operator + Q r x x , , 1 , so this decomposition is optimal in bondlink dimension m . Combining all the + Q r x x , , 1 together is straightforward now, since they are nearest-neighbor projectors diagonal in the reduced basis: doing so leads again to an MPO form of Q r like equation (21), but with simpler tensor blocks: j j [ ] as sketched in figure 4 (right panel). We know that this MPO representation is optimal in bondlink dimension m because it uses the minimal bondlink to represent faithfully the Schmidt ranks of the matrices + Q r x x , , 1 . Such representation is extremely versatile: we will exploit it, for instance to understand how QLM space dimensions (and thus computational costs) grow as a function of the total system size, in section 5.

Fast link-constrained time-evolution scheme
As mentioned, since the Hamiltonian commutes with every gauge or link symmetry in the original model, time-evolution of the QLM dynamics should theoretically preserve all the constraints. Unfortunately, in numerical frameworks systematic errors are generated which may have a dramatic and disruptive impact in the conservation of symmetries (if not addressed properly), e.g. in real-time evolution. The imaginary-time evolution does not suffer from this issue: in fact, since local gauge symmetries can not be spontaneously broken [63], convergence of the algorithm to the gauge-invariant ground state is guaranteed. However, even in this scenario addressing the gauge symmetry explicitly is computationally helpful: setting nongauge invariant states, which might be low-energy excitations, out of the variational picture can only speed-up the convergence rate to the ground state.
Moreover, the quasi-local constraints will allow us to significantly speed-up the timeevolution algorithms by performing all the linear algebraic operations in a computationally efficient block-wise fashion.

Enforcing link constraints over time
In this section we assume that we want to apply a (real or imaginary) time-evolution scheme of a nearest-neighbor, time-independent QLM Hamiltonian H onto a many-body (unnormalized) mixed state ρ: We also assume to have ρ expressed variationally in a matrix product density operator (MPDO) formulation, i.e. instead of numerically addressing ρ, we store the many-body operator X such that ρ = XX † , which always exists since ρ > 0 and we encode X as an MPO. The time evolution can be then carried out directly on X, because applying β β β ( ) e ( ) for real time, or e for imaginary time, Ht H 0 i¯0 0¯2 0 recovers exactly equation (29) via ρ = XX † [64]. Here we focus on nearest-neighbor Hamiltonians and thus it is convenient to evolve the state by time-evolved block decimation, a well-known procedure in DMRG contexts based on Suzuki-Trotter (ST) decomposition of H into odd-even site blocks and even-odd site blocks [27]. More precisely, ⎛ where p is known as the ST-order and the coefficients c t and d t are calculated via the Baker-Campbell-Hausdorff formula. To enforce the link constraint one might evolve the state via γ∑ More generally one might want to apply the link projector Q either before, after or before and after the evolution, since = Q Q2 . In this instance we chose to apply it after the evolution step. We now show that within the presented framework this is straightforward and requires no additional computational cost. We start by showing that

Link constraint computational speed-up
Here we show that the link constraint formulation allows us to gain a consistent advantage in both of the two elementary operations on the MPDO architecture required to apply [ , ] on X (remember that the many-body state is ρ = XX † ), namely: 1. the contraction and 2. the singular value decomposition (SVD)-truncated separation. These two operations between multilinear tensors are represented for the reader in figure 5. Let us first recall that the MPDO design stores the 'semi-state' X in the form The first term accounts for the cost of contracting the two X tensors together, the second term for assembling M. Exploiting the link symmetry in this procedure is achieved by considering that both physical label pairs in input ′ ″ j j ( , ) and in output + j j ( , ) x x 1 must satisfy the link constraint, i.e. the triplets ′ ″ j j q ( , , ), and ′ + j j q ( , , ) x x 1 must belong to Ω + x x , 1 for some q and ′ q . All the other pairs are identically zero both in input and output, and thus need not be considered in the computation. This remark reduces the computation as follows . (ii) SVD-truncated separation-the second operation is needed to maintain the MPDO structure. Thus, we split the Γ tensor back into two blocks Y via an SVD, so that x As usual, at every step one keeps the correlation bond link dimension m under control by discarding the smallest singular values. Since the cost of an SVD for a × s s ( )-dimensional square matrix scales like s 3 , the standard cost for this operation is ∼d b m 3 3 3 . This operation is quickened by exploiting the link constraint, observing that Γ is shaped with an internal block structure. In particular, an entire × bm bm ( )-dimensioned block Γ . Before performing the SVD, we reshuffle rows and columns of Γ blockwise (this operation clearly preserves the SVD decomposition). In particular, we reorder the rows so that + n x j ( , ) x is monotonically increasing while descending the rows, and we reorder the columns so that +  speed-up for the SVD procedure which is often the computational bottleneck of the timeevolution algorithm.

Dimension of QLM spaces: cellular automata
A fundamental issue that arises while numerically addressing a quantum many-body problem is the amount of computational resources required for the exact microscopic description scale with the total system size ℓ. This is a general problem which becomes even more relevant in the presence of symmetries: the constraints introduced by the additional integrals of motion often reduce the Hilbert dimension scaling with ℓ, up to the point where the numerical complexity might change dramatically. To understand how the one-dimensional QLM space dimension grows with the system size ℓ we work in the reduced basis, while assuming an open boundary conditions (OBC) setup so to proceed inductively by adding one site at a time, say from left to right. We consider a QLM chain of length ℓ; assume that we classified the 'physical' states, which are of the form Ψ 〉 Q l ( )| r rℓ , according to the rightmost link charge q ℓ . That is, we characterize a many-body orthogonal basis labeling the states via 〉 q k | , In all the cases we considered, the dimension growth reduction due to the link symmetry is not stringent enough to make the scaling polynomial. In fact, the scaling is still exponential α ∝ D l ( ) ℓ but the basis α is strictly smaller than the reduced local space dimension, i.e. α < d. Moreover, we found α to be even smaller than the total number of allowed local matter states 〉 ψ s | . Before showing some examples on how the cellular automata works in practice, and what insight it can provide, we wish to remind the reader that the present scheme is meant only for one-dimensional quantum link models. Indeed, higher-dimensionality lattices would require one to keep track of the intermediate charges for every open link when growing the lattice siteby-site, ultimately resulting in a more difficult treatment which can not be trivially translated into a cellular automata paradigm. This is nevertheless an interesting problem and it will constitute the focus for future research.  with support dimension χ = 5, or equivalently which requires a correlation bondlink m = 2, i.e. we have two nodes in the cellular automata q = 0 and q = 1. Regarding the automata connections we have: state 〉 |1 r connects q=1 to the left, and q = 0 to the right, so it is an arrow from q = 1 to q = 0. State 〉 |2 r is an arrow from q = 0 to q = 0, while state 〉 |3 r goes from q = 0 to q = 1. A visual representation of this cellular automaton is shown in figure 6. One can check immediately that the dimension of this Hilbert space grows with ℓ exactly as the Fibonacci sequence: in fact with φ being the golden ratio φ = + (1 5 ) 2. This tells us that for large sizes ℓ, the Hilbert dimension grows exponentially α ∝ D l ( ) ℓ , but instead of using an exponential basis the local space dimension d = 3 or the matter local dimension = ψ d 2, it is α = + ≃ (1 5 ) 2 1.618, i.e. the scaling is somewhat smoother.

Example: U (1), spinless matter fermion, multiple rishons
Again we explore the U (1) QLM scenario, but this time we fix the number of rishons per link to be N , while the matter is again a two-level quantum system. In this model there are = D N 2¯2 local states available, however only = + d N 2¯1 are allowed by the gauge constraint. We label them according to    with support dimension χ = + N 4¯1. The allowed charges q here go from 0 to N , so we have + N 1 nodes in the cellular automata. The arrows are defined as follows: odd index states and even index states connect the nodes as After a reordering of all the nodes, the cellular automata appear as sketched in figure 7. We calculated scaling of the total QLM Hilbert space dimension D l ( ) with the system size ℓ, for several different rishon number choices N . In every case considered starting from sizes of ∼ ℓ 10, D l ( ) matches an exponential scaling in ℓ. We fitted the exponential basis α N (¯) α ∝ D ℓ N ( ) (¯) ℓ in the interval ∈ ℓ [100, 1000]. A smooth, monotonic behavior of α as a function of N is observed and reported in figure 8. The calculated α values are never greater than 2, and saturate to 2 in a polynomial fashion with increasing numbers of rishons per link N , i.e. when approaching the Wilson limit (as shown in figure 8, inset). This study reveals that in Figure 7. Sketch of the cellular automata for the U (1) QLM with N rishons on the link, where + N 1 is the number of nodes in the picture. The rightmost node, the only one with a self-pointing arrow, is the one related to charge = ⌊ ⌋ q N 2 .   . The cellular automata for this setup is shown in figure 9, middle panel. Using the automata mechanism, we numerically calculated the effective Hilbert space dimensions D l ( ) for system sizes up to ∼ ℓ 850. Once again, an asymptotically exponential scaling α ∝ D l ( ) ℓ is detected: In figure 9, right panel, we show how the exponential curve (cyan line) fits the data points, (which have been enlarged on purpose not to be hidden by the fit curve). The exponential basis we estimated from the fit is α ≃ 2.2470. ( ) as a function of the 1D chain length ℓ, in the U (2) double rishon case, evaluated numerically (black dots). The cyan line is an exponential fit, revealing the exponential basis α ≃ 2.246 9796.

Conclusions
In this work we have merged the quantum link formalism with the tensor network framework and showed that in combination they allow one to efficiently describe both equilibrium and outof-equilibrium properties of lattice gauge theories in the Hamiltonian formulation. We showed how to efficiently combine gauge constraints and link constraints in matrix product operator formalism in 1D, which can be easily generalized to a projected entangled pair formalism in higher dimensions. This paradigm is instrumental to merge time-evolution schemes, native to tensor network architectures, with gauge-invariance constraints ultimately leading to a symmetry protected dynamics algorithm. Moreover, the local symmetries can be furthermore exploited to obtain a substantial enhancement in the algorithm performance. Finally, we adopted the tensor network picture and developed a cellular automata formalism to compute the scaling of the gauge-invariant subspace of the quantum link models, and thus the effective complexity of the model. This analysis might be useful to estimate the computational complexity of a simulation of a given model and to guide the search for mappings from the original model to simplified ones.
The framework introduced here will pave the way to the study of extremely interesting lattice gauge problems, ranging from high energy physics in low dimensions up to topological condensed matter models. Indeed, global symmetries (e.g. conserved particle numbers or total magnetization) might be combined with this approach to achieve even higher performances.