††* These authors contributed equally to this work.

# Scalable digital quantum simulation of lattice fermion theories with local
encoding

Marco Ballarin  Dipartimento di Fisica e Astronomia “G. Galilei”, Università
di Padova, I-35131 Padova, Italy.  Padua Quantum Technologies Research Center,
Università degli Studi di Padova  Istituto Nazionale di Fisica Nucleare
(INFN), Sezione di Padova, I-35131 Padova, Italy.  Giovanni Cataldi
Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova,
I-35131 Padova, Italy.  Padua Quantum Technologies Research Center, Università
degli Studi di Padova  Istituto Nazionale di Fisica Nucleare (INFN), Sezione
di Padova, I-35131 Padova, Italy.  Giuseppe Magnifico  Dipartimento di Fisica
e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova, Italy.  Padua
Quantum Technologies Research Center, Università degli Studi di Padova
Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131
Padova, Italy.  Dipartimento di Fisica, Università di Bari, I-70126 Bari,
Italy.  Daniel Jaschke  Dipartimento di Fisica e Astronomia “G. Galilei”,
Università di Padova, I-35131 Padova, Italy.  Padua Quantum Technologies
Research Center, Università degli Studi di Padova  Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy.  Institute
for Complex Quantum Systems, Ulm University, Albert-Einstein-Allee 11, 89069
Ulm, Germany  Ilaria Siloi  Dipartimento di Fisica e Astronomia “G. Galilei”,
Università di Padova, I-35131 Padova, Italy.  Padua Quantum Technologies
Research Center, Università degli Studi di Padova  Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy.  Simone
Montangero  Dipartimento di Fisica e Astronomia “G. Galilei”, Università di
Padova, I-35131 Padova, Italy.  Padua Quantum Technologies Research Center,
Università degli Studi di Padova  Istituto Nazionale di Fisica Nucleare
(INFN), Sezione di Padova, I-35131 Padova, Italy.  Pietro Silvi  Dipartimento
di Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova,
Italy.  Padua Quantum Technologies Research Center, Università degli Studi di
Padova  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova,
I-35131 Padova, Italy.

(August 31, 2024)

###### Abstract

We numerically analyze the feasibility of a platform-neutral, general strategy
to perform quantum simulations of fermionic lattice field theories under open
boundary conditions. The digital quantum simulator requires solely one- and
two-qubit gates and is scalable since integrating each Hamiltonian term
requires a finite (non-scaling) cost. The exact local fermion encoding we
adopt relies on auxiliary $\mathbb{Z}_{2}$ lattice gauge fields by adding a
pure gauge Hamiltonian term akin to the Toric Code. By numerically emulating
the quantum simulator real-time dynamics, we observe spin-charge separation
phenomena for a spin-$\frac{1}{2}$ Hubbard ladder. Additionally, we show that
such local fermion encoding is also useful in tensor network numerical
simulations.

Figure 1: Simulation of the Hubbard model in 2D with tensor networks
(equilibrium) and quantum circuits (out-of-equilibrium). (a) Schematic of the
tree-tensor network (TTN) installed on the 2D lattice: the fermionic degrees
of freedom are removed with the gauge defermionization and encoded in a
dressed site (dashed closed line) comprising fermionic matter and rishons.
With TTN, the ground state of the defermionized model is computed for up to
$4\times 4$ lattice. (b) To run the digital quantum simulation, dressed sites
are decomposed into qubits and Hamiltonian terms in Pauli strings. The
symmetries of the system are here enforced through the stabilizer formalism.
(c) (left) Quantum circuit for simulating out-of-equilibrium dynamics:
adiabatic preparation of the initial state (half-filling, repulsive $U$),
injection of charge (spin) excitation, and time evolution with Hubbard
Hamiltonian. (right) Spin-charge separation schematic: the spin excitation
propagates faster than the charge one.

The quest for quantum simulation of interacting fermionic models Fradkin
(2013); Auerbach (1994); Giamarchi (2003) is considered a necessity to reach a
novel understanding of collective phenomena both at low and high energies Lee
_et al._ (2006); Hartke  _et al._ (2023); Halperin and Jain (2020); Hémery
_et al._ (2023); Hofrichter  _et al._ (2016), but it is hindered by
fundamental challenges Gattringer and Langfeld (2016); Troyer and Wiese
(2005); Altman  _et al._ (2021). While analog quantum simulation, e.g. with
optical lattices, has advanced greatly in recent decades Esslinger (2010);
Bloch  _et al._ (2012); Tarruell and Sanchez-Palencia (2018); Bohrdt  _et al._
(2021); Cheuk  _et al._ (2015); Duarte  _et al._ (2015); Edge  _et al._
(2015); Greif  _et al._ (2013); Haller  _et al._ (2015); Hart  _et al._
(2015); Hofstetter  _et al._ (2002); Jördens  _et al._ (2008); Messer  _et
al._ (2015); Murmann  _et al._ (2015); Omran  _et al._ (2015); Parsons  _et
al._ (2015); Schneider  _et al._ (2008); Taie  _et al._ (2012); Uehlinger  _et
al._ (2013), it presents limitations in tailoring exotic interactions. And
while fermionic digital quantum processors are still at an early stage of
development Zache  _et al._ (2023); González-Cuadra  _et al._ (2023), the
well-established conventional digital quantum simulation platforms (e.g.
superconducting qubits, trapped ions, Rydberg arrays, quantum dots) Arute  _et
al._ (2020); Barends  _et al._ (2015); Salathé  _et al._ (2015); O’Malley  _et
al._ (2016); Stanisic  _et al._ (2022) are built on distinguishable, spatially
localized qubits (or qudits). In this framework, ‘Fermion Encoding’ is the
analytical process of exactly converting a fermionic algebra (mutually-
anticommuting operations) into a genuinely local algebra (mutually-commuting
operations) of qudits.

Such encoding can not be carried out free of cost. Traditional strategies
focused on encoding $N$ Dirac orbitals into $N$ qubits. By construction, these
strategies can not preserve locality and, as a result, end up encoding the
ubiquitous two-body interactions into cumbersome $W$-body interactions (with
$W$ often addressed as _Pauli weight_). Alongside recent efforts, which were
able to reduce $W$ from linear to logarithmic in $N$ Jordan and Wigner (1928);
Bravyi and Kitaev (2002), a separate sector of strategies arose: encodings
attempting to preserve locality Fradkin  _et al._ (1980); Srednicki (1980);
Verstraete and Cirac (2005); Chen  _et al._ (2018). They exhibit a flat Pauli
weight ($W$ does not scale with $N$) at the price of requiring a number of
qubits larger (but still linear scaling) than the number $N$ of fermion
orbitals. For non-gauge lattice fermion theories, Ref. Kitaev (2006) showed
that it is sufficient to add one qubit for each lattice bond to achieve local
encoding (a generalization to lattice gauge fermionic theories was also
recently introduced Zohar and Cirac (2018, 2019)). The extra qubits play the
role of an effective, discrete lattice gauge field with pure gauge dynamics
akin to the Toric Code Hamiltonian Kitaev (2003a). The formal mapping is well
understood Chen  _et al._ (2023), but evidence of practical feasibility of the
local encoding for numerical simulation or digital quantum simulation is still
subject of active research Pardo  _et al._ (2023); Irmejs  _et al._ (2023a,
b).

In this work, we perform a twofold study on the feasibility of the local
fermion encoding for classical and quantum simulation, as summarized in Fig.
1: on one hand side, we investigate its performance for ground-state numerical
simulations based on tensor network ansatz states. On the other hand, we
numerically test its applicability and scalability for digital quantum
simulation on a generic platform (considering standard 1-qubit and 2-qubit
programmable dynamical resources), by numerically emulating the quantum
simulation processing in real time. Equipped with these methods, we
investigate the properties of the spin-$\frac{1}{2}$ Hubbard model on a two-
dimensional square lattice, both at equilibrium and out-of-equilibrium, via
tensor networks. At zero temperature, we can access system sizes that allow us
to identify the transition between liquid and insulating (crystalline spin-
lattice) phases. In real-time dynamics on a ladder, we observe a sharp
signature of spin-charge separation, an effect that has been observed in 1D
Hubbard chains Arute  _et al._ (2020); Vijayan  _et al._ (2020).

The paper is organized as follows: in Sec. I, we revisit the gauge-based
encoding to locally remove the fermionic matter, or gauge defermionization. We
devote Sec. II to the implementation of the encoding on a tree-tensor network
(TTN) ansatz and show ground state simulation results of the 2D Hubbard model.
Finally, in Sec. III, we detail our prescription for digital quantum
simulation of the 2D Hubbard model, based on gauge defermionization, and
numerically emulate the digital quantum simulation to observe spin-charge
separation effects in two-dimensional lattice samples.

##  I Gauge Defermionization

In this section, we revisit and extend a known technique, based on lattice
gauge theories, to eliminate fermionic matter from 2D lattice models Zohar and
Cirac (2018, 2019). The strategy consists of developing an analytical mapping
between an input fermion lattice Hamiltonian (here we consider it to be a pure
non-gauge fermion theory) and a lattice gauge Hamiltonian, with
$\mathbb{Z_{2}}$ gauge symmetry, equipped with appropriate gauge constraints.
Furthermore, the resulting gauge theory can be manipulated so that the
fermionic parity at each (dressed) site is protected, thus resulting in a
theory where operator algebras at different sites always commute, i.e. a local
spin theory.

Our formulation of the defermionization technique can be applied to any 2D
fermion lattice Hamiltonian, regardless of the lattice system (it works on
Honeycomb and Kagome lattices as well as irregular lattices), as long as few
general constraints are satisfied: _(i)_ the global parity of fermions is
protected; _(ii)_ Hamiltonian terms that flip the local fermion parity must be
nearest-neighbor; _(iii)_ the full system must be under open boundary
conditions, i.e. every closed path on the lattice must be topologically
shrinkable to a point. These conditions are often satisfied in solid state
Hamiltonians and cold atoms systems Lieb (2013).

Therefore, without any loss of generality, we consider a finite 2D lattice
$\Lambda=(L_{x},L_{y})$ whose sites and links are respectively identified by
the couple ($\mathbf{j},{\bf\it\mu}$), where $\mathbf{j}=(j_{x},j_{y})$ is any
2D site, while ${\bf\it\mu}$ is one of the two positive lattice unit vectors:
${\bf\it\mu}_{x}=(1,0)$, ${\bf\it\mu}_{y}=(0,1)$. Then, a general fermion
lattice Hamiltonian reads:

$\begin{split}H_{0}=&\sum_{\mathbf{j}\in\Lambda}\sum_{{\bf\it\mu}}\sum_{f,f^{\prime}}\quantity(h_{f,f^{\prime}}\,\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}+\text{H.c.}).\\\
&+\sum_{\mathbf{j}\in\Lambda}\sum_{{\bf\it\mu}}\sum_{f,f^{\prime}}\quantity(\Delta_{f,f^{\prime}}\,\psi_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}+\mbox{H.c.}\vphantom{h_{f,f^{\prime}}\psi^{\dagger}_{\mathbf{j},f}})\\\
&+V\left(\\{n_{\mathbf{j},f}\\}\right),\end{split}$ (1)

where the fermion fields satisfy the usual Dirac anti-commutation relations

$\displaystyle\quantity{\psi^{\dagger}_{\mathbf{j},f},\psi_{\mathbf{j^{\prime}},f^{\prime}}}$
$\displaystyle=\delta_{\mathbf{j},\mathbf{j^{\prime}}}\delta_{f,f^{\prime}}$
$\displaystyle\quantity{\psi_{\mathbf{j},f},\psi_{\mathbf{j^{\prime}},f^{\prime}}}$
$\displaystyle=0,$ (2)

