Lattice gauge tensor networks

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, $\{{{\psi }_{x}},\psi _{y}^{\dagger }\}={{\delta }_{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)$ 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 }_{x}}}^{ab}$ exist on the links of the lattice $\langle x,x+{{\mu }_{x}}\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 }_{x}}}^{ab}=c_{x+{{\mu }_{x}},-{{\mu }_{x}}}^{a\dagger }c_{x,+{{\mu }_{x}}}^{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}^{\prime }}^{a\dagger }]}_{\pm }}={{\delta }_{a,b}}{{\delta }_{x,y}}{{\delta }_{{{\mu }_{x}},\mu _{y}^{\prime }}}.$ 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 ${{\mathcal{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 $[{{\mathcal{N}}_{x,x+{{\mu }_{x}}}},U_{y,y+{{\mu }_{y}}}^{ab}]=0$, from which it follows that $[{{\mathcal{N}}_{x,x+\mu }},G_{y}^{\nu }]=[{{\mathcal{N}}_{x,x+\mu }},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 ${{\mathcal{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 $\bar{N}$ one selects. In any case, we restrict the Hilbert space to the 'physical' states $|{{\varphi }_{{\rm phys}}}\rangle$ which satisfy ${{\mathcal{N}}_{x,x+{{\mu }_{x}}}}|{{\varphi }_{{\rm phys}}}\rangle =|{{\varphi }_{{\rm phys}}}\rangle {{\bar{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 ${{\bar{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.

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,
  $$\begin{eqnarray}&&\tilde{U}_{x,x+{{\mu }_{x}}}^{ab}={{{\rm{e}} }^{{\rm{i}} \theta {{E}_{x,x+{{\mu }_{x}}}}}}U_{x,x+{{\mu }_{x}}}^{ab}{{{\rm{e}} }^{-{\rm{i}} \theta {{E}_{x,x+{{\mu }_{x}}}}}}={{{\rm{e}} }^{{\rm{i}} \theta }}U_{x,x+{{\mu }_{x}}}^{ab},\end{eqnarray} \tag{ 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 Ω_ak acting on the left or on the right of the field as
  $$\begin{eqnarray}\begin{array}{rcl} \tilde{U}_{x,x+{{\mu }_{x}}}^{ab} & = & {{{\rm{e}} }^{{\rm{i}} \mathop{\sum }\limits_{\nu }{{\theta }^{\nu }}L_{x,x+{{\mu }_{x}}}^{\nu }}}U_{x,x+{{\mu }_{x}}}^{ab}{{{\rm{e}} }^{-{\rm{i}} \mathop{\sum }\limits_{\nu }{{\theta }^{\nu }}L_{x,x+{{\mu }_{x}}}^{\nu }}}=\mathop{\sum }\limits_{k}{{\Omega }_{ak}}U_{x,x+{{\mu }_{x}}}^{kb} \\ \bar{U}_{x,x+{{\mu }_{x}}}^{ab} & = & {{{\rm{e}} }^{{\rm{i}} \mathop{\sum }\limits_{\nu }{{\theta }^{\nu }}R_{x,x+{{\mu }_{x}}}^{\nu }}}U_{x,x+{{\mu }_{x}}}^{ab}{{{\rm{e}} }^{-{\rm{i}} \mathop{\sum }\limits_{\nu }{{\theta }^{\nu }}R_{x,x+{{\mu }_{x}}}^{\nu }}}=\mathop{\sum }\limits_{k}U_{x,x+{{\mu }_{x}}}^{ak}\Omega _{bk}^{*}, \\ \end{array}\end{eqnarray} \tag{ 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}]{\mkern 1mu} ={{\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 }}]=2i{{f}_{\mu \nu \omega }}{{\lambda }^{\omega }}$, with ${{f}_{\mu \nu \omega }}$ the structure constants of the SU(N) algebra and ${\rm 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{eqnarray}\begin{array}{rcl} {{G}_{x}} & = & \psi _{x}^{\dagger }{{\psi }_{x}}+\mathop{\sum }\limits_{{{\mu }_{x}}}\left( {{E}_{x-{{\mu }_{x}},x}}-{{E}_{x,x+{{\mu }_{x}}}} \right) \\ {} & = & \psi _{x}^{\dagger }{{\psi }_{x}}+\mathop{\sum }\limits_{{{\mu }_{x}}}\left( {{n}_{x,-{{\mu }_{x}}}}+{{n}_{x,+{{\mu }_{x}}}} \right)-\frac{{{{\bar{N}}}_{x-{{\mu }_{x}},x}}+{{{\bar{N}}}_{x,x+{{\mu }_{x}}}}}{2}, \\ \end{array}\end{eqnarray} \tag{ 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 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
  $$\begin{eqnarray}&&[\psi _{x}^{\dagger }{{\psi }_{x}}+\mathop{\sum }\limits_{{{\mu }_{x}}}({{n}_{x,-{{\mu }_{x}}}}+{{n}_{x,+{{\mu }_{x}}}})]|{{\varphi }_{{\rm phys}}}\rangle =|{{\varphi }_{{\rm phys}}}\rangle \times {\rm constant},\end{eqnarray} \tag{ 4 }$$
  which is equivalent to Gauss' law, and the physical states $|{{\varphi }_{{\rm 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 }]={\rm{i}} {{\epsilon }^{\mu \nu \omega }}G_{x}^{\omega }{{\delta }_{x,y}}$ with
  $$\begin{eqnarray}\begin{array}{rcl} G_{x}^{\nu } & = & \mathop{\sum }\limits_{ab}\psi _{x}^{a\dagger }\frac{\lambda _{ab}^{\nu }}{2}\psi _{x}^{b}+\mathop{\sum }\limits_{{{\mu }_{x}}}\left( R_{x-{{\mu }_{x}},x}^{\nu }+L_{x,x+{{\mu }_{x}}}^{\nu } \right) \\ {} & = & \mathop{\sum }\limits_{ab}\left[ \psi _{x}^{a\dagger }\frac{\lambda _{ab}^{\nu }}{2}\psi _{x}^{b}+\mathop{\sum }\limits_{{{\mu }_{x}}}\left( c_{x,-{{\mu }_{x}}}^{a\dagger }\frac{\lambda _{ab}^{\nu }}{2}c_{x,-{{\mu }_{x}}}^{b}+c_{x,+{{\mu }_{x}}}^{a\dagger }\frac{\lambda _{ab}^{\nu }}{2}c_{x,+{{\mu }_{x}}}^{b} \right) \right] \\ \end{array}\end{eqnarray} \tag{ 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 }_{{\rm 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 ψ_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

$$\begin{eqnarray}&&|{{j}_{x}}{{\rangle }_{r}}=\mathop{\sum }\limits_{{{s}_{\psi }},{{{\vec{s}}}_{{{\mu }_{x}}}}}A_{{{s}_{\psi }},{{{\vec{s}}}_{{{\mu }_{x}}}}}^{[x]j}|{{s}_{\psi }},{{\vec{s}}_{{{\mu }_{x}}}}{{\rangle }_{x}},\end{eqnarray} \tag{ 6 }$$

once a canonical basis for the matter field $|{{s}_{\psi }}{{\rangle }_{x}}$ and one for each rishon field $|{{s}_{{{\mu }_{x}}}}{{\rangle }_{x,{{\mu }_{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.

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. $[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{eqnarray}\begin{array}{rcl} {{H}_{{\rm pure}}} & = & {{H}_{{\rm electric}}}+{{H}_{{\rm magn}}} \\ {} & = & \mathop{\sum }\limits_{x,{{\mu }_{x}}}\{g_{{\rm abel}}^{2}{{({{E}_{x,x+{{\mu }_{x}}}})}^{2}}+g_{{\rm non}-{\rm ab}}^{2}\mathop{\sum }\limits_{\nu }[{{(L_{x,x+{{\mu }_{x}}}^{\nu })}^{2}}+{{(R_{x,x+{{\mu }_{x}}}^{\,\nu })}^{2}}]\} \\ {} & {} & -\frac{1}{g_{{\rm magn}}^{2}}\mathop{\sum }\limits_{x,{{\mu }_{x}},{{\mu }_{y}}}[{\rm Tr}({{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}})+{\rm{h}} .{\rm{c}} .]. \\ \end{array}\end{eqnarray} \tag{ 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_{{\rm abel}}^{2}$, $g_{{\rm non}-{\rm ab}}^{2}$ and $g_{{\rm 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_{{\rm magn}}^{2}=4a{{g}^{2}}$, and $g_{{\rm abel}}^{2}=g_{{\rm non}-{\rm 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}_{{\rm abel}}}\simeq {{g}_{{\rm non}-{\rm ab}}}\gg 1/{{g}_{{\rm 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,

$$\begin{eqnarray}&&{{H}_{{\rm coup}}}=\mathop{\sum }\limits_{x,{{\mu }_{x}}}{{J}_{x,{{\mu }_{x}}}}(\psi _{x}^{\dagger }{{U}_{x,x+{{\mu }_{x}}}}{{\psi }_{x+{{\mu }_{x}}}}+{\rm{h}} .{\rm{c}} .)+\mathop{\sum }\limits_{x}{{m}_{x}}\psi _{x}^{\dagger }{{\psi }_{x}}\end{eqnarray} \tag{ 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.

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

$$\begin{eqnarray}&&H=J\mathop{\sum }\limits_{x}\left( \psi _{x}^{\dagger }{{U}_{x,x+1}}{{\psi }_{x+1}}+{\rm{h}} .{\rm{c}} . \right)+{{g}^{2}}\mathop{\sum }\limits_{x}{{\left( {{E}_{x,x+1}} \right)}^{2}}+m\mathop{\sum }\limits_{x}{{\left( -1 \right)}^{x}}\psi _{x}^{\dagger }{{\psi }_{x}}\end{eqnarray} \tag{ 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² 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 }_{{\rm phys}}}\rangle =({{n}_{x+1,-}}+{{n}_{x,+}})|{{\varphi }_{{\rm phys}}}\rangle =\bar{N}\;|{{\varphi }_{{\rm phys}}}\rangle$
• (ii)  
  Gauss' law on every vertex: $\left( \psi _{x}^{\dagger }{{\psi }_{x}}+{{n}_{x,-}}+{{n}_{x,+}} \right)|{{\varphi }_{{\rm phys}}}\rangle =|{{\varphi }_{{\rm phys}}}\rangle \left( \bar{N}-\frac{1+{{\left( -1 \right)}^{x}}}{2} \right)$

The factor $(\bar{N}-\frac{1+{{\left( -1 \right)}^{x}}}{2})$ 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 $\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 ${{\mathcal{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 ϑ 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$,
  $$\begin{eqnarray}\begin{array}{rcl} \frac{1}{\sqrt{\frac{{\bar{N}}}{2}(\frac{{\bar{N}}}{2}+1)}}U|\bar{N},E\rangle & \to & |\bar{N},E+1\rangle ;\;\frac{1}{\frac{{\bar{N}}}{2}(\frac{{\bar{N}}}{2}+1)}[U,{{U}^{\dagger }}]\to 0; \\ \qquad \qquad \frac{1}{\sqrt{\frac{{\bar{N}}}{2}(\frac{{\bar{N}}}{2}+1)}}[E,U] & = & \frac{1}{\sqrt{\frac{{\bar{N}}}{2}(\frac{{\bar{N}}}{2}+1)}}U \\ \end{array}\end{eqnarray} \tag{ 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

$$\begin{eqnarray}\begin{array}{rcl} G_{x}^{\nu } & = & R_{x}^{\nu }+\mathop{\sum }\limits_{a,b}\left( \psi _{x}^{a\dagger }\frac{\sigma _{ab}^{\nu }}{2}\psi _{x}^{b} \right)+L_{x}^{\nu } \\ {} & {} & =\mathop{\sum }\limits_{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) \\ \end{array}\end{eqnarray} \tag{ 11 }$$

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

$$\begin{eqnarray}\begin{array}{rcl} H & = & \frac{1}{2}\mathop{\sum }\limits_{x}\left( g_{{\rm{a}} }^{2}\;E_{x}^{2}+g_{{\rm na}}^{2}\mathop{\sum }\limits_{\nu }\left[ {{(R_{x}^{\nu })}^{2}}+{{(L_{x}^{\nu })}^{2}} \right] \right)+m\mathop{\sum }\limits_{x,a}{{\left( -1 \right)}^{x}}\psi _{x}^{a\dagger }\psi _{x}^{a} \\ {} & {} & +t\mathop{\sum }\limits_{x,a,b}\left[ \psi _{x}^{a\dagger }U_{x,x+1}^{ab}\psi _{x+1}^{b}+{\rm{h}} .{\rm{c}} . \right] \\ \end{array}\end{eqnarray} \tag{ 12 }$$

The $g_{{\rm{a}} }^{2}$ and $g_{{\rm 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

$$\begin{eqnarray}&&{{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},\end{eqnarray} \tag{ 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},{{\mathcal{N}}_{x,x+\mu }}]=0$. The link sector selected ($\bar{N}=1$, i.e. $\mathcal{N}|{{\varphi }_{{\rm phys}}}\rangle =|{{\varphi }_{{\rm phys}}}\rangle$), recovers exactly a two-level system on every link, as shown in figure 2.

Zoom In Zoom Out Reset image size

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}{\mkern 1mu} +\sigma _{x,x+{{\mu }_{y}}}^{z}+\sigma _{x-{{\mu }_{x}},x}^{z}+\sigma _{x-{{\mu }_{y}},x}^{z})|{{\varphi }_{{\rm phys}}}\rangle =-|{{\varphi }_{{\rm 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:

$$\begin{eqnarray}&&{{H}_{{\rm 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}}.\end{eqnarray} \tag{ 14 }$$

Effectively this interaction projects the Hilbert space to the zero magnetization subspace ${{G}_{x}}|{{\varphi }_{{\rm 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 }_{{\rm 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}}$.

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.

$$\begin{eqnarray}&&[{{n}_{x,-}},G_{{{x}^{\prime }}}^{\nu }]=[{{n}_{x,+}},G_{{{x}^{\prime }}}^{\nu }]=0.\end{eqnarray} \tag{ 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}}{{|}_{{\rm phys}}}\to |{{j}_{x}}{{\rangle }_{r}}$, precisely:

$$\begin{eqnarray}&&{{n}_{x,-}}|{{j}_{x}}{{\rangle }_{r}}=|{{j}_{x}}{{\rangle }_{r}}\cdot {{\bar{n}}_{-}}(x,j)\quad {\rm and}\quad {{n}_{x,+}}|{{j}_{x}}{{\rangle }_{r}}=|{{j}_{x}}{{\rangle }_{r}}\cdot {{\bar{n}}_{+}}(x,j).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{{Q}_{r,x,x+1}}|{{j}_{x}}\;j_{x+1}^{\prime }{{\rangle }_{r}}=|{{j}_{x}}\;j_{x+1}^{\prime }{{\rangle }_{r}}\ {{\delta }_{{{{\bar{n}}}_{+}}(x,j)+{{{\bar{n}}}_{-}}(x+1,{{j}^{\prime }}),{{{\bar{N}}}_{x}}}}=|{{j}_{x}}\;j_{x+1}^{\prime }{{\rangle }_{r}}\mathop{\sum }\limits_{q=0}^{{{{\bar{N}}}_{x,x+1}}}V_{{{j}_{x}},q}^{[x]}\cdot Z_{j_{x+1}^{\prime },q}^{[x+1]}\end{eqnarray} \tag{ 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:

$$\begin{eqnarray}&&F_{{{j}_{x}},j_{x}^{\prime }}^{[x]{{q}_{x-1}},{{q}_{x}}}={{\delta }_{{{j}_{x}},j_{x}^{\prime }}}\cdot Z_{{{q}_{x-1}},{{j}_{x}}}^{[x]}\cdot V_{{{j}_{x}},{{q}_{x}}}^{[x]},\end{eqnarray} \tag{ 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}}$.

Zoom In Zoom Out Reset image size

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.

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

$$\begin{eqnarray}\begin{array}{rcl} \rho ({{t}_{0}}+t) & = & {{{\rm{e}} }^{{\rm{i}} \bar{H}t}}\rho (t){{{\rm{e}} }^{-{\rm{i}} \bar{H}t}}\quad {\rm{for}} \;{\rm{real}} -{\rm{time,}} \;{\rm{or}} \\ \rho ({{\beta }_{0}}+\beta ) & = & {{{\rm{e}} }^{-\bar{H}\beta /2}}\rho ({{\beta }_{0}}){{{\rm{e}} }^{-\bar{H}\beta /2}}\qquad {\rm{for}} \;{\rm{imaginary}} -{\rm{time}} . \\ \end{array}\end{eqnarray} \tag{ 29 }$$

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 $\rho =X{{X}^{\dagger }}$, which always exists since $\rho \gt 0$ and we encode X as an MPO. The time evolution can be then carried out directly on X, because applying

$$\begin{eqnarray}\begin{array}{rcl} X({{t}_{0}}+t) & = & {{{\rm{e}} }^{{\rm{i}} \bar{H}t}}X({{t}_{0}})\quad {\rm{for}} \;{\rm{real}} -{\rm{time,}} \;{\rm{or}} \\ X({{\beta }_{0}}+\beta ) & = & {{{\rm{e}} }^{-\bar{H}\beta /2}}X({{\beta }_{0}})\qquad {\rm{for}} \;{\rm{imaginary}} -{\rm{time}} , \\ \end{array}\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}\begin{array}{rcl} {\rm exp} \left( \gamma \mathop{\sum }\limits_{x}H_{x,x+1}^{[r]} \right) & = & \left( \mathop{\otimes }\limits_{x}{{{\rm{e}} }^{{{c}_{1}}\gamma H_{2x-1,2x}^{[r]}}} \right)\left( \mathop{\otimes }\limits_{x}{{{\rm{e}} }^{{{d}_{1}}\gamma H_{2x,2x+1}^{[r]}}} \right) \\ {} & {} & \times \left( \mathop{\otimes }\limits_{x}{{{\rm{e}} }^{{{c}_{2}}\gamma H_{2x-1,2x}^{[r]}}} \right)\ldots \left( \mathop{\otimes }\limits_{x}{{{\rm{e}} }^{{{d}_{p}}\gamma H_{2x,2x+1}^{[r]}}} \right)+O({{\gamma }^{p}}) \\ \end{array}\end{eqnarray} \tag{ 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 $\bar{Q}{\rm exp} (\gamma {{\sum }_{x}}{{\bar{H}}_{x,x+1}})$. More generally one might want to apply the link projector $\bar{Q}$ either before, after or before and after the evolution, since ${{\bar{Q}}^{2}}=\bar{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 $\left[ {{Q}_{r,x,x+1}},{{H}_{r,{{x}^{\prime }},{{x}^{\prime }}+1}} \right]=0$, i.e. even local Hamiltonian terms commute with the link constraints. Indeed,

$$\begin{eqnarray}\begin{array}{rcl} \left[ {{Q}_{r,x,x+1}},{{H}_{r,{{x}^{\prime }},{{x}^{\prime }}+1}} \right] & = & \left[ \bar{A}{{Q}_{x,x+1}}{{{\bar{A}}}^{\dagger }},\bar{A}{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}{{{\bar{A}}}^{\dagger }} \right] \\ {} & = & \bar{A}\left( {{Q}_{x,x+1}}\bar{P}{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}-{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}\bar{P}{{Q}_{x,x+1}} \right){{{\bar{A}}}^{\dagger }} \\ {} & = & \bar{A}\left( \bar{P}{{Q}_{x,x+1}}{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}-{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}{{Q}_{x,x+1}}\bar{P} \right){{{\bar{A}}}^{\dagger }} \\ {} & = & \bar{A}\left[ {{Q}_{x,x+1}},{{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}} \right]{{{\bar{A}}}^{\dagger }}, \\ \end{array}\end{eqnarray} \tag{ 32 }$$

as $\bar{A}\bar{P}=\bar{A}$. In the original basis ${{Q}_{x,x+1}}$ and ${{H}_{{{x}^{\prime }},{{x}^{\prime }}+1}}$ act on common degrees of freedom only if $x={{x}^{\prime }}$, 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{eqnarray}\begin{array}{rcl} \bar{Q}{\rm exp} \left( \gamma \mathop{\sum }\limits_{x}{{{\bar{H}}}_{x,x+1}} \right) & = & \left( \mathop{\otimes }\limits_{x}M_{2x-1,2x}^{[r,1]} \right)\left( \mathop{\otimes }\limits_{x}M_{2x,2x+1}^{[r,1]} \right) \\ {} & {} & \times \left( \mathop{\otimes }\limits_{x}M_{2x-1,2x}^{[r,2]} \right)\ldots \left( \mathop{\otimes }\limits_{x}M_{2x,2x+1}^{[r,p]} \right)+O({{\gamma }^{p}}) \\ \end{array}\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}&&M_{2x-1,2x}^{[r,\nu ]}={{Q}_{r,2x-1,2x}}\ {\rm exp} \left( {{c}_{\nu }}\gamma H_{2x-1,2x}^{[r]} \right),\quad M_{2x,2x+1}^{[r,\nu ]}={{Q}_{r,2x,2x+1}}\ {\rm exp} \left( {{d}_{\nu }}\gamma H_{2x,2x+1}^{[r]} \right).\end{eqnarray} \tag{ 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.

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 =X{{X}^{\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

$$\begin{eqnarray}&&X=\sum 1\ldots {{j}_{L}}=1dj\sum 1\ldots {{k}_{L}}=1bk\sum 1\ldots {{w}_{L-1}}=1mwX_{{{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}}}\ldots X_{{{j}_{L}},{{k}_{L}}}^{[L]{{w}_{L-1}}}\ |{{j}_{1}}\ldots {{j}_{L}}{{\rangle }_{r}}\langle {{k}_{1}}\ldots {{k}_{L}}{{|}_{b}},\end{eqnarray} \tag{ 35 }$$

Zoom In Zoom Out Reset image size

where the correlation bondlink dimension m and the bath bondlink dimension b are both arbitrary (although $b\leqslant {{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}^{\prime }},q)}_{x,x+1}}$ so that the corresponding state $|{{j}_{x}},j_{x+1}^{\prime }\rangle$, belongs to the support of ${{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}^{\prime }})$. All these triplets are collected into the set ${{\Omega }_{x,x+1}}$, and their number is $\chi =\#{{\Omega }_{x,x+1}}$, clearly with $\chi \lt {{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

  $$\begin{eqnarray}&&\Gamma _{{{w}_{x+1}},{{k}_{x+1}},{{j}_{x+1}}}^{{{w}_{x-1}},{{k}_{x}},{{j}_{x}}}=\sum xmw\sum ^{\prime} ,j^{\prime\prime} dj\left( M_{x,x+1}^{[r,\nu ]} \right)_{{{j}_{x}},{{j}_{x+1}}}^{{{j}^{\prime }},j^{\prime\prime} }X_{{{j}^{\prime }},{{k}_{x}}}^{[x]{{w}_{x-1}}{{w}_{x}}}X_{j^{\prime\prime} ,{{k}_{x+1}}}^{[x+1]{{w}_{x}}{{w}_{x+1}}},\end{eqnarray} \tag{ 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}^{\prime }},j^{\prime\prime} )$ and in output $({{j}_{x}},{{j}_{x+1}})$ must satisfy the link constraint, i.e. the triplets $({{j}^{\prime }},j^{\prime\prime} ,q)$, and $({{j}_{x}},{{j}_{x+1}},{{q}^{\prime }})$ must belong to ${{\Omega }_{x,x+1}}$ for some q and ${{q}^{\prime }}$. 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

  $$\begin{eqnarray}&&{\rm cost}\sim {{d}^{2}}{{b}^{2}}\;{{m}^{3}}+{{d}^{4}}{{b}^{2}}\;{{m}^{2}}\qquad {\rm reduced}\;{\rm to}\to \qquad {\rm cost}\sim \chi \;{{b}^{2}}\;{{m}^{3}}+{{\chi }^{2}}{{b}^{2}}\;{{m}^{2}}\end{eqnarray} \tag{ 37 }$$

  .
• (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

  $$\begin{eqnarray}&&\Gamma _{{{w}_{x+1}},{{k}_{x+1}},{{j}_{x+1}}}^{{{w}_{x-1}},{{k}_{x}},{{j}_{x}}}\simeq \sum xmwY_{{{j}_{x}},{{k}_{x}}}^{[x]{{w}_{x-1}}{{w}_{x}}}Y_{{{j}_{x+1}},{{k}_{x+1}}}^{[x+1]{{w}_{x}}{{w}_{x+1}}}.\end{eqnarray} \tag{ 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³, 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 Γ 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+1}}}^{{{w}_{x-1}},{{k}_{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 Γ 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 Γ 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 $(b\;m\;\xi (q)\times b\;m\;\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 Γ matrix, we perform separate SVDs on each distinct q block. This approach reduces the computational costs as follows

  $$\begin{eqnarray}&&{\rm cost}\sim {{d}^{3}}{{b}^{3}}\;{{m}^{3}}\qquad {\rm reduced}\;{\rm to}\to \qquad {\rm cost}\sim {{b}^{3}}\;{{m}^{3}}\mathop{\sum }\limits_{q}^{{{{\bar{N}}}_{x}}}{{\xi }^{3}}(q).\end{eqnarray} \tag{ 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.

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{eqnarray}&&\begin{array}{l} |1{{\rangle }_{r}}=|0,1,1{{\rangle }_{o}}\quad |2{{\rangle }_{r}}=|1,0,1{{\rangle }_{o}}\quad |3{{\rangle }_{r}}=|1,1,0{{\rangle }_{o}}\quad {\rm{on}} \;{\rm{odd}} \;{\rm{sites}} \\ |1{{\rangle }_{r}}=|1,0,0{{\rangle }_{e}}\quad |2{{\rangle }_{r}}=|0,1,0{{\rangle }_{e}}\quad |3{{\rangle }_{r}}=|0,0,1{{\rangle }_{e}}\quad {\rm{on}} \;{\rm{even}} \;{\rm{sites}} {\rm{.}} \\ \end{array}\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{{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}},\end{eqnarray} \tag{ 42 }$$

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

$$\begin{eqnarray}&&{{V}_{j,q}}=({{\delta }_{j,1}}+{{\delta }_{j,2}}){{\delta }_{q,1}}+{{\delta }_{j,3}}\;{{\delta }_{q,2}}\quad {\rm and}\quad {{Z}_{q,j}}={{\delta }_{q,1}}({{\delta }_{j,2}}+{{\delta }_{j,3}})+{{\delta }_{q,2}}\;{{\delta }_{j,1}},\end{eqnarray} \tag{ 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 ℓ 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 φ being the golden ratio $\varphi =(1+\sqrt{5})/2$. This tells us that for large sizes ℓ, 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\simeq 1.618$, i.e. the scaling is somewhat smoother.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&|2k+1{{\rangle }_{r}}=|k,1,\bar{N}-k{{\rangle }_{o}},\quad |2k{{\rangle }_{r}}=|k,0,\bar{N}+1-k{{\rangle }_{o}}\end{eqnarray} \tag{ 44 }$$

on odd sites, while on even sites

$$\begin{eqnarray}&&|2k+1{{\rangle }_{r}}=|\bar{N}-k,0,k{{\rangle }_{e}},\quad |2k{{\rangle }_{r}}=|\bar{N}-k,1,k-1{{\rangle }_{e}},\end{eqnarray} \tag{ 45 }$$

where $|j{{\rangle }_{r}}$ spans within $1\leqslant j\leqslant 2\bar{N}+1$. Like in the previous example ${{Q}_{r,x,x+1}}$ is homogeneous and reads

$$\begin{eqnarray}\begin{array}{rcl} {{Q}_{r,x,x+1}} & = & |2\bar{N}+1,1\rangle \langle 2\bar{N}+1,1{{|}_{r}}+\mathop{\sum }\limits_{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}}) \\ \end{array}\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{{q}_{\ell }}=k\mathop{\to }\limits^{|2k+1\rangle }{{q}_{\ell +1}}=\bar{N}-k\quad {\rm and}\quad {{q}_{\ell }}=k\mathop{\to }\limits^{|2k\rangle }{{q}_{\ell +1}}=\bar{N}+1-k.\end{eqnarray} \tag{ 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 ℓ, 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 ℓ. 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 α as a function of $\bar{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 $\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-$\frac{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-$\frac{1}{2}$ model with constraints, and these constraints vanish in the Wilson limit, where one should recover the corresponding unconstrained spin-$\frac{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-$\frac{1}{2}$ model [20].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this example, which corresponds to the scenario we already introduced in section 2.0.4, both matter and link fields host spin-$\frac{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 $\bar{N}$ of rishons on a lattice bond in this example to $\bar{N}=1$. With the aforementioned choices, the reduced local space is four-dimensional:

$$\begin{eqnarray}\begin{array}{rcl} |1{{\rangle }_{r}} & = & \frac{1}{\sqrt{2}}\left( |0,\uparrow ,\downarrow \rangle -|0,\downarrow ,\uparrow \rangle \right),\qquad |2{{\rangle }_{r}}=\frac{1}{\sqrt{2}}\left( |\uparrow ,0,\downarrow \rangle -|\downarrow ,0,\uparrow \rangle \right), \\ |3{{\rangle }_{r}} & = & \frac{1}{\sqrt{2}}\left( |\uparrow ,\downarrow ,0\rangle -|\downarrow ,\uparrow ,0\rangle \right),\qquad |4{{\rangle }_{r}}=|0,\phi ,0\rangle , \\ \end{array}\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&\begin{array}{l} {{Q}_{r,x,x+1}}=\left( |1\rangle \langle 1|+|2\rangle \langle 2| \right)\otimes \left( |1\rangle \langle 1|+|4\rangle \langle 4| \right){\mkern 1mu} + \\ \quad \left( |3\rangle \langle 3|+|4\rangle \langle 4| \right)\otimes \left( |2\rangle \langle 2|+|3\rangle \langle 3| \right). \\ \end{array}\end{eqnarray} \tag{ 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 $\bar{D}(\ell )={{2}^{\ell +1}}$ and ultimately $\alpha =2$. This again implies that a mapping of this class of models to a spin-$\frac{1}{2}$ system might exist, even though it has not presently been proved to the best of our knowledge.

Zoom In Zoom Out Reset image size

The last scenario we discuss explicitly is again the $U(2)$ case with spin-$\frac{1}{2}$ particles, but now we set $\bar{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

$$\begin{eqnarray}\begin{array}{rcl} |1{{\rangle }_{r}} & = & |\phi ,\phi ,0{{\rangle }_{o}},\quad |2{{\rangle }_{r}}=\ \frac{1}{\sqrt{2}}\left( |\phi ,\uparrow ,\downarrow {{\rangle }_{o}}-|\phi ,\downarrow ,\uparrow {{\rangle }_{o}} \right), \\ |3{{\rangle }_{r}} & = & \ \frac{1}{\sqrt{2}}\left( |\uparrow ,\phi ,\downarrow {{\rangle }_{o}}-|\downarrow ,\phi ,\uparrow {{\rangle }_{o}} \right),\quad |4{{\rangle }_{r}}=\ \frac{1}{\sqrt{2}}\left( |\uparrow ,\downarrow ,\phi {{\rangle }_{o}}-|\downarrow ,\uparrow ,\phi {{\rangle }_{o}} \right), \\ |5{{\rangle }_{r}} & = & \ |\phi ,0,\phi {{\rangle }_{o}},\quad |6{{\rangle }_{r}}=\ |0,\phi ,\phi {{\rangle }_{o}}, \\ \end{array}\end{eqnarray} \tag{ 50 }$$

while on even sites

$$\begin{eqnarray}\begin{array}{rcl} |1{{\rangle }_{r}} & = & |0,0,\phi {{\rangle }_{e}},\quad |2{{\rangle }_{r}}=\ \frac{1}{\sqrt{2}}\left( |0,\uparrow ,\downarrow {{\rangle }_{e}}-|0,\downarrow ,\uparrow {{\rangle }_{e}} \right), \\ |3{{\rangle }_{r}} & = & \ \frac{1}{\sqrt{2}}\left( |\uparrow ,0,\downarrow {{\rangle }_{e}}-|\downarrow ,0,\uparrow {{\rangle }_{e}} \right),\quad |4{{\rangle }_{r}}=\ \frac{1}{\sqrt{2}}\left( |\uparrow ,\downarrow ,0{{\rangle }_{e}}-|\downarrow ,\uparrow ,0{{\rangle }_{e}} \right), \\ |5{{\rangle }_{r}} & = & \ |0,\phi ,0{{\rangle }_{e}},\quad |6{{\rangle }_{r}}=\ |\phi ,0,0{{\rangle }_{e}}. \\ \end{array}\end{eqnarray} \tag{ 51 }$$

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\lt {{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 \simeq 2.2470$.