# Lattice gauge tensor networks

Pietro. Silvi, Enrique. Rico Ortega, Tommaso. Calarco, Simone. Montangero

► To cite this version:

Pietro. Silvi, Enrique. Rico Ortega, Tommaso. Calarco, Simone. Montangero. Lattice gauge tensor networks. New Journal of Physics, 2014, 16 (10) (Article number: 103015), 10.1088/1367-2630/16/10/103015. [DOI: 10.1088/1367-2630/16/10/103015](https://doi.org/10.1088/1367-2630/16/10/103015)

HAL Id: hal-04918518

https://hal.science/hal-04918518v1

Submitted on 29 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
PAPER • OPEN ACCESS

Lattice gauge tensor networks

To cite this article: Pietro Silvi et al 2014 *New J. Phys.* **16** 103015

View the [article online](#) for updates and enhancements.

You may also like

- Tensor renormalization group for fermions
  Shinichiro Akiyama, Yannick Meurice and
  Ryo Sakai

- Group field theory and tensor networks;
  towards a Ryu–Takayanagi formula in full
  quantum gravity
  Goffredo Chirco, Daniele Oriti and Mingyi
  Zhang

- Coarse graining methods for spin net and
  spin foam models
  Bianca Dittrich, Frank C Eckert and
  Mercedes Martin-Benito
# Lattice gauge tensor networks

Pietro Silvi<fnref n="1" />, Enrique Rico<fnref n="1" />,<fnref n="2" />, Tommaso Calarco<fnref n="1" />,<fnref n="3" /> and Simone Montangero<fnref n="1" />,<fnref n="3" />

E-mail: pietro.silvi@uni-ulm.de

Received 5 May 2014, revised 14 August 2014

Accepted for publication 28 August 2014

Published 9 October 2014

New Journal of Physics 16 (2014) 103015

doi:10.1088/1367-2630/16/10/103015

## Abstract

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.

Keywords: lattice gauge invariant field theory, tensor networks, quantum many-body physics, classical simulation of quantum models

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Institut für Quanteninformationsverarbeitung, Universität Ulm, D-89069 Ulm, Germany</footnote>

<footnote n="2">²IPCMS (UMR 7504) and ISIS (UMR 7006), Université de Strasbourg and CNRS, Strasbourg, France</footnote>

<footnote n="3">³Center for Integrated Quantum Science and Technology (IQST), Universities of Ulm/Stuttgart and MPI for Solid State Research, Germany</footnote>
# 1. 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–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–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–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–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–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–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 increasing 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.

## 2. 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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) The commutation relations $[H, G_x^\nu] = 0$ guarantee that the gauge invariant subspace, i.e. the trivial irreducible representation subspace for every lattice gauge subgroup, is dynamically decoupled from the rest of the Hilbert space. (b) The nontrivial support of every lattice gauge generator is a single matter field site $\psi_x$ and all the gauge field links $U_{x,x+\mu_k}$ connected to it. (c) Typical coupling Hamiltonian terms involve two matter sites $\psi_x$ and $\psi_{x+\mu_k}$ and the gauge boson connecting them $U_{x,x+\mu_k}$. (d) In the QLM formulation, the gauge boson is split into a pair of rishons, linked together by a $U(1)$ symmetry constraint.</figcaption></figure>

quantum link model language of lattice gauge theories follows without additional loss of generality [50–52]. We define the gauge invariant model of interests by defining three elements:

* The **local degrees of freedom** [$\psi_x^a$, $U_{x,x+\mu_k}^{ab}$, $E_{x,x+\mu_k}$]—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^\nu$]—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.

## 2.1. 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 $\psi_x$ are located on the vertices of the lattice $x$. They are usually fermionic fields that describe the ‘quarks’ of the model, $\{\psi_x^a, \psi_y^\dagger\} = \delta_{x,y}$. They can also be bosonic fields describing, for instance the Higgs field. In non-abelian models, fermions $\psi_x^a$ carry color degrees of freedom $a$. For example, in $U(2)$ or $SU(2)$ models $a \in \{\uparrow, \downarrow\}$, in $U(3)$ or $SU(3)$ models $a \in \{b, g, r\}$

* Gauge field $U_{x, x+\mu_k}^{ab}$ exist on the links of the lattice $\langle x, x+\mu_k \rangle$. 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: $U_{x, x+\mu_k}^{ab} = c_{x+\mu_k, -\mu_k}^{a\dagger} c_{x, +\mu_k}^b$, as sketched in figure 1, panel d.
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 $\langle x, + \mu_x \rangle$ and $\langle x + \mu_x, - \mu_x \rangle$. 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 $[c_{x,\mu_x}^b, c_{y,\mu'_y}^{a\dagger}]_{\pm} = \delta_{a,b} \delta_{x,y} \delta_{\mu_x,\mu'_y}$. The statistics of the rishon fields $c_{x, \mu_x}^a$ are completely arbitrary and do not change the statistics of the original gauge bosons $U_{x, x + \mu_x}^{ab}$, since the rishon operators $c_{x, \mu_x}^a$ always appear in pairs related to the same link. Notice that due to the hard-core nature of fermionic statistics, fermionic rishons pose limits to the maximal number of rishons per link. The total number of rishons $N_{x,x+\mu_x} = n_{x+\mu_x, -\mu_x} + n_{x, +\mu_x} = \sum_a (c_{x+\mu_x, -\mu_x}^{a\dagger} c_{x+\mu_x, -\mu_x}^a + c_{x, +\mu_x}^{a\dagger} c_{x, +\mu_x}^a)$ on every link is a conserved quantity. This is due to the fact that the rishon degrees of freedom $c_{x, \mu_x}^a$ appear both in the gauge symmetry operators $G_x^\nu$ and in the Hamiltonian $H$ only via $U_{x, x + \mu_x}^{ab}$ and by construction $[N_{x,x+\mu_x}, U_{y, y + \mu_y}^{ab}] = 0$, from which it follows that $[N_{x,x+\mu_x}, G_y^\nu] = [N_{x,x+\mu_x}, H] = 0$. 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 $N_{x,x+\mu_x}$ (regardless of the symmetry group generated by $G_x^\nu$ which may as well be non-abelian). There are different representations of the same symmetry depending on the number of rishons per link $\tilde{N}$ one selects. In any case, we restrict the Hilbert space to the ‘physical’ states $|\varphi_{\text{phys}}\rangle$ which satisfy $N_{x,x+\mu_x} |\varphi_{\text{phys}}\rangle = |\varphi_{\text{phys}}\rangle \tilde{N}_{x,x+\mu_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^\nu$ instead (see next paragraph). With a little abuse of notation, in cases where the total number of rishons on a link $\tilde{N}_{x,x+\mu_x}$ is independent of the link itself (i.e. homogeneous and isotropic QLM) we will sometimes omit the link label subscript.

## 2.2. Local generators of the gauge symmetry, and gauge constraint (Gauss' law)

The gauge symmetry is defined via the set of its generators $G_x^\nu$: they all commute with the Hamiltonian $[H, G_x^\nu] = 0$, and have localized support. To properly characterize the generators $G_x^\nu$, 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. $E_{x,x+\mu_x} = \frac{1}{2}(n_{x+\mu_x, -\mu_x} - n_{x, +\mu_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,

$$ \tilde{U}_{x, x + \mu_x}^{ab} = e^{i \theta E_{x,x+\mu_x}} U_{x, x + \mu_x}^{ab} e^{-i \theta E_{x,x+\mu_x}} = e^{i \theta} U_{x, x + \mu_x}^{ab}, \quad (1) $$

or infinitesimally $[E_{x,x+\mu_x}, U_{x, x + \mu_x}^{ab}] = U_{x, x + \mu_x}^{ab}$.

* **Non-abelian $SU(N)$:** in this scenario, the corresponding non-abelian version of the electric field has a left component $L_{x, x + \mu_x}^\nu = \sum_{ab} c_{x, + \mu_x}^{a\dagger} \frac{\lambda_{ab}^\nu}{2} c_{x, + \mu_x}^b$ and a right component
$R_{x, x + \mu_x}^\nu = \sum_{ab} c_{x + \mu_x, -\mu_x}^{a\dagger} \frac{\lambda_{ab}^\nu}{2} c_{x + \mu_x, -\mu_x}^b$ operator, depending on whether their action changes the bosonic gauge field $U_{x, x + \mu_x}^{ab}$ with a unitary $\Omega_{ak}$ acting on the left or on the right of the field as

$$ \begin{aligned} \tilde{U}_{x, x + \mu_x}^{ab} &= e^{i \sum_{\nu} \theta^{\nu} L_{x, x + \mu_x}^{\nu}} U_{x, x + \mu_x}^{ab} e^{-i \sum_{\nu} \theta^{\nu} L_{x, x + \mu_x}^{\nu}} = \sum_k \Omega_{ak} U_{x, x + \mu_x}^{kb} \\ \tilde{U}_{x, x + \mu_x}^{ab} &= e^{i \sum_{\nu} \theta^{\nu} R_{x, x + \mu_x}^{\nu}} U_{x, x + \mu_x}^{ab} e^{-i \sum_{\nu} \theta^{\nu} R_{x, x + \mu_x}^{\nu}} = \sum_k U_{x, x + \mu_x}^{ak} \Omega_{bk}^* \end{aligned} \quad (2) $$

or infinitesimally $[L_{x, x + \mu_x}^\nu, U_{x, x + \mu_x}^{ab}] = -\sum_k \lambda_{ak}^\nu U_{x, x + \mu_x}^{kb}$ and $[R_{x, x + \mu_x}^\nu, U_{x, x + \mu_x}^{ab}] = \sum_k U_{x, x + \mu_x}^{ak} \lambda_{kb}^\nu$, where $\lambda^\nu$ are the Hermitian generators of SU(N) which obey $[\lambda^\mu, \lambda^\nu] = 2if_{\mu\nu\omega}\lambda^\omega$, with $f_{\mu\nu\omega}$ the structure constants of the SU(N) algebra and $\text{Tr}(\lambda^\mu\lambda^\nu) = 2\delta^{\mu\nu}$.

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

$$ \begin{aligned} G_x &= \psi_x^\dagger \psi_x + \sum_{\mu_x} (E_{x-\mu_x,x} - E_{x,x+\mu_x}) \\ &= \psi_x^\dagger \psi_x + \sum_{\mu_x} (n_{x,-\mu_x} + n_{x,+\mu_x}) - \frac{\bar{N}_{x-\mu_x,x} + \bar{N}_{x,x+\mu_x}}{2}, \end{aligned} \quad (3) $$

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 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

$$ \left[ \psi_x^\dagger \psi_x + \sum_{\mu_x} (n_{x,-\mu_x} + n_{x,+\mu_x}) \right] |\varphi_{\text{phys}}\rangle = |\varphi_{\text{phys}}\rangle \times \text{constant}, \quad (4) $$

which is equivalent to Gauss' law, and the physical states $|\varphi_{\text{phys}}\rangle$ are those which satisfy it.

* The generators for the non-abelian SU(N) part of the gauge transformations fulfill the usual algebra $[G_x^\mu, G_y^\nu] = i\epsilon^{\mu\nu\omega}G_x^\omega \delta_{x,y}$ with

$$ \begin{aligned} G_x^\nu &= \sum_{ab} \psi_x^{a\dagger} \frac{\lambda_{ab}^\nu}{2} \psi_x^b + \sum_{\mu_x} (R_x^\nu - \mu_x, x + L_x^\nu, x + \mu_x) \\ &= \sum_{ab} \left[ \psi_x^{a\dagger} \frac{\lambda_{ab}^\nu}{2} \psi_x^b + \sum_{\mu_x} \left( c_{x-\mu_x, x}^{a\dagger} \frac{\lambda_{ab}^\nu}{2} c_{x-\mu_x, x}^b + c_{x+\mu_x, x}^{a\dagger} \frac{\lambda_{ab}^\nu}{2} c_{x+\mu_x, x}^b \right) \right] \end{aligned} \quad (5) $$

Again, in QLM formulation $G_x^\nu$ 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^\nu$, i.e. the singlet subspace: $G_x^\nu |\varphi_{\text{phys}}\rangle = 0$. 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^\nu$ 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 $\psi_x^a$ and $z$ gauge field sites $U_{x,x+\mu_x}^{ab}$ (or rishon sites $c_{x,\mu_x}^a$ 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^\nu$ 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

$$|j_x\rangle_r = \sum_{s_{\psi}, \vec{s}_{\mu_x}} A_{s_{\psi}, \vec{s}_{\mu_x}}^{[x]} |s_{\psi}, \vec{s}_{\mu_x}\rangle_x, \qquad (6)$$

once a canonical basis for the matter field $|s_\psi\rangle_x$ and one for each rishon field $|s_{\mu_x}\rangle_x$, have been chosen. This mapping defines a 'reduced' local basis $|j_x\rangle_r$ which spans exactly and solely the local gauge-invariant subspace. In the next sections, the set of states $|j_x\rangle_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.

## 2.3. Gauge invariant dynamics

The last element that has to be 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. $[H, G_x^\nu] = [H, \mathcal{N}_{x,x+\mu_x}] = 0$. 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:

$$ \begin{aligned} H_{\text{pure}} = H_{\text{electric}} + H_{\text{magn}} & = \sum_{x, \mu_x} \left\{ g_{\text{abel}}^2 (E_{x,x+\mu_x})^2 + g_{\text{non-ab}}^2 - \sum_{\nu} \left[ (L_{x,x+\mu_x}^\nu)^2 + (R_{x,x+\mu_x}^\nu)^2 \right] \right\} \\ & \quad - \frac{1}{g_{\text{magn}}^2} \sum_{x, \mu_x, \mu_y} \left[ \text{Tr} \left( U_{x,x+\mu_x} U_{x+\mu_x,x+\mu_x+\mu_y} U_{x+\mu_x+\mu_y,x+\mu_y} U_{x+\mu_y,x} \right) + \text{h.c.} \right]. \end{aligned} \qquad (7) $$

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_{\text{abel}}^2$, $g_{\text{non-ab}}^2$ and $g_{\text{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 $g_{\text{magn}}^2 = 4ag^2$, and $g_{\text{abel}}^2 = g_{\text{non-ab}}^2 = \frac{g^2}{2a}$, 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 $g_{\text{abel}} \approx g_{\text{non-ab}} \gg 1/g_{\text{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,

$$H_{\text{coup}} = \sum_{x, \mu_x} J_{x, \mu_x} (\psi_x^\dagger U_{x,x+\mu_x} \psi_{x+\mu_x} + \text{h. c.}) + \sum_x m_x \psi_x^\dagger \psi_x \quad (8)$$

where we have defined site dependence hopping constants $J_{x,\mu_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.

## 2.4. 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

$$H = J \sum_x (\psi_x^\dagger U_{x,x+1} \psi_{x+1} + \text{h.c.}) + g^2 \sum_x (E_{x,x+1})^2 + m \sum_x (-1)^x \psi_x^\dagger \psi_x \quad (9)$$

where the last term is a staggered chemical potential profile for the matter field which is a spinless fermion field {$\psi_x, \psi_y^\dagger$} = $\delta_{x,y}$. 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 $U_{x,x+1} = c_{x,+}c_{x+1, -}^\dagger$, which are bosonic in nature [$c_{x,a}, c_{y, b}^\dagger$] = $\delta_{x,y}\delta_{a,b}$.

The two independent local symmetries in this $U(1)$ QLM are:

(i) Constant number of rishons per link: $\mathcal{N}_{x,x+1}|\varphi_{\text{phys}}\rangle = (n_{x+1,-} + n_{x,+})|\varphi_{\text{phys}}\rangle = \bar{N} |\varphi_{\text{phys}}\rangle$

(ii) Gauss’ law on every vertex: $(\psi_x^\dagger \psi_x + n_{x,-} + n_{x,+})|\varphi_{\text{phys}}\rangle = |\varphi_{\text{phys}}\rangle(\bar{N} - \frac{1+(-1)^x}{2})$

The factor $(\bar{N} - \frac{1+(-1)^x}{2})$ appears because we introduce $\psi_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 $\bar{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 $N_{x,x+1} = \bar{N}$ and the electric flux $E_{x,x+1} = \frac{n_{x+1,-} - n_{x,+}}{2}$, i.e., $|n_+, n_-\rangle = |\bar{N}, E\rangle$, where we have omitted the labels of the link $\langle x, x+1\rangle$:

* $\bar{N} \gg 1$ (Wilson limit) [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 $\vartheta$ that fulfil the usual commutation relation of position and momentum $[E, \vartheta] = i$. Then by defining the link operator $U = e^{-i\vartheta}$, it is straightforward to check that $[U, U^\dagger] = 0$, $[E, U] = U$ or, in an eigenstate basis of the electric field operator, $U|E\rangle = |E+1\rangle$. In $U(1)$ QLM for general occupation $\bar{N}$, the link operator and the electric field fulfil $U|\bar{N}, E\rangle = \sqrt{\frac{\bar{N}}{2}(\frac{\bar{N}}{2}+1)} - E(E+1)|\bar{N}, E+1\rangle$ and $[U, U^\dagger] = E$. In the limit $\bar{N} \gg E$,

$$ \frac{1}{\sqrt{\frac{\bar{N}}{2}\left(\frac{\bar{N}}{2}+1\right)}} U |\bar{N}, E\rangle \rightarrow |\bar{N}, E+1\rangle; \quad \frac{1}{\frac{\bar{N}}{2}\left(\frac{\bar{N}}{2}+1\right)} [U, U^\dagger] \rightarrow 0; \\ \frac{1}{\sqrt{\frac{\bar{N}}{2}\left(\frac{\bar{N}}{2}+1\right)}} [E, U] = \frac{1}{\sqrt{\frac{\bar{N}}{2}\left(\frac{\bar{N}}{2}+1\right)}} U \qquad (10) $$

which is the usual definition of the Wilson type lattice theories if we identify $\frac{1}{\sqrt{\frac{\bar{N}}{2}(\frac{\bar{N}}{2}+1)}} U$ with a unitary operator or parallel transporter of a $U(1)$ gauge model.

* The other extreme limit is $\bar{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 $[G_x^\mu, G_y^\nu] = i\epsilon^{\mu\nu\omega}G_x^\omega\delta_{x,y}$ with

$$ G_x^\nu = R_x^\nu + \sum_{a,b} \left( \psi_x^{a\dagger} \frac{\sigma_{ab}^\nu}{2} \psi_x^b \right) + L_x^\nu \\ = \sum_{a,b} \left( c_{x, -}^{a\dagger} \frac{\sigma_{ab}^\nu}{2} c_{x, -}^b + \psi_x^{a\dagger} \frac{\sigma_{ab}^\nu}{2} \psi_x^b + c_{x, +}^{a\dagger} \frac{\sigma_{ab}^\nu}{2} c_{x, +}^b \right) \qquad (11) $$

The gauge invariant subspace corresponds to a singlet of this operator, i.e. $G_x^\nu|\varphi_{\text{phys}}\rangle = 0$. A $U(2)$ gauge invariant Hamiltonian can be written as

$$ H = \frac{1}{2} \sum_x \left( g_a^2 E_x^2 + g_{na}^2 \sum_\nu \left[ (R_x^\nu)^2 + (L_x^\nu)^2 \right] \right) + m \sum_{x,a} (-1)^x \psi_x^{a\dagger} \psi_x^a \\ + t \sum_{x,a,b} \left[ \psi_x^{a\dagger} U_{x,x+1}^{ab} \psi_x^b + \text{h. c.} \right] \qquad (12) $$

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 $\frac{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 $n_{x,\psi} + n_{x,-} + n_{x,+}$ and on a link $\bar{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 $\bar{N} = 1$, and $n_{x,\psi} + n_{x,-} + n_{x,+} = 2$. The local gauge invariant basis is four dimensional: $\{|\uparrow, \downarrow, 0\rangle, |\uparrow, 0, \downarrow\rangle, |0, \uparrow, \downarrow\rangle, |0, \phi, 0\rangle\}$, where $|\uparrow, \downarrow\rangle \equiv \frac{1}{\sqrt{2}}(|\uparrow, \downarrow\rangle - |\downarrow, \uparrow\rangle)$, and $|\phi\rangle$ is the doubly-occupied site, with the two spin-$\frac{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, $\bar{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-$\frac{1}{2}$ for simplicity) living on the links of a square lattice. The gauge symmetry generators are built upon one component of the Pauli matrices vector, say the third one $\sigma_{x, x+\mu}^z$. The spin-ice and dimer model share the same gauge symmetry generator, which reads

$$G_x = \sigma_{x, x+\mu_x}^z + \sigma_{x, x+\mu_y}^z + \sigma_{x-\mu_x, x}^z + \sigma_{x-\mu_y, x}^z, \qquad (13)$$

however, in the two cases a different symmetry sector (irreducible subspace) is selected. The QLM prescription splits the spin-$\frac{1}{2}$ in a pair of rishons, which are spinless fermions in both cases: we thus rewrite $\sigma_{x, x+\mu}^z = \frac{1}{2}(n_{x+\mu, -\mu} - n_{x, \mu})$, obviously yielding $[\sigma_{x, x+\mu}^z, N_{x,x+\mu}] = 0$. The link sector selected ($\bar{N} = 1$, i.e. $N|\varphi_{\text{phys}}\rangle = |\varphi_{\text{phys}}\rangle$), recovers exactly a two-level system on every link, as shown in figure 2.

In the quantum dimer model, the lattice is covered with dimer configurations on the links of the lattice. The dimer Hilbert space is characterized by the state $|\rangle$ (on the vertical links) or $|\rangle$ (on the horizontal links) with $\sigma_{x, x+\mu}^z = 1/2$ if the link is occupied, otherwise the state is $|\rangle$ on the link and $\sigma_{x, x+\mu}^z = -1/2$. 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 $(\sigma_{x, x+\mu_x}^z + \sigma_{x, x+\mu_y}^z + \sigma_{x-\mu_x, x}^z + \sigma_{x-\mu_y, x}^z)|\varphi_{\text{phys}}\rangle = -|\varphi_{\text{phys}}\rangle$. 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:

$$H_{\text{Ising}} = (\sigma_{x,x+\mu_x}^z + \sigma_{x,x+\mu_y}^z + \sigma_{x-\mu_x,x}^z + \sigma_{x-\mu_y,x}^z)^2. \qquad (14)$$

Effectively this interaction projects the Hilbert space to the zero magnetization subspace $G_x|\varphi_{\text{phys}}\rangle = (\sigma_{x, x+\mu_x}^z + \sigma_{x, x+\mu_y}^z + \sigma_{x-\mu_x, x}^z + \sigma_{x-\mu_y, x}^z)|\varphi_{\text{phys}}\rangle = 0$. The local gauge invariant space is reduced to configurations with two spins $|\uparrow\rangle$ and two spins $|\downarrow\rangle$ around a vertex, resulting in a local gauge vertex space dimension of 6 instead of $2^4$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Gauge generator supports in (a) the standard formulation and in (b) the quantum link formulation of LGTs of frustrated spin systems. Blue circles represent sites of the lattice, orange ones the link degrees of freedom, i.e. the spins. (a) The red square on the lattice highlights the degrees of freedom on which the local gauge invariant generator acts: a site and the connected links. (b) The red diamond shows the degrees of freedom on which the local gauge invariant generator acts for a QLM. The original link degree of freedom is split into two rishons, which are modeled by spinless fermions in this context.</figcaption></figure>

**3. 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_{\psi}(d_c)^z$, where $d_{\psi}$ 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_c = N + 1$ 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 $|\varphi_{\text{phys}}\rangle_x$ of the local gauge symmetry group at vertex $x$, identified by $G_x^\nu|\varphi_{\text{phys}}\rangle_x = 0$, with $G_x^\nu$ the group symmetry generators we defined in the previous section. Since gauge
symmetries on different vertices commute, i.e. $[G_x^\nu, G_{x'}^{\nu'}] = 0$, we can enforce the gauge requirement simultaneously on all vertices $x$. Typical examples are the SU(2) or SO(3) gauge group cases, where the restricted states $|\phi_{\text{phys}}\rangle$ are those vertex states which behave like a spin-0 under $G^\nu$. Now let $P_x$ be the projector upon the physical space related to vertex $x$, and $|j_x\rangle_r$ an orthonormal basis for its range (which coincides with its support, since $P_x = P_x^2 = P_x^\dagger$). The subscript $r$ indicates that we reduced the effective dimension to $d = \text{rk}(P)$, 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,

$$|j_x\rangle_r = \sum_{s_\psi, s_-, s_+} A_{s_-, s_\psi, s_+}^{[x]j} |s_-, s_\psi, s_+\rangle_x, \quad (15)$$

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_\psi, 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

$$\begin{aligned} \delta_{j,j'} &= \langle j | j' \rangle_r = \sum_{s_\psi, s_-, s_+} \sum_{s'_\psi, s'_-, s'_+} A_{s_-, s_\psi, s_+}^{[x]j} A_{s'_\psi, s'_\psi, s'_+}^{*[x]j'} \langle s'_-, s'_\psi, s'_+ | s_-, s_\psi, s_+ \rangle_x \\ &= \sum_{s_\psi, s_-, s_+} A_{s_-, s_\psi, s_+}^{[x]j} A_{s'_\psi, s'_\psi, s'_+}^{*[x]j'} \end{aligned} \quad (16)$$

or, in matrix representation, $AA^\dagger = \mathbb{1}_r$, i.e. $A^\dagger$ is an isometry. Similarly $A^\dagger A = P$, and thus $AP = A$. The reduced basis $|j_x\rangle_r$ defines the local computational basis for any type of simulation on QLMs, since it generates the full set of states fulfilling the gauge constraint.

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 $\mathcal{N}_{x,x+1} = n_{x,+} + n_{x+1,-}$. Here the operator $n_{x,\pm} = \sum_a c_{x,\pm}^{a\dagger} c_{x,\pm}^a$ counts the total number of rishons in the mode $\langle x, \pm \rangle$, disregarding their color $a$. By construction, the link group commutes with the Hamiltonian, i.e. $[H, \mathcal{N}_{x,x+1}] = 0$, as well as with the gauge group, i.e. $[\mathcal{N}_{x,x+1}, G_x^\nu] = 0$. The link constraint requires that the number of rishons on the link $\langle x, x+1 \rangle$ is fixed to an integer number $\bar{N}_{x,x+1}$, which means $\mathcal{N}_{x,x+1}|\phi_{\text{phys}}\rangle = \bar{N}_{x,x+1}|\phi_{\text{phys}}\rangle$.

The link constraint can be implemented by applying a projector $Q_{x,x+1} = Q_{x, x+1}^2 = Q_{x, x+1}^\dagger$, which is diagonal as every chosen rishon basis state $|s_\mu\rangle$ has a well defined occupation number $n_{x,\pm}|s_\pm\rangle_{x,\pm} = |s_\pm\rangle_{x,\pm}\bar{n}_{x,\pm}(s)$. In this case it reads

$$Q_{x,x+1} = \sum_{s_-, s_+, q} B_{s_+}^{[x]} C_{q, s_-}^{[x+1]} |s_+\rangle \langle s_+|_x \otimes |s_-\rangle \langle s_-|_{x+1}, \quad (17)$$

where we split $Q_{x,x+1}$ according to its left-right Schmidt rank, resulting in

$$B_{s_+}^{[x]} = \delta_{\langle s_+ | n_{x,+} | s_+ \rangle, q} \quad \text{and} \quad C_{q, s_-}^{[x+1]} = \delta_{\bar{N}_{x,x+1}-q, \langle s_- | n_{x+1,-} | s_- \rangle}. \quad (18)$$

Of course, the fact that $[\mathcal{N}_{x,x+1}, G_x^\nu] = 0$ also implies that $[Q_{x,x+1}, P_{x'}] = 0$. Now, since all the $P_x$ act on mutually disjointed degrees of freedom for different $x$ (and so do the $A_x$ and the $Q_{x,x+1}$) we can define
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Tensor network graphical diagram, representing the MPO formulation of the combined link constraint projector in the reduced basis space $\tilde{Q}_r$. This picture corresponds to equations (21) and (22).</figcaption></figure>

$$ \bar{P} = \otimes_{x=1}^{L} P_x, \quad \bar{A} = \otimes_{x=1}^{L} A_x \quad \text{and} \quad \bar{Q} = \otimes_{x=1}^{L-1} Q_{x,x+1}, \qquad (19) $$

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 $|\Psi\rangle$ we get

$$ \bar{A}\bar{Q}|\Psi\rangle = \bar{A}\bar{P}\bar{Q}|\Psi\rangle = \bar{A}\bar{Q}\bar{P}|\Psi\rangle = \bar{A}\bar{Q}\bar{A}^{\dagger}\bar{A}|\Psi\rangle = \bar{A}\bar{Q}\bar{A}^{\dagger}|\Psi\rangle_r = \bar{Q}_r|\Psi\rangle_r, \quad (20) $$

where $|\Psi\rangle_r$ is now a generic many-body state in the gauge-reduced space, and $\bar{Q}_r \equiv \bar{A}\bar{Q}\bar{A}^{\dagger}$ is the link constraint projector expressed in the reduced space. Notice that $\bar{Q}_r$ is again a projector, since $\bar{Q}_r^2 = \bar{A}\bar{Q}\bar{A}^{\dagger}\bar{A}\bar{Q}\bar{A}^{\dagger} = \bar{A}\bar{Q}^2\bar{A}^{\dagger} = \bar{A}\bar{Q}\bar{A}^{\dagger} = \bar{Q}_r$. Moreover it is possible to write $\bar{Q}_r$ as follows:

$$ \bar{Q}_r = \sum_{j_1 \dots j_L} \sum_{j_1' \dots j_L'} F_{j_1, j_1'}^{[1]q_1} F_{j_2, j_2'}^{[2]q_1 q_2} F_{j_3, j_3'}^{[3]q_2 q_3} \dots F_{j_L, j_L'}^{[L]q_{L-1}} |j_1 \dots j_L\rangle \langle j_1' \dots j_L'|_r \quad (21) $$

where

$$ F_{j_x, j_x'}^{[x]q_{x-1}, q_x} = \sum_{s_{-}, s_{\psi}, s_{+}} A_{s_{-}, s_{\psi}, s_{+}}^{[x]j_x} C_{q_{x-1}, s_{-}, B_{s_{+}, q_x}}^{[x]} A_{s_{-}, s_{\psi}, s_{+}}^{*[x]j_x'} \quad (22) $$

Equation (21), diagrammatically represented in figure 3, is the MPO formulation of the projector $\bar{Q}_r$, with the common index $q_\ell$ shared by two neighboring tensors, $F^{[\ell]}$ and $F^{[\ell+1]}$ and assuming $m = \bar{N}_{x,x+1} + 1$ distinct values (all integers from 0 to $\bar{N}_{x,x+1}$). 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 $\log m$.

By construction the effective Hamiltonian expressed within the reduced space will preserve the link symmetry as it did in the original formulation. In fact, let $H_r = \bar{A}H\bar{A}^{\dagger}$ be the reduced Hamiltonian, then it holds that

$$ [H_r, \bar{Q}_r] = [\bar{A}H\bar{A}^{\dagger}, \bar{A}\bar{Q}\bar{A}^{\dagger}] = \bar{A}[H, \bar{Q}]\bar{A}^{\dagger} = 0. \qquad (23) $$
<figure><img src="image_6.png" /><figcaption>Figure 4. Tensor network graphical diagram of the $\bar{Q}_r$ in the canonical link-gauge basis. Left: the diagonal projector $Q_{r,x,x+1}$ decomposed according to equation (27). Right: simplified MPO representation of the combined link constraint in the reduced space $\bar{Q}_r$.</figcaption></figure>

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 $\bar{Q}_r|\Psi_0\rangle_r$, where $\bar{Q}_r$ enforces the link constraint. The gauge-symmetric reduced Hamiltonian $H_r$ will then preserve the link constraint since

$$|\Psi(t)\rangle_r = U_r(t) \bar{Q}_r |\Psi_0\rangle_r = U_r(t) \bar{Q}_r^2 |\Psi_0\rangle_r = \bar{Q}_r U_r(t) \bar{Q}_r |\Psi_0\rangle_r = \bar{Q}_r |\Psi(t)\rangle_r, \quad (24)$$

where $U_r(t) \equiv \exp(itH_r) = \bar{A} \exp(itH)\bar{A}^\dagger$. Moreover, it is possible to apply the projector $\bar{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 $\bar{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 $\bar{N}_{x,x+\mu_x} + 1$.

### 3.1. Canonical link-gauge basis

As an additional remark, we will demonstrate that by introducing a particular basis $|j\rangle_r$ for the reduced space, the picture simplifies further: in the new basis $\bar{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^\nu$ conserve the number of rishons on their related links $n_{x,-}$ and $n_{x,+}$ separately, i.e.

$$[n_{x,-}, G_{x'}^{\nu}] = [n_{x,+}, G_{x'}^{\nu}] = 0. \qquad (25)$$

This means that there exists a basis $|\psi_j\rangle_x$ in the space defined by $|s_{-}, s_{\psi}, s_{+}\rangle_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 $|\psi_j\rangle_x|_{\text{phys}} \to |j_x\rangle_r$, precisely:

$$n_{x,-}|j_x\rangle_r = |j_x\rangle_r \cdot \bar{n}_{-}(x, j) \quad \text{and} \quad n_{x,+}|j_x\rangle_r = |j_x\rangle_r \cdot \bar{n}_{+}(x, j). \qquad (26)$$

For obvious reasons, we refer to this special local basis choice as the canonical gauge-link basis. In this framework, the reduced link constraint projector reads

$$Q_{r,x,x+1}|j_x j'_{x+1}\rangle_r = |j_x j'_{x+1}\rangle_r \delta_{\bar{n}_+(x,j)+\bar{n}_-(x+1,j'), \bar{n}_x} = |j_x j'_{x+1}\rangle_r \sum_{q=0}^{\bar{N}_{x,x+1}} V_{j_x, q}^{[x]} \cdot Z_{j_x+1, q}^{[x+1]} \qquad (27)$$
where we substituted $V_{j,q}^{[x]} = \delta_{\bar{n}_+(x,j),q}$ and $Z_{q,j}^{[x+1]} = \delta_{\bar{n}_{x,x+1}-q,\bar{n}_-(x+1,j)}$. Such simplified decomposition is sketched in figure 4 (left panel). Notice that $\bar{N}_{x,x+1} + 1$ 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 $\bar{Q}_r$ like equation (21), but with simpler tensor blocks:

$$F_{j_x, j'_x}^{[x]q_{x-1}, q_x} = \delta_{j_x, j'_x} \cdot Z_{q_{x-1}, j_x}^{[x]} \cdot V_{j_x, q_x}^{[x]}, \quad (28)$$

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.

## 4. 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 non-gauge 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 time-evolution algorithms by performing all the linear algebraic operations in a computationally efficient block-wise fashion.

### 4.1. 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 $\bar{H}$ onto a many-body (unnormalized) mixed state $\rho$:

$$\rho(t_0 + t) = e^{i\bar{H}t}\rho(t)e^{-i\bar{H}t} \quad \text{for real-time, or}$$

$$\rho(\beta_0 + \beta) = e^{-\bar{H}\beta/2}\rho(\beta_0)e^{-\bar{H}\beta/2} \quad \text{for imaginary-time.} \quad (29)$$

We also assume to have $\rho$ expressed variationally in a matrix product density operator (MPDO) formulation, i.e. instead of numerically addressing $\rho$, we store the many-body operator $X$ such that $\rho = XX^\dagger$, which always exists since $\rho > 0$ and we encode $X$ as an MPO. The time evolution can be then carried out directly on $X$, because applying
$$X(t_0 + t) = e^{i\hat{H}_t} X(t_0) \quad \text{for real-time, or}$$

$$X(\beta_0 + \beta) = e^{-\hat{H}\beta/2} X(\beta_0) \quad \text{for imaginary-time,} \quad (30)$$

recovers exactly equation (29) via $\rho = XX^\dagger$ [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 $\tilde{H}$ into odd–even site blocks and even–odd site blocks [27]. More precisely,

$$\exp\left(\gamma \sum_x H_{x, x+1}^{[r]}\right) = \left(\bigotimes_x e^{c_1 \gamma H_{2x-1, 2x}^{[r]}}\right) \left(\bigotimes_x e^{d_1 \gamma H_{2x, 2x+1}^{[r]}}\right) \\ \times \left(\bigotimes_x e^{c_2 \gamma H_{2x-1, 2x}^{[r]}}\right) \dots \left(\bigotimes_x e^{d_p \gamma H_{2x, 2x+1}^{[r]}}\right) + O(\gamma^p) \quad (31)$$

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 $\tilde{Q} \exp(\gamma \sum_x \tilde{H}_{x,x+1})$. More generally one might want to apply the link projector $\tilde{Q}$ either before, after or before and after the evolution, since $\tilde{Q}^2 = \tilde{Q}$. 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 $[Q_{r,x,x+1}, H_{r,x',x'+1}] = 0$, i.e. even local Hamiltonian terms commute with the link constraints. Indeed,

$$\begin{aligned} [Q_{r,x,x+1}, H_{r,x',x'+1}] &= [\bar{A} Q_{x,x+1} \bar{A}^\dagger, \bar{A} H_{x',x'+1} \bar{A}^\dagger] \\ &= \bar{A} (Q_{x,x+1} \bar{P} H_{x',x'+1} - H_{x',x'+1} \bar{P} Q_{x,x+1}) \bar{A}^\dagger \\ &= \bar{A} (\bar{P} Q_{x,x+1} H_{x',x'+1} - H_{x',x'+1} \bar{P} Q_{x,x+1}) \bar{A}^\dagger \\ &= \bar{A} [Q_{x,x+1}, H_{x',x'+1}] \bar{A}^\dagger, \end{aligned} \quad (32)$$

as $\bar{A}\bar{P} = \bar{A}$. In the original basis $Q_{x,x+1}$ and $H_{x',x'+1}$ act on common degrees of freedom only if $x = x'$, but in this case the commutator is zero, since the local Hamiltonian term of equation (8), respects the link symmetry on the inner bond. Finally,

$$\begin{aligned} \tilde{Q} \exp\left(\gamma \sum_x \tilde{H}_{x,x+1}\right) &= \left(\bigotimes_x M_{2x-1, 2x}^{[r, 1]}\right) \left(\bigotimes_x M_{2x, 2x+1}^{[r, 1]}\right) \\ &\quad \times \left(\bigotimes_x M_{2x-1, 2x}^{[r, 2]}\right) \dots \left(\bigotimes_x M_{2x, 2x+1}^{[r, p]}\right) + O(\gamma^p) \end{aligned} \quad (33)$$

where

$$M_{2x-1, 2x}^{[r, \nu]} = Q_{r, 2x-1, 2x} \exp(c_\nu \gamma H_{2x-1, 2x}^{[r]}) , \quad M_{2x, 2x+1}^{[r, \nu]} = Q_{r, 2x, 2x+1} \exp(d_\nu \gamma H_{2x, 2x+1}^{[r]}) . \quad (34)$$

This formulation ensures that the link symmetry is always protected without increasing the computational cost. We will see now that actually one can reduce such cost by exploiting the constraint, and gain a significant speed-up of the algorithm.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Pictorial tensor-network representation of the two basic operations which can be made faster by exploiting the link constraint: (a) contraction; (b) SVD-truncated separation.</figcaption></figure>

## 4.2. 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 $M_{x, x+1}^{[r, \nu]}$ on $X$ (remember that the many-body state is $\rho = XX^\dagger$), 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

$$X = \sum_{j_1 \dots j_L = 1}^{d} \sum_{k_1 \dots k_L = 1}^{b} \sum_{w_1 \dots w_{L-1} = 1}^{m} X_{j_1, k_1}^{[1] w_1} X_{j_2, k_2}^{[2] w_1 w_2} X_{j_3, k_3}^{[3] w_2 w_3} \dots X_{j_L, k_L}^{[L] w_{L-1}} |j_1 \dots j_L\rangle_r \langle k_1 \dots k_L|_b, \quad (35)$$

where the correlation bondlink dimension $m$ and the bath bondlink dimension $b$ are both arbitrary (although $b \le m^2 d$), and determine the computational costs and the final numerical precision of the simulation. On the other hand, we now order all possible triplets of labels $(j, j', q)_{x,x+1}$ so that the corresponding state $|j_x, j'_{x+1}\rangle$, belongs to the support of $\mathcal{Q}_{r,x,x+1}$, and $q$ is the ‘intermediate charge’ of the pair i.e. $q = \bar{n}_+(x, j) = \bar{N}_x - \bar{n}_-(x+1, j')$. All these triplets are collected into the set $\Omega_{x,x+1}$, and their number is $\chi=\#\Omega_{x,x+1}$, clearly with $\chi < d^2$. After these initial remarks, we can study the two operations separately:

(i) **Contraction**—the goal of this operation is to calculate entry-wise the tensor

$$\Gamma_{w_{x+1}, k_{x+1}, j_{x+1}}^{w_{x-1}, k_x, j_x} = \sum_{w_x} \sum_{j', j''}^{m} (M_{x, x+1}^{[r, \nu]})_{j_x, j_{x+1}}^{j', j''} X_{j', k_x}^{[x] w_{x-1} w_x} X_{j'', k_{x+1}}^{[x+1] w_x w_{x+1}}, \quad (36)$$

whose cost normally scales (without considering fast matrix-multiplication schemes) as $\sim d^2 b^2 m^3 + d^4 b^2 m^2$. 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_x, j_{x+1})$ must satisfy the link constraint, i.e. the triplets $(j', j'', q)$, and $(j_x, j_{x+1}, q')$ must belong to $\Omega_{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
$$ \text{cost} \sim d^2 b^2 m^3 + d^4 b^2 m^2 \quad \text{reduced to} \quad \text{cost} \sim \chi^2 b^2 m^2 \quad (37) $$

(ii) SVD-truncated separation—the second operation is needed to maintain the MPDO structure. Thus, we split the $\Gamma$ tensor back into two blocks $Y$ via an SVD, so that

$$ \Gamma_{w_{x+1}, k_{x+1}, j_{x+1}}^{w_{x-1}, k_x, j_x} \approx \sum_{w_x}^{m} Y_{j_x, k_x}^{[x]w_{x-1}w_x} Y_{j_{x+1}, k_{x+1}}^{[x+1]w_xw_{x+1}}. \quad (38) $$

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 \times s)$-dimensional square matrix scales like $s^3$, the standard cost for this operation is $\sim d^3 b^3 m^3$. This operation is quickened by exploiting the link constraint, observing that $\Gamma$ is shaped with an internal block structure. In particular, an entire $(bm \times bm)$-dimensioned block $(\Gamma_{j_{x+1}}^{j_x, w_{x-1}, k_x} \Gamma_{j_{x+1}}^{w_{x-1}, k_x, w_x})$ is zero unless $\Omega_{x,x+1}$ contains a triplet $(j_x, j_{x+1}, q)$. Before performing the SVD, we reshuffle rows and columns of $\Gamma$ 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 $n_-(x+1, j_{x+1})$ is monotonically decreasing while moving right in the columns. Having done that, the resulting $\Gamma$ is block diagonal, with a number of blocks equal to the number of intermediate charges $q$, usually (and always up to) $\bar{N}_x + 1$. A single block, e.g. that related to intermediate charge $q$, has dimension $(bm \xi(q) \times bm \xi(q))$ where $\xi(q)$ is the number of triplets in $\Omega_{x,x+1}$ containing intermediate charge $q$. Finally, instead of performing SVD on the whole $\Gamma$ matrix, we perform separate SVDs on each distinct $q$ block. This approach reduces the computational costs as follows

$$ \text{cost} \sim d^3 b^3 m^3 \quad \text{reduced to} \quad \text{cost} \sim b^3 m^3 \sum_{q}^{\bar{N}_x} \xi^3(q). \quad (39) $$

In the best case scenario, where the blocks have all roughly the same size $\xi = d/\bar{N}_x$, the reduced cost ultimately scales as $\sim (bmd/\bar{N}_x)^3 \bar{N}_x \sim d^3 b^3 m^3 \bar{N}_x^{-2}$, resulting in a net $\bar{N}_x^{-2}$ speed-up for the SVD procedure which is often the computational bottleneck of the time-evolution algorithm.

## 5. 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 $\ell$. 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 $\ell$, 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 $\ell$ 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 $\ell$; assume that we classified the ‘physical’ states, which are of the form $\bar{Q}_r(\ell)|\Psi\rangle_{r,\ell}$ according to the rightmost link charge $q_\ell$. That is, we characterize a many-body orthogonal basis labeling the states via $|q_\ell, k\rangle_{r,\ell}$ so that
$n_{\ell,+} |q_{\ell}, k\rangle_{r,\ell} = q_{\ell} |q_{\ell}, k\rangle_{r,\ell}$, and the degeneracy label $k \in \{1, D_q(\ell)\}$ spans within this charge sector. Clearly, the total Hilbert dimension is given by $\bar{D}(\ell) = \sum_q D_q(\ell)$.

When adding a site to the previous picture, every (link-gauge) reduced basis state $|j_{\ell+1}\rangle_r$ connects only with those states $|q, k\rangle_{r,\ell}$ such that $q = \bar{N}_{\ell+1} - \bar{n}_-(\ell+1, j)$. Moreover, every state of this form will have a well defined rightmost link charge $q' = \bar{n}_+(\ell+1, j)$ and so they can be labeled again according to the rightmost link charge sectors. Such inductive steps will produce a new orthonormal complete basis of the type $|q', k'\rangle_{r,\ell+1}$. By construction the new sector dimensions read

$$D_{q'}(\ell + 1) = \sum_q \sum_j D_q(\ell) \delta_{q, \bar{N}_{\ell+1} - \bar{n}_-(\ell+1, j)} \delta_{\bar{n}_+(\ell+1, j), q'} = \sum_j D[\bar{N}_{\ell+1} - \bar{n}_-(\ell+1, j)](\ell) \delta_{\bar{n}_+(\ell+1, j), q'}. \quad (40)$$

It can be useful for a clearer understanding, to encode this recursive formula for calculating dimensions into a cellular automaton. The automata works according the following steps:

(i) Draw a node for each link charge $q$ allowed.

(ii) Associate to every node the sector dimension $D_q(\ell)$. Starting step (zero sites): $D_q(0) = 1$ for every $q$.

(iii) For every local reduced basis state $|j_{\ell+1}\rangle_r$ evaluate $q = \bar{N}_{\ell+1} - \bar{n}_-(\ell+1, j)$ and $q' = \bar{n}_+(\ell+1, j)$. Then draw an arrow from node $q$ to node $q'$.

(iv) The new sector dimensions $D_{q'}(\ell + 1)$ are obtained by the sum, over all arrows that point to $q'$, of the dimension $D_q(\ell)$ of the node $q$ where that arrow starts from.

(v) Return to point 2 and iterate.

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 $\bar{D}(\ell) \propto \alpha^\ell$ but the basis $\alpha$ is strictly smaller than the reduced local space dimension, i.e. $\alpha < d$. Moreover, we found $\alpha$ to be even smaller than the total number of allowed local matter states $|s_{\psi}\rangle$. 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 site-by-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.

## 5.1. Example: U(1), spinless matter fermion, single rishon

This example corresponds to the QLM class introduce in paragraph 2.0.4, with the number of rishons per bond fixed to $N = 1$. Here, both local matter and local gauge fields are two-level systems. The gauge constraint allows only $d = 3$ states out of the $D = 8$ original ones to survive. Precisely, written as $|s_-, s_\psi, s_+\rangle$, they read

$$ \begin{aligned} |1\rangle_r &= |0, 1, 1\rangle_o & |2\rangle_r &= |1, 0, 1\rangle_o & |3\rangle_r &= |1, 1, 0\rangle_o & \text{on odd sites} \\ |1\rangle_r &= |1, 0, 0\rangle_e & |2\rangle_r &= |0, 1, 0\rangle_e & |3\rangle_r &= |0, 0, 1\rangle_e & \text{on even sites.} \end{aligned} \quad (41) $$

With the previous labeling of reduced basis states, the link constraint $Q_{r,x,x+1}$ becomes translationally-invariant (i.e. independent of $x$), and ultimately reads as
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Example of the cellular automata strategy to calculate the effective Hilbert dimension growth. This sketch corresponds to the two-level abelian-U(1) QLM, see section 5.0.1. The number written on node q at the stage ℓ of the automata corresponds to $D_q(\ell)$, i.e. the total number of many-body states on ℓ sites, satisfying all the gauge and link constraints, whose rightmost link charge is q.</figcaption></figure>

$$Q_{r,x,x+1} = |12\rangle\langle 12|_r + |13\rangle\langle 13|_r + |22\rangle\langle 22|_r + |23\rangle\langle 23|_r + |31\rangle\langle 31|_r, \quad (42)$$

with support dimension $\chi = 5$, or equivalently

$$V_{j,q} = (\delta_{j,1} + \delta_{j,2})\delta_{q,1} + \delta_{j,3}\delta_{q,2} \quad \text{and} \quad Z_{q,j} = \delta_{q,1}(\delta_{j,2} + \delta_{j,3}) + \delta_{q,2}\delta_{j,1}, \quad (43)$$

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\rangle_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\rangle_r$ is an arrow from $q = 0$ to $q = 0$, while state $|3\rangle_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 $\ell$ exactly as the Fibonacci sequence: in fact $D_1(\ell) = D_0(\ell - 1)$ while $D_0(\ell) = D_0(\ell - 1) + D_1(\ell - 1) = D_0(\ell - 1) + D_0(\ell - 2)$, and finally $\bar{D}(\ell) = D_0(\ell + 1)$. In conclusion we have $\bar{D}(\ell) = (\varphi^{\ell+3} - (1 - \varphi)^{\ell+3})/\sqrt{5}$ with $\varphi$ being the golden ratio $\varphi = (1 + \sqrt{5})/2$. This tells us that for large sizes $\ell$, the Hilbert dimension grows exponentially $\bar{D}(\ell) \propto \alpha^\ell$, but instead of using an exponential basis the local space dimension $d=3$ or the matter local dimension $d_\psi = 2$, it is $\alpha = (1 + \sqrt{5})/2 \approx 1.618$, i.e. the scaling is somewhat smoother.

## 5.2. 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 $\bar{N}$, while the matter is again a two-level quantum system. In this model there are $D = 2\bar{N}^2$ local states available, however only $d = 2\bar{N} + 1$ are allowed by the gauge constraint. We label them according to

$$|2k + 1\rangle_r = |k, 1, \bar{N} - k\rangle_o, \quad |2k\rangle_r = |k, 0, \bar{N} + 1 - k\rangle_o \quad (44)$$

on odd sites, while on even sites

$$|2k + 1\rangle_r = |\bar{N} - k, 0, k\rangle_e, \quad |2k\rangle_r = |\bar{N} - k, 1, k - 1\rangle_e, \quad (45)$$

where $|j\rangle_r$ spans within $1 \le j \le 2\bar{N} + 1$. Like in the previous example $Q_{r,x,x+1}$ is homogeneous and reads
<figure><img src="image_9.png" /><figcaption>**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]$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Fitted exponential basis $\alpha$ of the Hilbert space dimension growth rate, as a function of the selected number of rishons per link $\bar{N}$, in the U (1) QLM scenario. Inset: distance of $\alpha$ from 2, plotted as a function of $\bar{N}$, in double logarithmic scale.</figcaption></figure>

$$Q_{r,x,x+1} = |2\bar{N} + 1, 1\rangle\langle 2\bar{N} + 1, 1|_r + \sum_{k=1}^{\bar{N}} (|2k-1\rangle\langle 2k-1|_r + |2k\rangle\langle 2k|_r) \\ \otimes (|2\bar{N} - 2k + 3\rangle\langle 2\bar{N} - 2k + 3|_r + |2\bar{N} - 2k + 2\rangle\langle 2\bar{N} - 2k + 2|_r) \quad (46)$$

with support dimension $\chi = 4\bar{N} + 1$. The allowed charges $q$ here go from 0 to $\bar{N}$, so we have $\bar{N} + 1$ nodes in the cellular automata. The arrows are defined as follows: odd index states and even index states connect the nodes as

$$q_{\ell} = k \xrightarrow{|2k+1\rangle} q_{\ell+1} = \bar{N} - k \quad \text{and} \quad q_{\ell} = k \xrightarrow{|2k\rangle} q_{\ell+1} = \bar{N} + 1 - k. \quad (47)$$

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 $\bar{D}(\ell)$ with the system size $\ell$, for several different rishon number choices $\bar{N}$. In every case considered starting from sizes of $\ell \sim 10$, $\bar{D}(\ell)$ matches an exponential scaling in $\ell$. We fitted the exponential basis $\alpha(\bar{N})$ $\bar{D}(\ell) \propto \alpha^\ell(\bar{N})$ in the interval $\ell \in [100, 1000]$. A smooth, monotonic behavior of $\alpha$ as a function of $\bar{N}$ is observed and reported in figure 8. The calculated $\alpha$ values are never greater than 2, and saturate to 2 in a polynomial fashion with increasing numbers of rishons per link $\bar{N}$, i.e. when approaching the Wilson limit (as shown in figure 8, inset). This study reveals that in
the proximity of the thermodynamical limit, the U (1) quantum link model will never require more computational resources or store more quantum information than an unconstrained spin-1/2 model, thus representing a comparative bound on the algebraic complexity of the quantum link model. This result seems to imply that, in some cases, a U (1) quantum link model (with spinless fermionic matter) could in principle be mapped into a spin-1/2 model with constraints, and these constraints vanish in the Wilson limit, where one should recover the corresponding unconstrained spin-1/2 model. Although our argument based on Hilbert dimension scaling does not provide any information whether such a mapping is actually possible for a given U (1) gauge theory, some examples where the mapping exists are known. For instance, it was shown that the Schwinger model can be mapped to a long range interacting spin-1/2 model [20].

### 5.3. Example: U(2), spin-1/2 fermions, single rishon

In this example, which corresponds to the scenario we already introduced in section 2.0.4, both matter and link fields host spin-1/2 fermionic excitations. The U (2) gauge constraint fixes the available vertex states $|j\rangle_r$ to be both in a SU (2) spin singlet and in a defined U (1) occupation number: here we consider the case $n_{x,\psi} + n_{x,-} + n_{x,+} = 2$. Moreover, the link constraint fixes the number $\tilde{N}$ of rishons on a lattice bond in this example to $\tilde{N} = 1$. With the aforementioned choices, the reduced local space is four-dimensional:

$$ 
\begin{aligned}
|1\rangle_r &= \frac{1}{\sqrt{2}} (|0, \uparrow, \downarrow\rangle - |0, \downarrow, \uparrow\rangle), & |2\rangle_r &= \frac{1}{\sqrt{2}} (|\uparrow, 0, \downarrow\rangle - |\downarrow, 0, \uparrow\rangle), \\
|3\rangle_r &= \frac{1}{\sqrt{2}} (|\uparrow, \downarrow, 0\rangle - |\downarrow, \uparrow, 0\rangle), & |4\rangle_r &= |0, \phi, 0\rangle,
\end{aligned} 
\quad (48) $$

where $|\uparrow\rangle \equiv c_{\uparrow}^{\dagger} |0\rangle$, $|\downarrow\rangle \equiv c_{\downarrow}^{\dagger} |0\rangle$ and $|\phi\rangle \equiv c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |0\rangle$. The reduced link projector $Q_{r,x,x+1}$ has support dimension $\chi = 8$ and reads

$$ Q_{r,x,x+1} = (|1\rangle\langle 1| + |2\rangle\langle 2|) \otimes (|1\rangle\langle 1| + |4\rangle\langle 4|) + \\ (|3\rangle\langle 3| + |4\rangle\langle 4|) \otimes (|2\rangle\langle 2| + |3\rangle\langle 3|). \quad (49) $$

The corresponding cellular automata appears as shown in figure 9, left panel. The Hilbert space dimension scaling in this example is exactly an exponential, specifically $\tilde{D}(\ell) = 2^{\ell+1}$ and ultimately $\alpha = 2$. This again implies that a mapping of this class of models to a spin-1/2 system might exist, even though it has not presently been proved to the best of our knowledge.

### 5.4. Example: U(2), spin-1/2 fermions, double rishon

The last scenario we discuss explicitly is again the U (2) case with spin-1/2 particles, but now we set $\tilde{N} = 2$ rishons on every link. The effective link constraint in the reduced formulation becomes translationally invariant when we start from a staggered original fermion filling in the vertices: the choice that describes more rich physics is given by $n_{x,\psi} + n_{x,-} + n_{x,+} = 3 - (-1)^x$, which results in $d = 6$ reduced states. Indeed, on odd sites we have
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Cellular automata scheme for the U (2) QLM with spin-$\frac{1}{2}$ fermions and single rishon per link (a) and for the double rishon per link (b) scenario. (c) Scaling of the total Hilbert space dimension $\bar{D}(\ell)$ as a function of the 1D chain length $\ell$, in the U (2) double rishon case, evaluated numerically (black dots). The cyan line is an exponential fit, revealing the exponential basis $\alpha \approx 2.246 9796$.</figcaption></figure>

$$ 
\begin{align} 
|1\rangle_r &= |\phi, \phi, 0\rangle_o, & |2\rangle_r &= \frac{1}{\sqrt{2}}(|\phi, \uparrow, \downarrow\rangle_o - |\phi, \downarrow, \uparrow\rangle_o), \nonumber \\ 
|3\rangle_r &= \frac{1}{\sqrt{2}}(|\uparrow, \phi, \downarrow\rangle_o - |\downarrow, \phi, \uparrow\rangle_o), & |4\rangle_r &= \frac{1}{\sqrt{2}}(|\uparrow, \downarrow, \phi\rangle_o - |\downarrow, \uparrow, \phi\rangle_o), \nonumber \\ 
|5\rangle_r &= |\phi, 0, \phi\rangle_o, & |6\rangle_r &= |0, \phi, \phi\rangle_o, \tag{50}
\end{align} 
$$

while on even sites

$$ 
\begin{align} 
|1\rangle_r &= |0, 0, \phi\rangle_e, & |2\rangle_r &= \frac{1}{\sqrt{2}}(|0, \uparrow, \downarrow\rangle_e - |0, \downarrow, \uparrow\rangle_e), \nonumber \\ 
|3\rangle_r &= \frac{1}{\sqrt{2}}(|\uparrow, 0, \downarrow\rangle_e - |\downarrow, 0, \uparrow\rangle_e), & |4\rangle_r &= \frac{1}{\sqrt{2}}(|\uparrow, \downarrow, 0\rangle_e - |\downarrow, \uparrow, 0\rangle_e), \nonumber \\ 
|5\rangle_r &= |0, \phi, 0\rangle_e, & |6\rangle_r &= |\phi, 0, 0\rangle_e. \tag{51}
\end{align} 
$$

As there are three possible occupations on the link degree of freedom, the automata has three charges $q$. The intermediate charge $q = 0$ connects $|1\rangle_r$ (to the left) to $|6\rangle_r$ (to the right), charge $q = 1$ connects $\{|2\rangle_r, |3\rangle_r\}$ to $\{|3\rangle_r, |4\rangle_r\}$, and finally $q = 2$ connects $\{|4\rangle_r, |5\rangle_r, |6\rangle_r\}$ to $\{|1\rangle_r, |2\rangle_r, |5\rangle_r\}$, for a total reduced link projector support of dimension $\chi = 14 < d^2 = 36$.

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 $\bar{D}(\ell)$ for system sizes up to $\ell \sim 850$. Once again, an asymptotically exponential scaling $\bar{D}(\ell) \propto \alpha^\ell$ 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 $\alpha \approx 2.2470$.
## 6. 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 out-of-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. Finite temperature, open system dynamics [64–66], richer tensor structures [67] and optimally controlled dynamics might be studied in the future [68, 69].

## Acknowledgments

Authors acknowledge support from EU through SIQS, ERC-St Grant ColdSIM (no. 307688), EOARD, UdS via Labex NIE and IdEX, and the German Research Foundation (DFG) via the SFB/TRR21. Authors thank M Dalmonte and G Pupillo for stimulating discussions.

## References

[1] DeGrand T and DeTar C 2006 *Lattice Methods for Quantum Chromodynamics* (Singapore: World Scientific)

[2] Rothe H J 2012 *Lattice Gauge Theories: An Introduction* (Singapore: World Scientific)

[3] Wilson K 1974 *Phys. Rev. D* **10** 2445

[4] Kogut J and Susskind L 1975 *Phys. Rev. D* **11** 395-408

[5] Lee P A, Nagaosa N and Wen X-G 2006 *Rev. Mod. Phys.* **78** 17

[6] Lacroix C, Mendels P and Mila F (ed) 2011 *Introduction to Frustrated Magnetism* (Springer Series in Solid-State Sciences vol 164) (Berlin: Springer)

[7] Fradkin E 2013 *Field Theories of Condensed Matter Physics* (Cambridge: Cambridge University Press)

[8] Kitaev A 2003 *Ann. Phys.* **303** 2

[9] Nayak C, Simon S H, Stern A, Freedman M and das Sarma S 2008 *Rev. Mod. Phys.* **80** 1083

[10] Trabesinger A 2012 *Nat. Phys. Insight* **8** 263

[11] Zohar E, Cirac J I and Reznik B 2012 *Phys. Rev. Lett.* **109** 125302

[12] Banerjee D, Dalmonte M, Muller M, Rico E, Stebler P, Wiese U-J and Zoller P 2012 *Phys. Rev. Lett.* **109** 175302

[13] Tagliacozzo L, Celi A, Zamora A and Lewenstein M 2013 *Ann. Phys.* **330** 160
[14] Marcos D, Rabl P, Rico E and Zoller P 2013 *Phys. Rev. Lett.* **111** 110504

[15] Hauke P, Marcos D, Dalmonte M and Zoller P 2013 *Phys. Rev. X* **3** 041018

[16] Wiese U-J 2013 *Ann. Phys* **525** 777

[17] Zohar E, Cirac J I and Reznik B 2013 *Phys. Rev. A* **88** 023617

[18] Stannigel K, Hauke P, Marcos D, Hafezi M, Diehl S, Dalmonte M and Zoller P 2014 *Phys. Rev. Lett.* **112** 120406

[19] Shimizu Y and Kuramashi Y 2014 *Phys. Rev. D* **90** 014508

[20] Rico E, Pichler T, Dalmonte M, Zoller P and Montangero S 2014 *Phys. Rev. Lett.* **112** 201601

[21] White S R 1992 *Phys. Rev. Lett.* **69** 2863

[22] Fannes M, Nachtergaele B and Werner R F 1992 *Comm. Math. Phys.* **144** 443

[23] Östlund S and Rommer S 1995 *Phys. Rev. Lett.* **75** 3537

[24] Dukelsky J, Martin-Delgado M A, Nishino T and Sierra G 1998 *Europhys. Lett.* **43** 457

[25] Vidal G 2003 *Phys. Rev. Lett.* **91** 147902

[26] Verstraete F, Porras D and Cirac J I 2004 *Phys. Rev. Lett.* **93** 227205

[27] Schollwöck U 2005 *Rev. Mod. Phys.* **77** 259

[28] Perez-Garcia D, Verstraete F, Wolf M M and Cirac J I 2007 *Quantum Inf. Comput.* **7** 401

[29] de Chiara G, Rizzi M, Rossini D and Montangero S 2008 *J. Comput. Theor. Nanosci.* **5** 1277

[30] Verstraete F, Wolf M M, Perez-Garcia D and Cirac J I 2006 *Phys. Rev. Lett.* **96** 220601

[31] Murg V, Legeza Ö, Noack R M and Verstraete F 2010 *Phys. Rev. B* **82** 205105

[32] Vidal G 2007 *Phys. Rev. Lett.* **99** 220405

[33] Verstraete F and Cirac J I 2010 *Phys. Rev. Lett.* **104** 190405

[34] Singh S, Pfeifer R N C and Vidal G 2010 *Phys. Rev. A* **82** 050301

[35] Singh S, Pfeifer R N C and Vidal G 2011 *Phys. Rev. B* **83** 115125

[36] Singh S and Vidal G 2012 *Phys. Rev. B* **86** 195114

[37] Singh S and Vidal G 2013 *Phys. Rev. B* **88** 121108

[38] Bauer B, Corboz P, Orus R and Troyer M 2011 *Phys. Rev. B* **83** 125106

[39] Weichselbaum A 2012 *Ann. Phys.* **327** 2972

[40] Corboz P, Lajkó M, Penc K, Mila F and Läuchli A M 2013 *Phys. Rev. B* **87** 195113

[41] Wegner F J 1971 *J. Math. Phys.* **12** 2259

[42] Ardonne E, Fendley P and Fradkin E 2004 *Ann. Phys.* **310** 493

[43] Martin-Delgado M A and Sierra G 1999 *Phys. Rev. Lett.* **83** 1514

[44] Byrnes T M R, Sriganesh P, Bursill R J and Hamer C J 2002 *Phys. Rev. D* **66** 013002

[45] Buyens B, Haegeman J, van Acoleyen K, Verschelde H and Verstraete F 2014 *Phys. Rev. Lett.* **113** 091601

[46] Haegeman J, Cirac J I, Osborne T J, Verschelde H and Verstraete F 2010 *Phys. Rev. Lett.* **105** 251601

[47] Tagliacozzo L and Vidal G 2011 *Phys. Rev. B* **83** 115127

[48] Banuls M C, Cichy K, Jansen K and Cirac J I 2013 *J. High Energy Phys.* **11** 158

[49] Silvi P, de Chiara G, Calarco T, Morigi G and Montangero S 2013 *Ann. Phys.* **525** 827

[50] Horn D 1981 *Phys. Lett. B* **100** 149

[51] Orland P and Rohrlich D 1990 *Nucl. Phys. B* **338** 647

[52] Chandrasekharan S and Wiese U-J 1997 *Nucl. Phys. B* **492** 455

[53] Brower R, Chandrasekharan S and Wiese U-J 1999 *Phys. Rev. D* **60** 094502

[54] Brower R, Chandrasekharan S, Riederer S and Wiese U-J 2004 *Nucl. Phys. B* **693** 149-75

[55] Misguich G and Lhuillier C 2003 *Frustrated Spin Systems* 2nd edn (Singapore: World Scientific) pp 235-319

[56] Bramwell S T and Gingras M J P 2001 *Science* **294** 1495-501

[57] Castelnovo C, Moessner R and Sondhi S L 2008 *Nature* **451** 42-45

[58] Balents L 2010 *Nature* **464** 199-208

[59] Banerjee D, Bögli M, Hofmann C P, Jiang F-J, Widmer P and Wiese U-J arXiv:1406.2077

[60] Glaetzle A W, Dalmonte M, Nath R, Rousochatzakis I, Moessner R and Zoller P arXiv:1404.5326

[61] Kogut J 1979 *Rev. Mod. Phys.* **51** 659
[62] Ignacio Cirac J 2012 arXiv:1205.3742

[63] Elitzur S 1975 *Phys. Rev. D* **12** 3978

[64] Verstraete F, Garcia-Rippol J J and Cirac J I 2004 *Phys. Rev. Lett.* **93** 207204

[65] Zwolak M and Vidal G 2004 *Phys. Rev. Lett.* **93** 207205

[66] Jaschke D, Silvi P, Calarco T and Montangero S in preparation

[67] Rizzi M, Montangero S and Vidal G 2008 *Phys. Rev. A* **77** 052328

[68] Caneva T, Calarco T and Montangero S 2011 *Phys. Rev. A* **84** 022326

[69] Doria P, Calarco T and Montangero S 2011 *Phys. Rev. Lett.* **106** 190501