and can exhibit an internal degree of freedom, here labeled by the flavor
index $f$.

The interaction potential $V$ term, which solely is a function of the fermion
densities $n_{\mathbf{j},f}=\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f}$,
can in principle have any range and shape, and even include chemical
potentials. Conversely, the hopping $h_{f,f^{\prime}}$ and the double
creation/annihilation $\Delta_{f,f^{\prime}}$ processes (the latter often seen
in superconductors) are the terms that break local fermion parity: their
single-site algebras do not commute between distant sites (i.e. they are not
_genuinely_ local), and this is the source of all the numerical difficulty
when simulating lattice fermions.

While in one spatial dimension, the Jordan-Wigner transformation provides an
easy solution to this problem, tackling the mapping in higher dimensions
requires sophisticated, often cumbersome techniques Verstraete and Cirac
(2005); Corboz  _et al._ (2010). In this perspective, converting a fermion
algebra into a genuinely local (spin-like) algebra using introducing a lattice
gauge field is an elegant strategy, which is also practical from a numerical
simulation standpoint, as we discuss later on.

###  I.1 Step 1: Mapping into a $\mathbb{Z_{2}}$ gauge theory

We now perform a set of exact algebraic manipulations to the Hamiltonian of
Eq. 1 until we reach a defermionized form. The very first step is to promote
the total fermion parity operator at site $\mathbf{j}$ to a gauge
transformation (for the matter sites $M$), namely

$G_{\mathbf{j}}^{[M]}=\text{exp}\bigg{(}i\pi{\textstyle{\sum_{f}}}\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f}\bigg{)}.$
(3)

It is indeed a parity operator, or $\mathbb{Z}_{2}$ group, since

$G_{\mathbf{j}}^{[M]}=G_{\mathbf{j}}^{[M]\dagger}=(G_{\mathbf{j}}^{[M]})^{-1}.$
(4)

Under $G_{\mathbf{j}}^{[M]}$, operators that preserve total fermion parity are
left invariant, such as densities
$G_{\mathbf{j}}^{[M]}n_{\mathbf{j},f}G_{\mathbf{j}}^{[M]}=n_{\mathbf{j},f}$.
Conversely, operators that flip the total fermion parity acquire a sign, that
is
$G_{\mathbf{j}}^{[M]}\psi_{\mathbf{j},f}G_{\mathbf{j}}^{[M]}=-\psi_{\mathbf{j},f}$.

At this stage, we promote the gauge transformation into a gauge symmetry. This
task is performed by adding an auxiliary quantum lattice field, the gauge
field, on the _bonds_ $(\mathbf{j},\mathbf{j}+{\bf\it\mu})$ of the lattice.
The local Hilbert space for a gauge site should correspond to the regular
representation space for the $\mathbb{Z}_{2}$ Zohar and Burrello (2015), so it
should be a two-level system, or a qubit, equipped with the (genuinely local)
algebra of Pauli matrices ${\vec{\sigma}}$.

Since the left- $G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[L]}$ and right-
$G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[R]}$ groups of transformations of the
gauge field must commute and square to the identity, without loss of
generality, we can set them both to
$G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[L]}=G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[R]}=\sigma^{z}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}$.

Thanks to this auxiliary field, we equip the gauge-violating terms of the
Hamiltonian with a parallel transporter operator. Namely, we replace

$\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}\quad\longrightarrow\quad\psi^{\dagger}_{\mathbf{j},f}U_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}},$
(5)

that transforms covariantly under the gauge field groups:

$G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[L]}U_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[L]}{=}G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[R]}U_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[R]}{=}-U_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}$
(6)

Again, without loss of generality, we set
$U_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}=\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}$.

Overall, this procedure translates into modifying the fermion parity-flipping
terms from Eq. 1 according to

$\displaystyle
h_{f,f^{\prime}}\,\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}$
$\displaystyle\longrightarrow\quad
h_{f,f^{\prime}}\,\psi^{\dagger}_{\mathbf{j},f}\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}$
(7)
$\displaystyle\Delta_{f,f^{\prime}}\,\psi_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}$
$\displaystyle\longrightarrow\quad\Delta_{f,f^{\prime}}\,\psi_{\mathbf{j},f}\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}\,.$

With this modification, the Hamiltonian satisfies a gauge symmetry around
every site, precisely
$G_{\mathbf{j}}^{[\text{total}]}=G_{\mathbf{j}}^{[M]}\prod_{\mu}G_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}^{[L/R]}$.
Consequently, we have to choose a symmetry sector (for each of these
symmetries) to play the role of physical space. As it is common for lattice
gauge theories, we consider the quantum states _invariant_ under all the
$G_{\mathbf{j}}^{[\text{total}]}$, i.e. which satisfy

$\left[\exp\left(i\pi{\textstyle{\sum_{f}}}\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f}\right)\prod_{{\bf\it\mu}}\sigma^{z}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\right]\ket{\Psi_{\text{phys}}}{=}\ket{\Psi_{\text{phys}}},$
(8)

in order to be the physical states $\ket{\Psi_{\text{phys}}}$. This equation
plays the role of effective Gauss’ Law of the resulting gauge theory, and the
condition can be satisfied at every site _only_ if the total number of
fermions in the system is even. To simulate an odd number of fermions, one or
more virtual bonds, going out of the system, must be added to keep track of
the appropriate parity gauge flux.

###  I.2 Step 2: Restoring equivalence with Plaquette Operators

It is then possible to check that, under the Gauss’ Law of Eq. 8, the gauge-
theory modification Eq. 7 does not change the actual dynamics given by the
original Hamiltonian as in Eq. 1, _but only_ when the lattice is a non-cyclic
graph (such as a 1D chain, or a dendrimer lattice). For each (product state)
fermion configuration, the extra degrees of freedom introduced by the gauge
fields are completely locked by Gauss’ Law; thus all residual degeneracy is
removed, and all the matrix elements of $H_{0}$ are unaltered.

Conversely, when lattice cycles are present, the modification changes the
dynamics. In fact, for each elementary closed cycle of sites, or _plaquette_
$\square\in\Lambda$, one Gauss’ Law operator becomes linearly dependent on
others, thus contributing to degeneracy with a 2-fold space. Moreover, when a
fermion winds around a closed cycle, all of the gauge fields in its path
undergo a $\sigma^{x}$ flip, thus the final state may differ from the original
one. If we want to restore the dynamics of the original model, we have to add
more physical contents to the plaquette, either in terms of an additional
_pure gauge field_ Hamiltonian or in terms of _extra gauge symmetries_.

Fortunately, a simple recipe to do either is inspired by Kitaev’s Toric Code
Kitaev (2003b). The idea is to add a plaquette Hamiltonian term with a
$\sigma^{x}$ for each gauge link of the plaquette, that is

$\begin{split}H_{p}&=-\sum_{\square\in\Lambda}\prod_{\langle\mathbf{j},\mathbf{j^{\prime}}\rangle\in\square}\sigma^{x}_{\mathbf{j},\mathbf{j^{\prime}}}=-\sum_{\square\in\Lambda}\quantity(\begin{array}[]{ccc}\ulcorner\hfil&\sigma^{x}&\urcorner\\\
\sigma^{x}&&\sigma^{x}\\\ \llcorner&\sigma^{x}&\lrcorner\\\
\end{array}).\end{split}$ (9)

Notably, each of these plaquette operators commutes with the lattice gauge
Hamiltonian and can be equivalently cast as a (gauge) symmetry. Moreover,
plaquettes also commute with Gauss’ laws, because a vertex and a plaquette
always share either zero or two bonds, regardless of the lattice system, and
$[\sigma^{x}_{1}\sigma^{x}_{2},\sigma^{z}_{1}\sigma^{z}_{2}]=0$. Ultimately,
each plaquette operator completely resolves the 2-fold degeneracy introduced
by the plaquette itself, and each closed path of moving fermions returns to
its original state with the correct amplitude and phase.

In the end, we have mapped the original lattice fermion Hamiltonian, Eq. 1, to
a dynamically-equivalent, $\mathbb{Z}_{2}$-invariant lattice gauge model,
which reads

$\begin{split}H_{1}&=\sum_{\mathbf{j},{\bf\it\mu}}\sum_{f,f^{\prime}}\left(h_{f,f^{\prime}}\,\psi^{\dagger}_{\mathbf{j},f}\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}\right.\\\
&\left.+\Delta_{f,f^{\prime}}\,\psi_{\mathbf{j},f}\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}+\mbox{H.c.}\vphantom{h_{f,f^{\prime}}\psi^{\dagger}_{\mathbf{j},f}}\right)+V\left(\\{n_{\mathbf{j},f}\\}\right),\end{split}$
(10)

with the constraints

$\left\\{\begin{aligned}
&\left(e^{i\pi\sum_{f}\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f}}\prod_{\mu}^{\text{vertex}}\sigma^{z}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\right)\ket{\Psi_{\text{phys}}}=\ket{\Psi_{\text{phys}}},\\\
&\left(\prod_{\left\langle\mathbf{j},\mathbf{j^{\prime}}\right\rangle}^{\text{plaquette}}\sigma^{x}_{\mathbf{j},\mathbf{j^{\prime}}}\right)\ket{\Psi_{\text{phys}}}=\ket{\Psi_{\text{phys}}}.\end{aligned}\right.$
(11)

A rigorous proof of the equivalence between Eq. 1 and Eqs. (10)-(11) was
provided for a number-conserving theory on the square lattice in Zohar and
Cirac (2018); the generalization to other cases is fairly straightforward.

Notice that the degeneracy is not fully removed in periodic boundary
conditions (PBC), where each winding dimension introduces an additional closed
cycle. This case requires the additional constraint of a full $\sigma^{x}$
string along each winding dimension, but we will not treat this case here
because the theory becomes non-local.

###  I.3 Step 3: Fermionic Rishons

So far, it seems that we have increased the complexity of the fermion lattice
theory. In this last step, we will locally manipulate the gauge fields and
achieve defermionization explicitly, ending with an algebra of genuinely local
operators only. Carrying out this last step is advantageous for tensor
networks, but not so much for the digital quantum simulator (see further
discussion in Sec. II.1 and Sec. III).

To do so, we split each gauge field, living on the
$(\mathbf{j},\mathbf{j}+{\bf\it\mu})$ lattice link, in a pair of spinless
fermion modes (rishons) $c_{\mathbf{j},{\bf\it\mu}}^{(\dagger)}$ and
$c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}^{(\dagger)}$, equipped with another
symmetry. These new Dirac fermion operators belong to dressed sites
$\mathbf{j}$ and $\mathbf{j}+{\bf\it\mu}$ respectively and satisfy the usual
anti-commutation relations
$\\{c^{\dagger}_{\mathbf{j},{\bf\it\mu}},c_{\mathbf{j}^{\prime},{\bf\it\mu}^{\prime}}\\}=\delta_{\mathbf{j}j^{\prime}}\delta_{\mu\mu^{\prime}}$
and
$\\{c_{\mathbf{j},{\bf\it\mu}},c_{\mathbf{j}^{\prime},{\bf\it\mu}^{\prime}}\\}=0$.
The combined 4-dimensional space of the two modes is then reduced back to the
2-dimensional space of a qubit by imposing that the total number of rishon
fermions on a bond must be an even number. This constraint can again be cast
as a gauge symmetry, requiring that

$\exp[i\pi(c_{\mathbf{j},{\bf\it\mu}}^{\dagger}c_{\mathbf{j},{\bf\it\mu}}{+}c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}^{\dagger}c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}})]{\ket{\Psi_{\text{phys}}}}{=}\ket{\Psi_{\text{phys}}}.$
(12)

We can now convert the Pauli algebra of the gauge fields into an operator
algebra acting on the Rishon pair:

$\displaystyle\sigma^{z}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}$
$\displaystyle\to
1-2c_{\mathbf{j},{\bf\it\mu}}^{\dagger}c_{\mathbf{j},{\bf\it\mu}}=1-2c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}^{\dagger}c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}$
(13) $\displaystyle\sigma^{x}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}$
$\displaystyle\to
i\gamma_{\mathbf{j},{\bf\it\mu}}\gamma_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}},$

where now
$\gamma_{\mathbf{j},{\bf\it\mu}}{=}c_{\mathbf{j},{\bf\it\mu}}{+}c^{\dagger}_{\mathbf{j},{\bf\it\mu}}{=}\gamma_{\mathbf{j},{\bf\it\mu}}^{\dagger}$
is a Majorana operator on the rishon fermion, which squares to the identity
$\gamma_{\mathbf{j},{\bf\it\mu}}^{2}{=}\mathbb{1}$ but still anti-commutes
with other fermions
$\\{\gamma_{\mathbf{j},{\bf\it\mu}},\gamma_{\mathbf{j}^{\prime},{\bf\it\mu}^{\prime}}\\}{=}2\delta_{\mathbf{j},\mathbf{j^{\prime}}}\delta_{{\bf\it\mu},{\bf\it\mu^{\prime}}}$.
These operators defined at Eq. 13 respect the fermion parity symmetry on the
bond, and on the even parity sector, they act exactly as Pauli matrices.

We can then plug this exact manipulation into the gauge theory Hamiltonian,
resulting in

$\begin{split}H&{}_{2}=\sum_{\mathbf{j},{\bf\it\mu}}\sum_{f,f^{\prime}}\quantity(ih_{f,f^{\prime}}\,\psi^{\dagger}_{\mathbf{j},f}\gamma_{\mathbf{j},{\bf\it\mu}}\gamma_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}+\text{H.c.})\\\
&+\sum_{\mathbf{j},{\bf\it\mu}}\sum_{f,f^{\prime}}\quantity(i\Delta_{f,f^{\prime}}\,\psi_{\mathbf{j},f}\gamma_{\mathbf{j},{\bf\it\mu}}\gamma_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f^{\prime}}+\mbox{H.c.}\vphantom{h_{f,f^{\prime}}\psi^{\dagger}_{\mathbf{j},f}})\\\
&+V\quantity(\\{n_{\mathbf{j},f}\\}),\end{split}$ (14)

where we have to add the new bond parity symmetry to the constraints, thus

$\left\\{\begin{aligned}
&\left[e^{i\pi\sum_{f}\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f}}\prod_{\mu}^{\text{vertex}}e^{i\pi
c^{\dagger}_{\mathbf{j},{\bf\it\mu}}c_{\mathbf{j},{\bf\it\mu}}}\sigma^{z}_{\mathbf{j},\mathbf{j}+{\bf\it\mu}}\right]\ket{\Psi_{\text{phys}}}=\ket{\Psi_{\text{phys}}}\\\
&\left(\prod_{\mathbf{j},{\bf\it\mu}}^{\text{plaquette}}i\gamma_{\mathbf{j},{\bf\it\mu}}\gamma_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}\right)\ket{\Psi_{\text{phys}}}=\ket{\Psi_{\text{phys}}}\\\
&e^{i\pi c_{\mathbf{j},\mu}^{\dagger}c_{\mathbf{j},\mu}}e^{i\pi
c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}^{\dagger}c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}}\ket{\Psi_{\text{phys}}}=\ket{\Psi_{\text{phys}}},\end{aligned}\right.$
(15)

which must hold for every vertex, every plaquette, and every bond
respectively.

We have finally reached the final form of our model. Indeed now we have that
every term of the Hamiltonian and every gauge symmetry _protects total fermion
parity_ at each dressed site, counting together matter and rishon fermions
(for plaquettes, remember that each dressed site contributes with two rishon
modes in a closed lattice path). Therefore, it is possible to think of each
dressed site as a large spin, and both Hamiltonian and gauge symmetries can be
written in terms of _genuinely local_ operator algebras, which commute on
different dressed sites. The theory has been effectively defermionized.

###  I.4 The price of defermionization

Defermionization of the lattice model using lattice gauge theory, i.e. going
from Eq. 1 to Eqs. (14)-(15) is an exact mapping. And while it provides the
clear benefits of eliminating fermionic operator algebras, it does not come
free of costs.

First of all, the transformation increases the local space dimension, only to
shrink the local dimension again once the gauge symmetries are installed. If
the model has $N$ sites, $f$ flavors and coordination number $v$, we increase
the dimension from $2^{fN}$ to $2^{(f+v)N}$.

Secondly, we impose a new interaction term in the form of the plaquettes.
While for typical 2D lattices plaquettes are rather small, it is still an
interaction involving three to six sites depending on the lattice system. If
the original interaction $V$ was on-site, or nearest-neighbor, the
defermionization effectively increased the interaction supports. Moreover,
there are no more components in the Hamiltonian that are quadratic in the
Fermi operators (except chemical potentials within $V$). This means that
Green’s function perturbative approaches, which start from the free Fermi gas
propagator of the quadratic theory, are no longer viable.

Finally, introducing an auxiliary field whose pure gauge dynamics is analogous
to a Toric Code Hamiltonian has the drawback of increasing the entanglement.
Even for product states of fermions, where fermions are locked in a specific
integer filling configuration, and the original entanglement is zero (e.g. for
a Jordan-Wigner encoding), adding a Toric Code field on top of that raises the
entanglement entropy to an exact area-law Kitaev and Laumann (2009), carrying
one e-bit of entanglement per plaquette that is cut by the bi-partition. This
observation has substantial implications for tensor network simulations of
said models, as we will discuss in detail later on.

##  II Hubbard Model, Defermionized for tensor networks

Despite its Hamiltonian’s apparent simplicity, the Hubbard model has eluded
physicists for decades Hubbard (1963); Arovas  _et al._ (2021). While exact
solutions are available for both one-dimensional Lieb and Wu (1968) and
infinite-dimensional cases Metzner and Vollhardt (1989); Müller-Hartmann
(1989), addressing finite-size systems in higher dimensions at arbitrary
temperature has required the development of various computational techniques.
Although Monte Carlo methods can handle substantially large systems, they
suffer from the well-known sign problem when computing physically significant
quantities within specific regimes. Comprehensive reviews detailing and
comparing the accomplished computational outcomes for the 2D case can be found
in the following references Qin  _et al._ (2021); LeBlanc  _et al._ (2015).

In this section, we join the effort by performing tensor network simulations
of the Hubbard model at equilibrium at zero temperature. Here, we briefly
review the model and manipulate its defermionized formulation to be ready for
tensor network simulation. Subsequently, we will focus on a 2D _square
lattice_ geometry, fermionic matter with 2 flavors (e.g. $\uparrow$ and
$\downarrow$), and open boundary conditions. On a rectangle lattice $\Lambda$
of $|\Lambda|=L_{x}\times L_{y}$ sites, the Hubbard Hamiltonian reads:

$\begin{split}H_{\text{Hub}}=&-t\sum_{\mathbf{j}\in\Lambda}\sum_{{\bf\it\mu}}\sum_{f=\uparrow,\downarrow}\quantity[\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j}+{\bf\it\mu},f}+\text{H.c.}]\\\
&+U\sum_{\mathbf{j}\in\Lambda}\quantity(n_{\mathbf{j},\uparrow}-\frac{1}{2})\quantity(n_{\mathbf{j},\downarrow}-\frac{1}{2}),\end{split}$
(16)

where, apart from an energy rescaling, its ground state properties only depend
on the dimensionless parameter $U/t$. The Hubbard model is regarded as the
simplest theory of strongly correlated electrons, where band electrons
interact via a two-body repulsive Coulomb interaction Fradkin (2013). This
model enables the description of a wide range of phenomena including metal-
insulator transitions, superconductivity, and magnetism Giamarchi (2003). The
Hubbard Hamiltonian in Eq. 16 belongs to the Hamiltonian class of Eq. 1, with
hopping terms transparent to flavor
($h_{f,f^{\prime}}=-t\,\delta_{f,f^{\prime}}$) and no pair creation processes
($\Delta_{f,f^{\prime}}=0$), so it can be readily defermionized.

The form of Eq. 16 includes the $-\frac{1}{2}$ terms in the interaction
component, which are equivalent to setting a specific uniform chemical
potential, to explicitly reveal the rich symmetry content of the Hubbard
model. In fact, besides the obvious lattice symmetries (translation by 1 site
in $x$ and $y$, $\frac{\pi}{2}$ rotations, vertical and horizontal reflection,
and compositions thereof) Hubbard dynamics exhibits two useful _glocal_
symmetries, where the transformation is global but comprised of separate
single-site operations, both of them non-Abelian.

The first symmetry, an $SU(2)$ group, represents rotational invariance in the
flavors. Since the hopping is flavor-transparent, and the interaction is based
on double occupancy, the model has to be flavor-invariant. This symmetry is a
Lie group, generated by the operator algebra

$\vec{S}_{\text{tot}}=\sum_{\mathbf{j}\in\Lambda}\vec{S}_{\mathbf{j}},\quad\mbox{where}\quad\vec{S}_{\mathbf{j}}=\frac{1}{2}\sum_{f,f^{\prime}}\vec{\sigma}_{f,f^{\prime}}\psi^{\dagger}_{\mathbf{j},f}\psi_{\mathbf{j},f^{\prime}}$
(17)

are the single-site flavor (spin) operators, acting non-trivially only on the
singly-occupied sites.

The second symmetry, assembling total fermion number conservation and
particle-hole inversion, is an $O(2)$ group, the symmetry group of the circle
(rotations of a scalar angle plus one reflection). Despite being a continuous
group, it is not Lie, but any element can be written as a rotation
${R}_{\text{tot}}(\theta)=\prod_{\mathbf{j}}{R}_{\mathbf{j}}(\theta)$, with
${R}_{\mathbf{j}}(\theta)=e^{i\theta(n_{\mathbf{j},\uparrow}+n_{\mathbf{j},\downarrow})}$,
eventually followed by a reflection
${F}_{\text{tot}}=\prod_{\mathbf{j}}{F}_{\mathbf{j}}$. It is indeed possible
to write the particle-hole transformation $F$ in a way that commutes with the
flavor-rotations $\vec{S}$, thus the two symmetries are indeed independent.
Doing so yields precisely

$\Bigg{\\{}\begin{aligned}
R_{\mathbf{j}}(\theta)\,\psi_{\mathbf{j},f}\,R_{\mathbf{j}}^{\dagger}(\theta)&=e^{-i\theta}\psi_{\mathbf{j},f}\\\
F_{\mathbf{j}}\,\psi_{\mathbf{j},f}\,F_{\mathbf{j}}^{\dagger}&=(-1)^{j}{\textstyle{\sum_{f^{\prime}}}}\sigma^{y}_{f,f^{\prime}}\psi_{\mathbf{j},f^{\prime}}^{\dagger}\end{aligned}.$
(18)

Other symmetries, such as the pseudospin conservation Lieb (1989); Arovas  _et
al._ (2021), are indeed present, but not practically useful for our numerical
simulation purposes.

When defermionizing the Hubbard model for tensor network simulation, we find
it convenient to enforce the _vertex gauge constraint_ as an exact, actual
symmetry. By the prescription discussed in Ref. Silvi  _et al._ (2014), this
constraint will allow us to select a reduced canonical basis for the dressed
site, made by only and all the vertex-gauge invariant states. Conversely, the
_link gauge constraint_ and especially the _plaquette gauge constraint_ are
cumbersome to treat as exact symmetries (their local algebras do not commute)
and instead are included in the Hamiltonian as penalty terms, that increase
the energy of gauge symmetry-violating sectors. In conclusion, we have

$\begin{split}H^{\prime}_{\text{Hub}}&=-t\sum_{f}\sum_{\mathbf{j},{\bf\it\mu}}\quantity[i\psi^{\dagger}_{\mathbf{j},f}\gamma_{\mathbf{j},{\bf\it\mu}}\gamma_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}\psi_{\mathbf{j}+{\bf\it\mu},f}]\\\
&+U\sum_{\mathbf{j}}\left(n_{\mathbf{j},\uparrow}-\frac{1}{2}\right)\left(n_{\mathbf{j},\downarrow}-\frac{1}{2}\right),\end{split}$
(19)

plus a penalty

$\begin{split}H^{\prime}_{\text{pen}}&=+\alpha_{b}\sum_{\mathbf{j},{\bf\it\mu}}\quantity((-1)^{c^{\dagger}_{\mathbf{j},{\bf\it\mu}}c_{\mathbf{j},{\bf\it\mu}}+c^{\dagger}_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}c_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}}-1)\\\
&+\alpha_{p}\sum_{\square\in\Lambda}\quantity(\begin{array}[]{ccc}\ulcorner\hfil&\gamma_{+{\bf\it\mu}_{x}}\gamma_{-{\bf\it\mu}_{x}}&\urcorner\\\
\gamma_{-{\bf\it\mu}_{y}}&&\gamma_{-{\bf\it\mu}_{y}}\\\
\gamma_{+{\bf\it\mu}_{y}}&&\gamma_{+{\bf\it\mu}_{y}}\\\
\llcorner&\gamma_{+{\bf\it\mu}_{x}}\gamma_{-{\bf\it\mu}_{x}}&\lrcorner\\\
\end{array}-1),\end{split}$ (20)

with separate penalty couplings $\alpha_{p}{>}0$ and $\alpha_{b}{>}0$ for the
plaquette and bond violations respectively. The added $-1$ constants ensure
that the correct gauge symmetry sector provides no energy contribution.

Once the vertex gauge symmetry is successfully installed, we are left with a
dressed-site vertex gauge-invariant canonical basis of dimension 32 (4 matter
states and four 2-dimensional rishon modes, divided in half by the vertex
constraint). Within this canonical basis, the algebra of 1-site operators is
genuinely local and can be expressed in terms of the following quadratic
operators

$\displaystyle Q_{\mathbf{j},{\bf\it\mu},f}$
$\displaystyle=\gamma_{\mathbf{j},{\bf\it\mu}}\psi_{\mathbf{j},f}$
$\displaystyle
C_{\mathbf{j},{\bf\it\mu}_{a},{\bf\it\mu}_{b}}=\gamma_{\mathbf{j},{\bf\it\mu}_{a}}\gamma_{\mathbf{j},{\bf\it\mu}_{b}}$
(21) $\displaystyle W_{\mathbf{j},{\bf\it\mu}}$
$\displaystyle=1-2c^{\dagger}_{\mathbf{j},{\bf\it\mu}}c_{\mathbf{j},{\bf\it\mu}},$

each one preserving local fermion parity by design. After this re-formatting,
we have a final expression for the defermionized Hubbard model, and it reads

$\begin{split}H^{\prime\prime}_{\text{Hub}}=&-t\sum_{f}\sum_{\mathbf{j}\in\Lambda}\sum_{{\bf\it\mu}}\quantity[iQ^{\dagger}_{\mathbf{j},{\bf\it\mu},f}Q_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu},f}+\text{H.c.}]\\\
&+U\sum_{\mathbf{j}\in\Lambda}\quantity(n_{\mathbf{j},\uparrow}-\frac{1}{2})\quantity(n_{\mathbf{j},\downarrow}-\frac{1}{2}),\end{split}$
(22)

plus a penalty

$\begin{split}H^{\prime\prime}_{\text{pen}}=&-\alpha_{b}\sum_{\mathbf{j}\in\Lambda}\sum_{{\bf\it\mu}}\quantity(W_{\mathbf{j},{\bf\it\mu}}W_{\mathbf{j}+{\bf\it\mu},-{\bf\it\mu}}-1)\\\
&-\alpha_{p}\sum_{\square\in\Lambda}\quantity(\begin{matrix}C_{\ulcorner}\hfil&C_{\urcorner}\\\
C_{\llcorner}&C_{\lrcorner}\end{matrix}-1).\end{split}$ (23)

Using tensor networks, we simulate the model as it exactly appears in these
expressions.

###  II.1 Numerical Results

Figure 2: Exact Diagonalization comparison on a $2\times 2$ lattice between
the ground state energy density $\varepsilon$ of the original 2D Hubbard model
and its defermionized version as a function of $U/t$, and for three values of
the particle density $\rho$: (a) below half-filling with $\rho=0.5$, (b) at
half-filling with $\rho=1.0$, (c) above half-filling with $\rho=1.5$. Figure
3: (a) Ground state energy density of the 2D Hubbard Hamiltonian at $\rho=1$
for different lattice sizes. Results concerning $2\times 2$, $3\times 2$, and
$3\times 3$ lattices are obtained via exact diagonalization (ED) of the
original Hubbard Hamiltonian in Eq. 16. Energies of the $4\times 4$ lattice
result from TTN simulations of the defermionized Hamiltonian in Eqs. (22)-(23)
with bond dimension $\chi=350$ and energy penalties
$\alpha_{p}=\alpha_{b}=15$. (b) Double occupancy $D$ as a function of $U/t$ is
computed with different numerical methods. Data for MRPHF and DMET are taken
from LeBlanc  _et al._ (2015).

In this section, we present the numerical results obtained using exact
diagonalization (ED) and tree tensor networks (TTN) algorithms Cataldi  _et
al._ (2021); Ferrari  _et al._ (2022); Cataldi  _et al._ (2023). Our numerical
computations concern only open boundary conditions (OBC).

To numerically check the validity of our mapping, we compare the original 2D
Hubbard model of Eq. 16 and its defermionized version described in Eqs.
(22)-(23). For both Hamiltonians, we perform ED on a $2\times 2$ lattice.

Since the two Hamiltonians do not fix a specific number of particles on the
lattice, we add to them an extra term of the form
$H_{N}=\tilde{\mu}\quantity(\sum_{\mathbf{j}\in\Lambda}\quantity(n_{\mathbf{j},\uparrow}+n_{\mathbf{j},\downarrow})-N)^{2}$,
where $\tilde{\mu}$ plays the role of a large penalty coefficient that
increases the energy of all the states with a number of particles differing
from $N$. In this way, by tuning $N$ and setting $\alpha_{p}=\alpha_{b}=20$
and $\tilde{\mu}=20$, we perform ED at fixed values of the particle density
$\rho=N/|\Lambda|$, comparing the ground-state energy densities of the two
models as a function of the ratio $U/t$.

As displayed in Fig. 2 for three values of the particle density, i.e.
$\rho\in\quantity{0.5,1.0,1.5}$, the relative distance between the energies of
the two models ($E_{H}$ for the original Hubbard model, $E_{def.H}$ for its
defermionized version) $\Delta
E=\absolutevalue{E_{H}-E_{def.H}}/\absolutevalue{E_{H}}<10^{-12}$ confirms the
exactness of our mapping. As expected, the above half-filling particles are
forced to form at least a single-site pair which determines a linear
dependence of the energy with $U/t$.

Aware of the equivalence between the two Hamiltonians, we focus on the
defermionized model for larger lattice sizes by using TTN simulations. Our TTN
algorithm for variational ground state search exploits the global $U(1)$
symmetry of the Hubbard model and the Krylov subspace expansion Silvi  _et
al._ (2019). Thus, the additional term $H_{N}$ is no longer needed, as the
algorithm directly conserves the total number of particles by encoding the
$U(1)$ symmetry sectors in the tensor networks.

Figure 4: Ground state characterization for the defermionized Hamiltonian with
$U/t=0.0,4.0$ numerically simulated via TTNs. Computed observables: (a),(d)
local configurations of fermionic matter for each lattice site $\left\langle
N_{\mathbf{j}}\right\rangle$, local spin along the $z$-axis
$S^{z}_{\mathbf{j}}$ represented by arrows in the center of the lattice sites
(values are close to zero), rishon modes occupation on the links
$n_{\text{rishon},\mathbf{j}}$; (b),(e) spin modulus squared
$S^{2}_{\mathbf{j}}$; (c),(f) variance of the matter occupation number
$\mathrm{Var}(N_{\mathbf{j}})$.

Given the defermionized Hubbard Hamiltonian of Eqs. (22)-(23), the ground
state at a specified bond dimension $\chi$ is determined by iteratively
optimizing each of the tensors in TTN, gradually reducing the energy
expectation value. This procedure is iterated several times to reach the
desired convergence. For a detailed and exhaustive description of the
algorithms, please see the technical reviews and textbooks Silvi  _et al._
(2019); Montangero (2018). In our TTN simulations, we use a maximum bond
dimension $\chi=350$, which is sufficient to reach a relative error on the
energy density of the order $10^{-6}$, ensuring the stability of our findings.

In Fig. 3(a), we show the TTN results concerning the defermionized Hamiltonian
on a $4\times 4$ lattice at half-filling, i.e. $\rho=1$. We also report the
corresponding energy densities obtained via ED for the original Hubbard
Hamiltonian at smaller lattice sizes, and the results of the extrapolation of
the infinite-size limit obtained by using iPEPS methods Corboz (2016). The
energies obtained with the TTN simulations for the defermionized Hamiltonian
are in agreement with the overall scaling shown by the exact energies as a
function of $U/t$ and the lattice sizes.

To further test the equivalence between the original Hubbard model and its
defermionized version, we compute the ground-state $D\equiv
1/|\Lambda|\sum_{\mathbf{j}}\left<n_{\mathbf{j},\uparrow}n_{\mathbf{j},\downarrow}\right>$.
In Fig. 3(b), we show the results obtained from the TTN simulations of the
defermionized model on the $4\times 4$ lattice. For the sake of comparison, we
report the data of the original Hubbard Hamiltonian, obtained in Ref. LeBlanc
_et al._ (2015) by using two independent numerical methods, i.e. the
multireference projected Hartree-Fock method (MRPHF) and the Density matrix
embedding theory (DMET). A good agreement and consistency of our results with
the reference data is visible for all the simulated values of $U/t$.

TTN simulations allow for an extra characterization of the ground state in
terms of local observables, such as the occupation of fermionic matter and its
variance

$\displaystyle\left\langle N_{\mathbf{j}}\right\rangle$
$\displaystyle=\left\langle
n_{\mathbf{j},\uparrow}+n_{\mathbf{j},\downarrow}\right\rangle$
$\displaystyle\mathrm{Var}(N_{\mathbf{j}})$
$\displaystyle=\langle\quantity(N_{\mathbf{j}})^{2}\rangle-\left\langle
N_{\mathbf{j}}\right\rangle^{2},$ (24)

the local spin along the $z$-axis and its square modulus,

$\displaystyle S^{z}_{\mathbf{j}}$ $\displaystyle=\frac{1}{2}\left\langle
n_{\mathbf{j},\uparrow}-n_{\mathbf{j},\downarrow}\right\rangle$ $\displaystyle
S^{2}_{\mathbf{j}}=\frac{3}{4}\left\langle(n_{\mathbf{j},\uparrow}-n_{\mathbf{j},\downarrow})^{2}\right\rangle,$
(25)

but also the rishon mode occupation on lattice links
$n_{\text{rishon},\mathbf{j}}=\langle
c^{\dagger}_{\mathbf{j},{\bf\it\mu}}c_{\mathbf{j},{\bf\it\mu}}\rangle$.

In Fig. 4, we show the results for $U/t=0$ and $U/t=4$. In both cases, as
shown in Fig. 4(a) and Fig. 4(d), the fermionic occupation number
$N_{\mathbf{j}}\sim 1$, while the local spin along $z$, represented by arrows
in the center of the lattice sites, is close to zero; correspondingly, all the
rishon modes occupations on the links are $n_{\text{rishon},\mathbf{j}}\sim
0.5$, satisfying the penalty constraints that encode the fermion parity. By
looking at Figs. 4(b) and 4(e), we observe that the squared spin modulus
$S^{2}_{\mathbf{j}}$ increases for all lattice sites by varying $U/t$ from
$0.0$ to $4.0$, whereas the variance of the fermionic occupation number
$\mathrm{Var}(N_{\mathbf{j}})$ decreases, as shown in Figs. 4(c) and 4(f).

These configurations agree with the expected ground state of the original
Hubbard model. Indeed, in the case of half-filling, the Hubbard model maps to
the Heisenberg model Auerbach (1994), and the ground state mimics an
antiferromagnetically long-range ordered state.

##  III Digital quantum simulation

The gauge ‘defermionization’ offers a natural encoding of fermionic degrees of
freedom to qubits, thus allowing for the digital quantum simulation of
fermionic models on quantum computers.

Here, we break down the procedure to express the spin-$\frac{1}{2}$ Hubbard
defermionized model on a square lattice in terms of qubits and Pauli operators
on a generic platform, see Sec. III.1, Sec. III.3, and Sec. III.2. This gauge-
field-based encoding is genuinely local, i.e. the lattice support of each
Hamiltonian term is not increased in the encoding, and thus each Pauli weight
$W$ is conserved and does not depend on the system size $N$. However, locality
comes at the cost of including auxiliary qubits representing the defermionized
modes, i.e. fermionic modes and $\mathbb{Z_{2}}$ gauge fields. There currently
exists a few local mappings Bravyi and Kitaev (2002); Verstraete and Cirac
(2005); Whitfield  _et al._ (2016); Chen  _et al._ (2018); Steudtner and
Wehner (2019); Setia  _et al._ (2019); Jiang  _et al._ (2019); Derby  _et al._
(2021) and the relation between them is analyzed in Chen  _et al._ (2023).
Fermion-to-qubit mappings need resources that are usually quantified in terms
of the number of qubits to simulate one fermion on average (qubit-fermion
ratio), the number of qubits to express the parity of a fermionic state
(parity weight), the length of the hopping operators expressed as Pauli
strings and the maximum length of Pauli strings for the stabilizers. The
resources required by the gauge-field-based encoding are comparable with
state-of-the-art methods: values for qubit-fermion ratio, fermion parity
weight, hopping weight, and stabilizer weight are $(3,1,6,6)$ respectively.
For an exhaustive comparison across the available encodings, see Table I in
Chen  _et al._ (2023).

Following these prescriptions, we look for nontrivial dynamical effects. The
spin-charge separation is a manifestation of the strongly correlated physics
in one dimension (1D), where spin and charge excitations become decoupled and
propagate at distinct velocities Giamarchi (2003). While the theory is well-
established Kollath  _et al._ (2005) and it has been recently observed both in
cold Fermi gases Vijayan  _et al._ (2020); Bohrdt  _et al._ (2021) and digital
quantum simulations Arute  _et al._ (2020), the persistence of spin-charge
separation in higher dimensions is still a subject of ongoing debate. In
recent decades, it has become evident that conclusive answers can only be
attained through numerical approaches Yang  _et al._ (2022); Chen  _et al._
(1994), with only a handful of experiments attempting to address this issue
Wang  _et al._ (2022). For these reasons, we simulate the digital dynamics of
spin and charge excitations over a half-filled ($\rho=1$) 2 $\times$ 4 system
in the antiferromagnetic phase $U/t=10$. In Sec. III.4, we lay out the
protocols to adiabatically prepare the antiferromagnetic ground state, and
then inject the excitations. Finally, in Sec. III.5, we observe signatures of
the spin-charge separation phenomenon. While protocols are designed for a
general digital quantum computer, the presented results are obtained from
tensor networks qma  to emulate the evolution of the quantum circuit.

###  III.1 Fermion to qubits mapping of the defermionized Hubbard Hamiltonian

We introduce a qubit for each flavor (up $u$ and down $d$) and one for each
rishon (north $n$, west $w$, east $e$, south $s$). To minimize the Pauli
weight of the Hamiltonian terms, a specific ordering for the qubits composing
each dressed site is defined:

$\displaystyle\begin{cases}\\{u,d,w,s,e,n\\}&\text{if $(-1)^{j_{x}+j_{y}}=1$
(even site)}\\\ \\{d,u,s,w,n,e\\}&\text{if $(-1)^{j_{x}+j_{y}}=-1$ (odd
site)}.\end{cases}$ (26)

All the operators in Eq. 22 and penalties in Eq. 23 are mapped to a
spin-$\frac{1}{2}$ algebra, while preserving all the commutation relations.
Following the order defined in Eq. 26, the operators acting on an even site
$j$ can be written as:

$\displaystyle\psi^{\dagger}_{\mathbf{j},\uparrow}=XZZZZZ,\quad$
$\displaystyle\psi^{\dagger}_{\mathbf{j},\downarrow}=IXZZZZ,$ (27)
$\displaystyle\psi_{\mathbf{j},\uparrow}=YZZZZZ,\quad$
$\displaystyle\psi_{\mathbf{j},\downarrow}=IYZZZZ,$ (28)
$\displaystyle\gamma_{\mathbf{j},r}=I_{u}I_{d}\bigotimes_{k=w}^{r-1}I_{k}$
$\displaystyle\otimes X_{r}\bigotimes_{k=r+1}^{n}Z_{k},$ (29)

where $X,Y,Z$ are Pauli matrices and $I$ is the identity. In the two equations
above we dropped the index of species and rishons for clarity. The string of
$Z$ operators applies to all the qubits forming a site. The index $r$ runs
over the different rishon species of a site, as defined in Eq. 26. The mapping
for an odd site is obtained by swapping the first two operators and using the
odd-ordered sites.

Now, we compute the explicit form of the Hamiltonian terms, hopping and on-
site interaction, in terms of Pauli operators. Following the notation defined
in Eq. 26, the horizontal hopping operator from an odd ($O$) to an even site
($E$) for flavor $up$ reads:

| u | d | w | s | e | n | d | u | s | w | n | e
---|---|---|---|---|---|---|---|---|---|---|---|---
$\psi_{\uparrow,E}^{\dagger}$ | X | Z | Z | Z | Z | Z | I | I | I | I | I | I
$\gamma_{eE}$ | I | I | I | I | X | Z | I | I | I | I | I | I
$\gamma_{wO}$ | I | I | I | I | I | I | I | I | I | X | Z | Z
$\psi_{\uparrow,O}$ | I | I | I | I | I | I | I | Y | Z | Z | Z | Z
hopping | X | Z | Z | Z | iY | I | I | Y | Z | -iY | I | I
(30)

Matrix multiplication goes from bottom to top. By computing also the hermitian
conjugate, the hopping term results:

$\displaystyle\psi_{\uparrow E}^{\dagger}\gamma_{eE}$
$\displaystyle\gamma_{wO}\psi_{\uparrow O}-\psi_{\uparrow
E}\gamma_{eE}\gamma_{wO}\psi_{\uparrow O}^{\dagger}=$ (31)
$\displaystyle\left(XZZZYIIYZYII-YZZZYIIXZYII\right).$

Terms for vertical hopping or hoppings involving different species can be
derived analogously.

Similarly, we derive the on-site interaction term. The number operator acting
on site $\mathbf{j}$ can be written as
$n_{\mathbf{j},\sigma}=\frac{1}{2}\left(1-Z_{\mathbf{j},\sigma}\right)$, then
the on-site interaction reads:

$\displaystyle\left(n_{\mathbf{j},\uparrow}-\frac{1}{2}\right)\left(n_{\mathbf{j},\downarrow}-\frac{1}{2}\right)=\frac{1}{4}Z_{\mathbf{j},\uparrow}Z_{\mathbf{j},\downarrow}.$
(32)

Once we obtain the Hamiltonian terms described as pauli strings, we decompose
the relative time propagator in single and two-qubit gates following Ref.
Whitfield  _et al._ (2016). A detailed description of the mapping from the
time propagator to the quantum circuit can be found in Appendix A.

###  III.2 Vertex and plaquette constraints

Following the prescription of the gauge defermionization, we have to constrain
the state to a specific subspace of the Hilbert space, a condition set by the
first two lines of Eq. 15. The gauge symmetries can be cast into a language
that directly translates into a quantum circuit scenario. Each constraint is
equivalent to a stabilizer of a quantum error correcting code; to respect the
gauge symmetry, the state is restrained to the $+1$ eigenstate of the
stabilizers Gottesman (1997). The Pauli representation of the stabilizers can
be obtained following the mapping presented in Eq. 27. As expected, all the
stabilizers commute with the Hamiltonian.

The vertex stabilizer term reads from the first line of Eq. 15. As the
exponentials are counting the parity of a given species/rishon, these
operators correspond to $Z$ operators in the qubit language. Thus, the
stabilizer acting on the lattice vertex $\mathbf{v}$ results in

$\displaystyle
S_{\mathbf{v}}=\bigotimes_{f=u,d}Z_{f}\bigotimes_{\mu\in\mathbf{v}}Z_{\mathbf{v}+{\bf\it\mu}},$
(33)

where $\mu$ runs on all rishons connected to the vertex $\mathbf{v}$. This
stabilizer must be applied to each vertex $\mathbf{v}$ in the system.

The plaquette stabilizer reads from the second row in Eq. 15. As the product
always runs on four terms, the imaginary unit $i$ drops. Following the
representation of the $\gamma$ operators reported in Eq. 29, the final form of
each plaquette term respects the qubits ordering of the dressed site. Thus,
the form of an even plaquette stabilizer is different from an odd one, where a
plaquette is said even (odd) if the lower-left corner of the plaquette lies on
an even (odd) site. We report the explicit computation for an even plaquette
in Table 1.

Finally, we report the size of the Hilbert space to perform a digital quantum
simulation of the defermionized Hubbard model over a $x\times y$ grid, and
compare it to the size of the accessible subspace:

$\displaystyle\text{dim}(\mathcal{H})$
$\displaystyle=2^{2xy+x(y-1)+y(x-1)}=2^{4xy-x-y}$ (34)
$\displaystyle\text{dim}(\mathcal{H}_{eff})$
$\displaystyle=2^{2xy+x(y-1)+y(x-1)-(x-1)(y-1)-xy}$ $\displaystyle=2^{2xy-1}.$
(35)

We notice that constraining the system to the physical subspace quadratically
decreases the dimension of the Hilbert space, showing that we are working with
a highly redundant space.

| u | d | w | s | e | n | d | u | s | w | n | e | d | u | s | w | n | e | u | d | w | s | e | n
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\gamma_{00e}$ | I | I | I | I | X | Z | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I
$\gamma_{00n}$ | I | I | I | I | I | X | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I
$\gamma_{01s}$ | I | I | I | I | I | I | I | I | X | Z | Z | Z | I | I | I | I | I | I | I | I | I | I | I | I
$\gamma_{01e}$ | I | I | I | I | I | I | I | I | I | I | I | X | I | I | I | I | I | I | I | I | I | I | I | I
$\gamma_{10w}$ | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | X | Z | Z | I | I | I | I | I | I
$\gamma_{10n}$ | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | X | Z | I | I | I | I | I | I
$\gamma_{11w}$ | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | X | Z | Z | Z
$\gamma_{11s}$ | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | I | X | Z | Z
stabilizer | I | I | I | I | X | iY | I | I | X | Z | Z | iY | I | I | I | X | iY | I | I | I | X | iY | I | I
Table 1: Explicit definition of the stabilizer for an even plaquette. This
calculation does not take into account the link symmetry: thus it is necessary
to do that further computation to obtain the final expression of the qubit
mapping.

###  III.3 Link constraint

The state of the two rishons sharing the same link is also subject to a
constraint set by the third line of Eq. 15. The only qubits states fulfilling
this constraint are $\ket{00}$ and $\ket{11}$. The effective Hilbert space for
a pair of rishons is $2$-dimensional; for this reason, both rishons can be
described with a single qubit, by mapping
$\ket{00}\rightarrow\ket{0},\ket{11}\rightarrow\ket{1}$. After this procedure,
we need to project all the operators acting on the two rishons on the new
subspace. Given an operator $A_{r_{1}r_{2}}$ we can write the new operator
living on the two-dimensional space as:

$\displaystyle A_{\overline{r}}=\text{Tr}_{\ket{01},\ket{10}}A_{r_{1}r_{2}},$
(36)

where we trace out the states violating the constraint. We report the new form
of some operators of interest as an example: $ZZ\rightarrow I$, $XX\rightarrow
X$. This way, each dressed site shares a rishon with its nearest neighbors,
and thus there is no longer a clear separation between the dressed sites. In
Fig. 5, we report the final form of the stabilizers and the Hamiltonian terms,
once the link constraint is applied.

Figure 5: Graphical representation of the Hamiltonian terms in the
defermionized Hubbard model and of the stabilizers for installing the gauge
symmetries. These terms result from applying Eq. 36, where the states of the
two rishons satisfying the link constraint are mapped to a single qubit for
each link. In (a), we depict the vertex stabilizer, while (b) corresponds to
the even plaquette stabilizer, and (c) represents the odd plaquette
stabilizer. Figures (d) and (e) illustrate the hopping term and the onsite
interaction respectively. As dictated by Eq. 26, in figure (d), the $u$ qubit
is the first on the even site and the second on the odd site.

###  III.4 Preparation of spin and charge excitations

Here, we lay out the protocols to adiabatically prepare the antiferromagnetic
ground state, and to inject spin- and charge-excitations.

#### Adiabatic ground state preparation

The adiabatic state preparation targets the ground state of the defermionized
Hubbard model over a $4\times 2$ lattice, at half-filling ($\rho=1$) with
$t=0.1$ and $U=1$. The adiabatic evolution starts from the ground state of Eq.
22 with $t=0$, here renamed $H_{0}=H^{\prime\prime}_{Hub}(t=0)$; the hopping
terms are slowly switched on to reach $H_{1}=H^{\prime\prime}_{Hub}$. Then,
the time-dependent Hamiltonian reads:

$\displaystyle H$ $\displaystyle=(1-\beta)H_{0}+\beta H_{1},$ (37)

where $\beta\in[0,1]$ is the adiabatic parameter that increases linearly for
$100$ steps; for each $\beta$, the system evolves for ten time-steps $d\tau$
to ensure a smooth convergence, for a total of $1000$ evolution steps. The
real-time evolution is simulated by decomposing the evolution operator via a
first-order Trotterization with time step $dt=d\tau=0.01$.

The initial state for the adiabatic process is the ground state of $H_{0}$ at
half-filling ($\rho=1$) and respects the spin-flip symmetry. While neglecting
the rishons’ state, we consider only the state of the matter qubits ($u$ and
$d$) for each dressed site over the $4\times 2$ lattice, where
$\ket{10}=\ket{\uparrow}$ and $\ket{01}=\ket{\downarrow}$. It results in

$\displaystyle\frac{1}{\sqrt{2}}\left(\left|\begin{matrix}\uparrow&\downarrow&\uparrow&\downarrow\\\
\downarrow&\uparrow&\downarrow&\uparrow\end{matrix}\right\rangle+\left|\begin{matrix}\downarrow&\uparrow&\downarrow&\uparrow\\\
\uparrow&\downarrow&\uparrow&\downarrow\end{matrix}\right\rangle\right).$ (38)

Notice that the state in Eq. 38 is a GHZ state, thus it can be easily prepared
using only Hadamard, NOT, and controlled-NOT gates. Finally, the rishons’
state is uniquely determined once the stabilizers are applied, see Fig. 5.

At the end of the adiabatic evolution, we checked that local spin projection
$S^{z}_{\mathbf{j}}$ and local charge $N_{\mathbf{j}}$ are stationary up to
oscillations compatible with the cut of the singular values in the tensor
network emulator ($10^{-8}$). For a discussion about convergence see Appendix
B.

Figure 6: Schematic of injecting spin- and charge-excitations on the $(0,0)$
corner. Both figures depict the qubit lattice, highlighting up and down qubits
within the matter site. Rishons are represented as black dots. All applied
operations commute with the stabilizers. (a) Spin excitation: To inject the
spin excitation, we initially measure the qubit states for up and down, post-
selecting only the $\ket{01}=\ket{\uparrow}$ state. Then, the $XX$ operator
flips the spin state. (b) Charge excitation: Injecting the charge excitation
involves measuring the up qubit and post-selecting the state $\ket{1}$. We
then apply the Pauli operators as depicted in the figure, equivalent to
causing the matter qubit $u$ to jump outside the lattice. To achieve this, an
additional rishon (the $s$ rishon in this case) must be added.

#### Excitations

Spin- and charge-excitations are injected in the ground state of $H_{1}$ using
local operations and classical communication (LOCC)s: this procedure is, in
principle, completely reproducible on state-of-the-art quantum computers. Both
protocols preserve the symmetries of the system, namely every constraint
previously defined is satisfied. In Fig. 6, we provide a graphical
representation of the two methods.

Spin excitation. Without losing any generality, the spin excitation is created
on the site $(0,0)$. Recalling that qubits $u,d$ represent up and down
flavors, we choose a general superposition of the two states; after a
projective measurement on both qubits, only the flavor up is post-selected,
i.e. the qubits state $\ket{10}=\ket{\uparrow}$. Then, a spin excitation is
introduced by flipping the spin state:

$\displaystyle\ket{\uparrow}$ $\displaystyle\longrightarrow\ket{\downarrow},$
$\displaystyle\ket{1_{u}0_{d}}$
$\displaystyle\xrightarrow{XX}\ket{0_{u}1_{d}}.$

This operation corresponds to applying the operator $XX$ on the qubits $u,d$.
Notice that $XX$ commutes with all the stabilizers, since it shares support
only with the vertex stabilizer in $(0,0)$, and they commute.

Charge excitation. In principle, the charge excitation can be created anywhere
on the lattice. However, for improved efficiency, we restrict the protocol to
the sites along the border and in particular to the site $(0,0)$. A projective
measurement is performed on the $u$ qubit, then post-selecting the state
$\ket{1}$. Being left with a single charge on the site, the excitation is
obtained by removing that very charge. However, this operation does not
commute with the vertex stabilizer in $(0,0)$, thus breaking that constraint.
To preserve Gauss’ law we introduce an additional qubit, labeled as an extra
rishon of site $(0,0)$, in this case, $s$ (or $w$). The excitation is
implemented by flipping the state of the qubits $u,s$:

$\displaystyle\ket{1_{u}0_{s}}$ $\displaystyle\stackrel{{\scriptstyle
XX}}{{\longrightarrow}}\ket{0_{u}1_{s}}.$

This procedure, however, does not commute with the stabilizers, and it would
bring the state outside the subspace of the physical states. To avoid this
issue, we implement this operation by applying the hopping operator
$i\gamma_{\mathbf{j},s}\psi_{\mathbf{j},\uparrow}$ with $\mathbf{j}=(0,0)$,
which makes the particle in $u$ hop through $s$ outside the lattice. The
operator reads $g_{charge}=X_{u}Z_{d}Z_{w}Y_{s}$, and it is equivalent to the
bit-flip just discussed. We stress that qubit $d$ is not projected during this
protocol.

###  III.5 Propagation of spin and charge excitations

By employing the local encoding, we simulate the dynamics of spin and charge
excitations over a half-filled ($\rho=1$) $2\times 4$ system in the
antiferromagnetic phase $U/t=10$. To this aim, we use quantum matcha TEA qma ,
an emulator of quantum circuits based on MPS. All the simulations converge
with bond dimension $\chi=1024$ \- below the maximum achievable for $27$
qubits, where the maximum bond dimension is bound by
$\chi_{\mathrm{max}}=8192$ \- and Trotter step $dt=0.01$. We find that the two
kinds of excitations propagate at different speeds, in agreement with the
spin-charge separation phenomenon.

Figure 7: Parity measurement to compute the correlator $\langle
Z_{\mathbf{j},u}Z_{\mathbf{j},d}\rangle(\tau)$ using an ancilla qubit,
depicted in blue. Only controlled-NOT gates and measurements are used. Figure
8: Digital quantum simulation of spin-charge separation: (a) Propagation of
the charge profile on each site at different times in the evolution
$\tau/U=1,2,3$. (b) Charge profile for selected sites on the lattice. (c)
Propagation of the spin profile on each site at different times in the
evolution $\tau/U=1,2,3$. (d) Spin profile for selected sites on the lattice.

The first column of Fig. 8 shows the evolution of the injected charge
excitation. By measuring the local charge $N_{\mathbf{j}}(\tau)$, one can
monitor in time the deviations from the initial condition, i.e. the excitation
injection:

$\displaystyle N_{\mathbf{j}}(\tau)=\langle
n_{\mathbf{j},\uparrow}\rangle(\tau)+\langle
n_{\mathbf{j},\downarrow}\rangle(\tau),$ (39)

where $\tau$ is the time dependence. This assessment requires only the
expectation values of local observables, since

$\displaystyle\langle
n_{\mathbf{j},\genfrac{}{}{0.0pt}{}{\uparrow}{\downarrow}}\rangle$
$\displaystyle=\frac{1-\langle
Z_{\mathbf{j},\genfrac{}{}{0.0pt}{}{u}{d}}\rangle}{2},$ where
$\displaystyle\langle Z_{\mathbf{j},\genfrac{}{}{0.0pt}{}{u}{d}}\rangle$ (40)

is the expectation value of the $u(d)$ qubits over the $Z$ basis. The second
column of Fig. 8 shows the evolution of the spin excitation. Similarly, one
can assess the deviation from the initial condition by measuring the local
spin modulus squared $S_{\mathbf{j}}^{2}(\tau)$:

$\displaystyle S_{\mathbf{j}}^{2}(\tau)=\langle
n_{\mathbf{j},\uparrow}\rangle(\tau)+\langle
n_{\mathbf{j},\downarrow}\rangle(\tau)-2\langle
n_{\mathbf{j},\uparrow}n_{\mathbf{j},\downarrow}\rangle(\tau).$ (41)

In this scenario, we also need to measure a two-point correlation function,
i.e. $\langle n_{\mathbf{j},\uparrow}n_{\mathbf{j},\downarrow}\rangle(\tau)$.
This quantity can be rewritten again in terms of expectation values over the
$u,d$ qubits using Eq. 40. The resulting correlator, $\langle
Z_{\mathbf{j},u}Z_{\mathbf{j},d}\rangle(\tau)$, is a parity measurement that
can be computed with the circuit reported in Fig. 7.

For three subsequent time steps, the first row of Fig. 8 reports a schematic
of the above-measured quantities (local charge and spin squared modulus) over
the lattice. With this plot, one can easily visualize the difference in the
propagation speed between the charge-excitation and the spin-excitation, see
Fig. 8(a) and Fig. 8(c), respectively. When the spin excitation reaches the
border, at $\tau/U=3$, the charge excitation has yet to explore most of the
lattice. To analyze this behavior more quantitatively we report, on the second
row, the time evolution up to $\tau/U=50$ for both quantities. The charge
excitation in Fig. 8(b) is oscillating much slower than the spin excitation in
Fig. 8(d). These results, i.e. a different propagation of the perturbation
when introducing spin- and charge-excitations are consistent with the
dynamical separation of the two types of degrees of freedom.

##  IV Conclusions

We have generalized a technique for local fermion encoding to any 2D lattice
configurations that eliminates the fermionic degrees of freedom by absorbing
them into an auxiliary $\mathbb{Z}_{2}$ gauge field (gauge defermionization).
We have successfully tested this method against the 2D spin-$\frac{1}{2}$
Hubbard model. The ground state properties have been computed for varying
particle densities utilizing both tree-tensor network ansatz and exact
diagonalization methods. Here, we have observed the expected transition from
the liquid to the anti-ferromagnetically ordered phase, accessing lattice
sizes up to $4\times 4$. Furthermore, gauge defermionization introduces a
fermion-to-qubit mapping with comparable resource requirements to state-of-
the-art encodings. We have shown that this mapping offers a scalable pathway
for the digital quantum simulation of fermionic theories. For the
defermionized 2D Hubbard Hamiltonian, hopping terms and gauge constraints,
here included as stabilizers, result in a maximum Pauli weight of 6. We have
then demonstrated the feasibility of this approach by simulating the digital
out-of-equilibrium dynamics of the defermionized 2D Hubbard Hamiltonian over a
$4\times 2$ lattice in the antiferromagnetic phase. Our protocol entails the
adiabatic preparation of the ground state (half-filling, $U/t=10$), the
injection of a charge (spin)-excitation in the system, and the time evolution
of the perturbed state. Finally, we have observed a slower propagation for the
charge excitation concerning the spin one, namely the spin-charge separation
effect.

Our collected results show numerical evidence that the local encoding is a
scalable and feasible pathway toward the digital quantum simulation of
fermionic lattice theories. We have highlighted that is indeed feasible to
provide physically relevant results with the current technology of digital
quantum processing platforms.

Tensor networks also showed positive results, although specifically the tree
tensor network ansatz state manifests some limitations to accommodate the
area-law of entanglement introduced by the auxiliary, resonant gauge fields.
We expect tensor network geometries capable of capturing a wider entanglement
distribution, such as ATTN Felser  _et al._ (2021) or iPEPS Jordan  _et al._
(2008), to yield even better results.

##  V Acknowledgments

We thank M. Di Liberto and E. Ercolessi for their helpful discussions. This
work is partially supported by the Italian PRIN2017 and PRIN2022 project
TANQU, the INFN project QUANTUM, the BMBF via the project QRydDemo, the EU via
QuantERA projects QuantHEP and T-NiSQ, the Horizon 2020 research and
innovation program (Quantum Flagship - PASQuanS2, EuRyQa, textarossa), the
European Union - NextGenerationEU project CN00000013 - Italian Research Center
on HPC, Big Data and Quantum Computing. We acknowledge computational resources
by the Cloud Veneto, CINECA, the BwUniCluster, the University of Padova
Strategic Research Infrastructure Grant 2017: “CAPRI: Calcolo ad
Alprestazionioni per la Ricerca e l’Innovazione”, and the WCRI-Quantum
Computing and Simulation Center of Padova University.

##  Appendix A Mapping the time propagator to a quantum circuit

To simulate the dynamics of our system on a quantum computer it is important
to map the evolution operator to a quantum circuit using single and two-qubit
gates since those are the available operations. First, we Trotterize the
evolution treating each Hamiltonian term separately. Then, it has been shown
that the propagator $e^{-\frac{i}{\hbar}JH_{i}dt}$ of a $H_{i}=ZZ\dots Z$
interaction with strength $J$ can be compiled combining a cascade of CNOTs
gates and a rotation along the $z$ axis Whitfield  _et al._ (2016):

$\displaystyle R_{z}^{\theta}=\begin{pmatrix}1&0\\\
0&e^{i\theta}\end{pmatrix},\quad\theta=Jdt.$ (42)

In general, we can treat an arbitrary Pauli string by moving to the $Z$ basis
along the CNOTs cascade. Practically, we apply the following single-qubit
basis change gate before and after the CNOT, respectively for the $X$ and $Y$
Pauli matrices:

$\displaystyle H=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix},\quad
R_{x}\left(-\frac{\pi}{2}\right)=\frac{1}{\sqrt{2}}\begin{pmatrix}1&i\\\
i&1\end{pmatrix}.$ (43)

In Figure 9 we report the form of the CNOTs cascade with the basis change,
showing an example for a general Pauli string.

Figure 9: Mapping of the time propagator of a generic Pauli string to a
quantum circuit using only single and two-qubit gates.

##  Appendix B Convergence checks

We perform different checks to ensure that the results obtained from the time
evolution of the MPS are meaningful. First, one should ensure that the
trotterization step is small enough. We perform a time evolution up to time
$t/U=20$ and then invert the time, letting the system evolve back. If we are
not making any trotter errors, we should go back to the very same state. We
estimate the error as the absolute value of the difference of the observable
$\langle n_{\uparrow,j}n_{\downarrow,j}\rangle$ at $\tau=0$ and at the end of
the evolution. In Fig. 10, we show the scaling of the error for different
trotter steps. We can safely state that the chosen trotter step of
$dt=10^{-2}$ is small enough to obtain reasonable results.

Figure 10: Trotter error as a function of the trotter step in a log-log scale.

The second source of error is the state truncation due to the finite bond
dimension of the system, which is set at $\chi=2^{10}=1024$. This bond
dimension is lower than the maximum bond dimension that a system of $27$ would
require, i.e. $\chi_{\mathrm{max}}=2^{13}=8192$. However, we can monitor the
approximation we make Jaschke and Montangero (2022): after each two-qubits
gate, we perform at most an error $E=\sum_{i}s_{i}^{2}$, where $s_{i}$ are the
singular values truncated by the MPS algorithm. A lower bound of the fidelity
$\mathcal{F}$ of the final state of the simulation is:

$\displaystyle\mathcal{F}=\prod_{k}(1-E_{k}),$ (44)

where the index $k$ runs over all the gates of the quantum circuit
implementing the time evolution. We obtained infidelity of
$1-\mathcal{F}_{s}=4\cdot 10^{-7}$ for the evolution after the application of
the spin excitation, and $1-\mathcal{F}_{c}=2\cdot 10^{-9}$ for the evolution
after the insertion of the charge excitation. We can thus state that all
sources of numerical errors are under control.

##  Appendix C Initial state for the adiabatic evolution

In Sec. III.4, we discussed the adiabatic preparation of the ground state of
the defermionized Hubbard Hamiltonian for a $4\times 2$ lattice at half
filling. However, we can still choose the initial state for such a procedure,
as long as it is a ground state of $H_{0}$. It is thus convenient to choose an
easy-to-prepare state that maximizes the overlap with the expected final
state. The following discussion will ignore the states of the rishons, which
are always considered to be in the correct state to ensure that all the
constraints are satisfied. Since we are interested in the strongly repulsive
regime of $U/t=10$, $U>0$, we can expect to find only one particle in a given
site, i.e., the up $u$ or down $d$ spin species. This is equivalent to saying
that the only qubit states we allow are $\ket{1_{u}0_{d}},\ket{0_{u}1_{d}}$.
We also notice that any combination of states that satisfy this condition is a
ground state of the defermionized Hubbard Hamiltonian in Eq. 22 with $t=0$.
First, we consider that, once the hopping is turned on, some configurations
will be favored: those that have a checkerboard configuration of sites, where
we alternate the up and down species, i.e., if the site $\mathbf{j}$ is in the
state $\ket{1_{u}0_{d}}$ all its neighboring sites will be in
$\ket{0_{u}1_{d}}$. We finally prepare a superposition of the checkerboard
states, such that the resulting state is symmetric under a global spin flip,
i.e., a flipping of the $u,d$ qubits since we expect this to be a symmetry of
the Hamiltonian. The protocol to prepare this state on a quantum computer is
as follows:

  1. 1.

Prepare the checkerboard state;

  2. 2.

Stabilize the state under the global spin flip stabilizer, i.e. an $X$
operator acting on all the $u,d$ qubits;

  3. 3.

Stabilize the state w.r.t. all the local stabilizer defined in Sec. III.

This state is an optimal choice of the initial state of the simulation.

## References

  * Fradkin (2013) E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Auerbach (1994) A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer New York, NY, 1994).
  * Giamarchi (2003) T. Giamarchi,  _Quantum physics in one dimension_ , Vol. 121 (Clarendon press, 2003).
  * Lee  _et al._ (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Hartke  _et al._ (2023) T. Hartke, B. Oreg, C. Turnbaugh, N. Jia, and M. Zwierlein, Science 381, 82 (2023).
  * Halperin and Jain (2020) B. I. Halperin and J. K. Jain,  _Fractional Quantum Hall Effects: New Developments_ (World Scientific, 2020).
  * Hémery  _et al._ (2023) K. Hémery, K. Ghanem, E. Crane, S. L. Campbell, J. M. Dreiling, C. Figgatt, C. Foltz, J. P. Gaebler, J. Johansen, M. Mills, S. A. Moses, J. M. Pino, A. Ransford, M. Rowe, P. Siegfried, R. P. Stutz, H. Dreyer, A. Schuckert, and R. Nigmatullin, “Measuring the Loschmidt amplitude for finite-energy properties of the Fermi-Hubbard model on an ion-trap quantum computer,”  (2023), arxiv:2309.10552 [quant-ph] .
  * Hofrichter  _et al._ (2016) C. Hofrichter, L. Riegger, F. Scazza, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Physical Review X 6, 021030 (2016).
  * Gattringer and Langfeld (2016) C. Gattringer and K. Langfeld, International Journal of Modern Physics A 31, 1643007 (2016), https://doi.org/10.1142/S0217751X16430077 .
  * Troyer and Wiese (2005) M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).
  * Altman  _et al._ (2021) E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, M. Greiner, K. R. Hazzard, R. G. Hulet, A. J. Kollár, B. L. Lev, M. D. Lukin, R. Ma, X. Mi, S. Misra, C. Monroe, K. Murch, Z. Nazario, K.-K. Ni, A. C. Potter, P. Roushan, M. Saffman, M. Schleier-Smith, I. Siddiqi, R. Simmonds, M. Singh, I. Spielman, K. Temme, D. S. Weiss, J. Vučković, V. Vuletić, J. Ye, and M. Zwierlein, PRX Quantum 2, 017003 (2021).
  * Esslinger (2010) T. Esslinger, Annual Review of Condensed Matter Physics 1, 129 (2010), https://doi.org/10.1146/annurev-conmatphys-070909-104059 .
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbène, Nature Physics 8, 267 (2012).
  * Tarruell and Sanchez-Palencia (2018) L. Tarruell and L. Sanchez-Palencia, Comptes Rendus Physique 19, 365 (2018).
  * Bohrdt  _et al._ (2021) A. Bohrdt, L. Homeier, C. Reinmoser, E. Demler, and F. Grusdt, Annals of Physics 435, 168651 (2021), special issue on Philip W. Anderson.
  * Cheuk  _et al._ (2015) L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Physical Review Letters 114, 193001 (2015).
  * Duarte  _et al._ (2015) P. M. Duarte, R. A. Hart, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, and R. G. Hulet, Physical Review Letters 114, 070403 (2015).
  * Edge  _et al._ (2015) G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, and J. H. Thywissen, Physical Review A 92, 063406 (2015).
  * Greif  _et al._ (2013) D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, and T. Esslinger, Science 340, 1307 (2013).
  * Haller  _et al._ (2015) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Nature Physics 11, 738 (2015).
  * Hart  _et al._ (2015) R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse, and R. G. Hulet, Nature 519, 211 (2015).
  * Hofstetter  _et al._ (2002) W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. D. Lukin, Physical Review Letters 89, 220407 (2002).
  * Jördens  _et al._ (2008) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature 455, 204 (2008).
  * Messer  _et al._ (2015) M. Messer, R. Desbuquois, T. Uehlinger, G. Jotzu, S. Huber, D. Greif, and T. Esslinger, Physical Review Letters 115, 115303 (2015).
  * Murmann  _et al._ (2015) S. Murmann, F. Deuretzbacher, G. Zürn, J. Bjerlin, S. M. Reimann, L. Santos, T. Lompe, and S. Jochim, Physical Review Letters 115, 215301 (2015).
  * Omran  _et al._ (2015) A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Physical Review Letters 115, 263001 (2015).
  * Parsons  _et al._ (2015) M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Physical Review Letters 114, 213002 (2015).
  * Schneider  _et al._ (2008) U. Schneider, L. Hackermüller, S. Will, Th. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008).
  * Taie  _et al._ (2012) S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi, Nature Physics 8, 825 (2012).
  * Uehlinger  _et al._ (2013) T. Uehlinger, G. Jotzu, M. Messer, D. Greif, W. Hofstetter, U. Bissbort, and T. Esslinger, Physical Review Letters 111, 185307 (2013).
  * Zache  _et al._ (2023) T. V. Zache, D. González-Cuadra, and P. Zoller, arXiv preprint arXiv:2303.08683 (2023).
  * González-Cuadra  _et al._ (2023) D. González-Cuadra, D. Bluvstein, M. Kalinowski, R. Kaubruegger, N. Maskara, P. Naldesi, T. V. Zache, A. M. Kaufman, M. D. Lukin, H. Pichler, _et al._ , arXiv preprint arXiv:2303.06985 (2023).
  * Arute  _et al._ (2020) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, A. Bengtsson, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, Y.-A. Chen, B. Chiaro, R. Collins, S. J. Cotton, W. Courtney, S. Demura, A. Derk, A. Dunsworth, D. Eppens, T. Eckl, C. Erickson, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, J. A. Gross, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. Huggins, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, M. Marthaler, O. Martin, J. M. Martinis, A. Marusczyk, S. McArdle, J. R. McClean, T. McCourt, M. McEwen, A. Megrant, C. Mejuto-Zaera, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Newman, M. Y. Niu, T. E. O’Brien, E. Ostby, B. Pató, A. Petukhov, H. Putterman, C. Quintana, J.-M. Reiner, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, P. Schmitteckert, M. Szalay, N. M. Tubman, A. Vainsencher, T. White, N. Vogt, Z. J. Yao, P. Yeh, A. Zalcman, and S. Zanker, “Observation of separated dynamics of charge and spin in the fermi-hubbard model,” (2020), arXiv:2010.07965 [quant-ph] .
  * Barends  _et al._ (2015) R. Barends, L. Lamata, J. Kelly, L. García-Álvarez, A. G. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, I.-C. Hoi, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, E. Solano, and J. M. Martinis, Nature Communications 6, 7654 (2015).
  * Salathé  _et al._ (2015) Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff, Phys. Rev. X 5, 021027 (2015).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Stanisic  _et al._ (2022) S. Stanisic, J. L. Bosse, F. M. Gambetta, R. A. Santos, W. Mruczkiewicz, T. E. O’Brien, E. Ostby, and A. Montanaro, Nature Communications 13, 5743 (2022).
  * Jordan and Wigner (1928) P. Jordan and E. P. Wigner, Z. Phys 47, 14 (1928).
  * Bravyi and Kitaev (2002) S. B. Bravyi and A. Y. Kitaev, Annals of Physics 298, 210 (2002).
  * Fradkin  _et al._ (1980) E. Fradkin, M. Srednicki, and L. Susskind, Phys. Rev. D 21, 2885 (1980).
  * Srednicki (1980) M. Srednicki, Phys. Rev. D 21, 2878 (1980).
  * Verstraete and Cirac (2005) F. Verstraete and J. I. Cirac, Journal of Statistical Mechanics: Theory and Experiment 2005, P09012 (2005).
  * Chen  _et al._ (2018) Y.-A. Chen, A. Kapustin, and Đ. Radičević, Annals of Physics 393, 234 (2018).
  * Kitaev (2006) A. Kitaev, Annals of Physics 321, 2 (2006), january Special Issue.
  * Zohar and Cirac (2018) E. Zohar and J. I. Cirac, Phys. Rev. B 98, 075119 (2018).
  * Zohar and Cirac (2019) E. Zohar and J. I. Cirac, Phys. Rev. D 99, 114511 (2019).
  * Kitaev (2003a) A. Y. Kitaev, Annals of Physics 303, 2 (2003a).
  * Chen  _et al._ (2023) Y.-A. Chen, Y. Xu, _et al._ , PRX Quantum 4, 010326 (2023).
  * Pardo  _et al._ (2023) G. Pardo, T. Greenberg, A. Fortinsky, N. Katz, and E. Zohar, Phys. Rev. Res. 5, 023077 (2023).
  * Irmejs  _et al._ (2023a) R. Irmejs, M. C. Banuls, and J. I. Cirac, “Quantum simulation of z2 lattice gauge theory with minimal resources,” (2023a), arXiv:2206.08909 [quant-ph] .
  * Irmejs  _et al._ (2023b) R. Irmejs, M.-C. Bañuls, and J. I. Cirac, Phys. Rev. D 108, 074503 (2023b).
  * Vijayan  _et al._ (2020) J. Vijayan, P. Sompet, G. Salomon, J. Koepsell, S. Hirthe, A. Bohrdt, F. Grusdt, I. Bloch, and C. Gross, Science 367, 186 (2020).
  * Lieb (2013) E. H. Lieb,  _Condensed matter physics and exactly soluble models: selecta of Elliott H. Lieb_ (Springer Science & Business Media, 2013).
  * Corboz  _et al._ (2010) P. Corboz, R. Orús, B. Bauer, and G. Vidal, Phys. Rev. B 81, 165104 (2010).
  * Zohar and Burrello (2015) E. Zohar and M. Burrello, Phys. Rev. D 91, 054506 (2015).
  * Kitaev (2003b) A. Kitaev, Annals of Physics 303, 2 (2003b).
  * Kitaev and Laumann (2009) A. Kitaev and C. Laumann, Exact methods in low-dimensional statistical physics and quantum computing , 101 (2009).
  * Hubbard (1963) J. Hubbard, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Arovas  _et al._ (2021) D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual Review of Condensed Matter Physics 13 (2021), 10.1146/annurev-conmatphys-031620-102024.
  * Lieb and Wu (1968) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
  * Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Müller-Hartmann (1989) E. Müller-Hartmann, Zeitschrift für Physik B Condensed Matter 74, 507 (1989).
  * Qin  _et al._ (2021) M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13 (2021), 10.1146/annurev-conmatphys-090921-033948.
  * LeBlanc  _et al._ (2015) J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).
  * Lieb (1989) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * Silvi  _et al._ (2014) P. Silvi, E. Rico, T. Calarco, and S. Montangero, New Journal of Physics 16, 103015 (2014).
  * Cataldi  _et al._ (2021) G. Cataldi, A. Abedi, G. Magnifico, S. Notarnicola, N. D. Pozza, V. Giovannetti, and S. Montangero, Quantum 5, 556 (2021).
  * Ferrari  _et al._ (2022) G. Ferrari, G. Magnifico, and S. Montangero, Phys. Rev. B 105, 214201 (2022).
  * Cataldi  _et al._ (2023) G. Cataldi, G. Magnifico, P. Silvi, and S. Montangero, “(2+1)D SU(2) Yang-Mills Lattice Gauge Theory at finite density via tensor networks,”  (2023), arxiv:2307.09396 [cond-mat, physics:hep-lat, physics:hep-th, physics:quant-ph] .
  * Silvi  _et al._ (2019) P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, SciPost Phys. Lect. Notes , 8 (2019).
  * Montangero (2018) S. Montangero,  _Introduction to Tensor Network Methods_ (Springer International Publishing, 2018).
  * Corboz (2016) P. Corboz, Phys. Rev. B 93, 045116 (2016).
  * Whitfield  _et al._ (2016) J. D. Whitfield, V. c. v. Havlíček, and M. Troyer, Phys. Rev. A 94, 030301 (2016).
  * Steudtner and Wehner (2019) M. Steudtner and S. Wehner, Phys. Rev. A 99, 022308 (2019).
  * Setia  _et al._ (2019) K. Setia, S. Bravyi, A. Mezzacapo, and J. D. Whitfield, Phys. Rev. Res. 1, 033033 (2019).
  * Jiang  _et al._ (2019) Z. Jiang, J. McClean, R. Babbush, and H. Neven, Phys. Rev. Appl. 12, 064041 (2019).
  * Derby  _et al._ (2021) C. Derby, J. Klassen, J. Bausch, and T. Cubitt, Phys. Rev. B 104, 035118 (2021).
  * Kollath  _et al._ (2005) C. Kollath, U. Schollwöck, and W. Zwerger, Phys. Rev. Lett. 95, 176401 (2005).
  * Yang  _et al._ (2022) L. Yang, I. Hamad, L. O. Manuel, and A. E. Feiguin, Phys. Rev. B 105, 195104 (2022).
  * Chen  _et al._ (1994) Y. C. Chen, A. Moreo, F. Ortolani, E. Dagotto, and T. K. Lee, Phys. Rev. B 50, 655 (1994).
  * Wang  _et al._ (2022) P. Wang, G. Yu, Y. H. Kwan, Y. Jia, S. Lei, S. Klemenz, F. A. Cevallos, R. Singha, T. Devakul, K. Watanabe, T. Taniguchi, S. L. Sondhi, R. J. Cava, L. M. Schoop, S. A. Parameswaran, and S. Wu, Nature 605, 57 (2022).
  * (82) “Quantum matcha tea,” https://baltig.infn.it/quantum_matcha_tea/py_api_quantum_matcha_tea.
  * Gottesman (1997) D. Gottesman,  _Stabilizer codes and quantum error correction_ (California Institute of Technology, 1997).
  * Felser  _et al._ (2021) T. Felser, S. Notarnicola, and S. Montangero, Phys. Rev. Lett. 126, 170603 (2021).
  * Jordan  _et al._ (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
  * Jaschke and Montangero (2022) D. Jaschke and S. Montangero, arXiv preprint arXiv:2205.12092 (2022).