# Quantum Gauge Networks:
A New Kind of Tensor Network

Kevin Slagle  Department of Electrical and Computer Engineering, Rice
University, Houston, Texas 77005 USA  Department of Physics, California
Institute of Technology, Pasadena, California 91125, USA  Institute for
Quantum Information and Matter and Walter Burke Institute for Theoretical
Physics, California Institute of Technology, Pasadena, California 91125, USA

###### Abstract

Although tensor networks are powerful tools for simulating low-dimensional
quantum physics, tensor network algorithms are very computationally costly in
higher spatial dimensions. We introduce _quantum gauge networks_ : a different
kind of tensor network ansatz for which the computation cost of simulations
does not explicitly increase for larger spatial dimensions. We take
inspiration from the gauge picture of quantum dynamics Slagle (2022), which
consists of a local wavefunction for each patch of space, with neighboring
patches related by unitary connections. A quantum gauge network (QGN) has a
similar structure, except the Hilbert space dimensions of the local
wavefunctions and connections are truncated. We describe how a QGN can be
obtained from a generic wavefunction or matrix product state (MPS). All
$2k$-point correlation functions of any wavefunction for $M$ many operators
can be encoded exactly by a QGN with bond dimension $O(M^{k})$. In comparison,
for just $k=1$, an exponentially larger bond dimension of $2^{M/6}$ is
generically required for an MPS of qubits. We provide a simple QGN algorithm
for approximate simulations of quantum dynamics in any spatial dimension. The
approximate dynamics can achieve exact energy conservation for time-
independent Hamiltonians, and spatial symmetries can also be maintained
exactly. We benchmark the algorithm by simulating the quantum quench of
fermionic Hamiltonians in up to three spatial dimensions.

Tensor network algorithms Orús (2019, 2014a); Kin-Lic Chan  _et al._ ; Cirac
_et al._ ; Ran  _et al._ (2020); Bridgeman and Chubb (2017) are very useful
for simulating strongly-correlated quantum physics. In one spatial dimension,
matrix product state (MPS) algorithms Orús (2014a); Kin-Lic Chan  _et al._ ;
Cirac  _et al._ are often the best available tool for this task. Although
still useful, tensor network algorithms in higher dimensions Orús (2014a);
Cirac  _et al._ ; Zaletel and Pollmann (2020); Hyatt and Stoudenmire ;
Haghshenas  _et al._ (2019); Tepaske and Luitz (2021); Vidal (2008); Evenbly
and Vidal (2014, 2009); Acuaviva  _et al._ (2022); Ferrari  _et al._ (2022)
typically suffer from computational costs that scale as a high power of the
bond dimension. Therefore, we are motivated to study a different kind of
tensor network ansatz that is more computationally efficient in many spatial
dimensions.

We take inspiration from the gauge picture of quantum dynamics Slagle (2022).
In the gauge picture, one first chooses a collection of possibly-overlapping
patches of space that cover space. For example, one could choose the patches
to be pairs of nearest-neighbor sites on a lattice. We use capital letters,
$I$, $J$, or $K$, to denote a spatial patch. The Hamiltonian is written as a
sum over terms $\hat{H}_{I}$ that act within a single patch $I$:

$\hat{H}=\sum_{I}\hat{H}_{I}$ (1)

Each patch is assigned a _local wavefunction_ $\ket{\Psi_{I}}$, and the
Hilbert spaces of neighboring patches are related by unitary connections
$\hat{U}_{IJ}$. The local wavefunctions and connections time-evolve according
to

$\displaystyle\partial_{t}\ket{\Psi_{I}}$ $\displaystyle=-i\hat{H}_{\langle
I\rangle}\ket{\Psi_{I}}$ (2) $\displaystyle\partial_{t}\hat{U}_{IJ}$
$\displaystyle=-i\hat{H}_{\langle
I\rangle}\hat{U}_{IJ}+i\hat{U}_{IJ}\hat{H}_{\langle J\rangle}$

where

$\hat{H}_{\langle I\rangle}=\sum_{J}^{J\cap
I\neq\emptyset}\hat{U}_{IJ}\,\hat{H}_{J}\,\hat{U}_{JI}$ (3)

is the sum of local Hamiltonian terms supported on patches that overlap with
patch $I$. The main feature of the gauge picture is that spatial locality is
explicit in the above equations of motion. Typically, we initialize
$\hat{U}_{IJ}(0)=\hat{\mathbbl{1}}$ and $\ket{\Psi_{I}(0)}=\ket{\Psi(0)}$ at
time $t=0$, where $\hat{\mathbbl{1}}$ is the identity operator and
$\ket{\Psi(t)}$ is the usual wavefunction in the Schrödinger picture. The
expectation value $\braket{\Psi}{\hat{A}_{I}}{\Psi}$ of a local operator
$\hat{A}_{I}$ that only acts within the patch $I$ can be evaluated in the
gauge picture as $\braket{\Psi_{I}}{\hat{A}_{I}}{\Psi_{I}}$. To calculate an
expectation value $\braket{\Psi}{\hat{A}_{I}\hat{B}_{J}}{\Psi}$ for a product
of operators acting on different patches, a connection $\hat{U}_{IJ}$ must be
inserted in the gauge picture, as in
$\braket{\Psi_{I}}{\hat{A}_{I}\hat{U}_{IJ}\hat{B}_{J}}{\Psi_{J}}$. Time
evolution preserves the following network of relations:

$\displaystyle\hat{U}_{IJ}\ket{\Psi_{J}}$ $\displaystyle=\ket{\Psi_{I}}$ (4)
$\displaystyle\hat{U}_{IJ}\hat{U}_{JK}$ $\displaystyle=\hat{U}_{IK}$

along with $\hat{U}_{IJ}^{\dagger}=\hat{U}_{JI}$ and
$\hat{U}_{II}=\hat{\mathbbl{1}}$.

In this work, we truncate the Hilbert spaces of the local wavefunctions and
connections so that approximate quantum simulations can be performed on a
computer. The utility of the gauge picture for approximating quantum mechanics
is that locality is explicit in both the time dynamics and the structure of
the local wavefunctions and connections. We call the resulting network of
truncated local wavefunctions and connections a _quantum gauge network_ (QGN).

An advantage of quantum gauge networks is that unlike traditional tensor
networks (e.g. PEPS Orús (2014a)), a QGN only involves matrices and vectors
regardless of the spatial dimension. As such, many QGN calculations only
require a computation time (e.g. CPU time) that scales as $O(\chi^{3})$,
regardless of the number of spatial dimensions, where $\chi$ is the dimension
of the truncated Hilbert spaces. $\chi$ can be viewed as the bond dimension of
the QGN. The $O(\chi^{3})$ computation time for a QGN is the same as MPS
algorithms, which are very efficient in one spatial dimension (1D). In three
spatial dimensions (3D), the most computationally efficient tensor network in
previous literature may be the isometric tensor network Zaletel and Pollmann
(2020); Hyatt and Stoudenmire ; Haghshenas  _et al._ (2019), for which the
computation time of the TEBD3 algorithm is $O(\chi^{12})$ Tepaske and Luitz
(2021).111More generally, the TEBD3 computation time is
$O(D^{10}\chi^{2})+O(D\chi^{6})$ when the bond dimension $\chi$ of the central
bonds is different from the bond dimension $D$ of the other bonds. Tepaske and
Luitz (2021) In 2D, the TEBD2 computation time scales as $\chi^{7}$. Zaletel
and Pollmann (2020) Thus, quantum gauge network algorithms can require
significantly less computation time for fixed bond dimension. This is useful
since larger bond dimensions allow more correlations to be transported through
the tensor network. The computation time per variational parameter also scales
favorable for QGN algorithms.222The number of variational parameters can
differ significantly for different tensor networks. Therefore, it is also
useful to compare the scaling of the computation time per variational
parameter. For a QGN, the number of parameters scales as $\chi^{2}$ (due to
the matrix-valued connections). Therefore, the computation time per parameter
scales as $\chi^{3}/\chi^{2}=\chi^{1.5}$ for a QGN. This is the same ratio as
MPS algorithms, which are very efficient in 1D. For the isometric tensor
network TEBD3 algorithm in 3D, this ratio scales as
$\chi^{12}/\chi^{6}=\chi^{2}$ Tepaske and Luitz (2021), which is remarkably
efficient but not as good as the ratio for a QGN or MPS.

Although most tensor networks typically directly encode a wavefunction or
density matrix, a QGN departs from this habit. Indeed, in the gauge picture,
only local correlation functions $\braket{\Psi_{I}}{\hat{A}_{I}}{\Psi_{I}}$
can be extracted from the local wavefunctions if we do not make use of the
connections $\hat{U}_{IJ}$. Long range entanglement is then encoded within the
connections. Unlike a generic PEPS Orús (2014a) but similar to a matrix
product state (MPS) or isometric tensor network Zaletel and Pollmann (2020);
Hyatt and Stoudenmire ; Haghshenas  _et al._ (2019), local expectation values
can be efficiently computed from a QGN. A disadvantage of quantum gauge
networks is that unphysical states can also be encoded. As such, using a QGN
to variationally optimize a ground state is not as straight-forward as for
tensor networks that directly encode a wavefunction. We leave QGN ground state
optimization algorithms to future work. In this work, we focus on QGN
fundamentals and time dynamics alorithms.

In Sec. I, we discuss basic properties of quantum gauge networks and how a QGN
can be constructed. We also show that for an arbitrary wavefunction (including
fermionic wavefunctions), all $2k$-point correlation functions of $M$ many
operators can be encoded exactly by a QGN with bond dimension $O(M^{k})$ [Eq.
(30)], while an MPS of qubits can require an exponentially larger bond
dimension $2^{M/6}$ for $k=1$. In Sec. II, we present a QGN algorithm for
approximately simulating quantum dynamics in the gauge picture. We benchmark
the algorithm using simulations of fermionic Hamiltonians in spatial
dimensions up to three.

##  I Quantum Gauge Networks

To define a quantum gauge network (QGN), we first choose a collection of
possibly-overlapping patches of space that cover space. For example, one could
choose the patches to be nearest-neighbor sites on a lattice, or each patch
could consist of just a single lattice site. A QGN then consists of (1) a
_local wavefunction_ $\ket{\psi_{I}}$ for each spatial patch; (2) non-unitary
_connections_ $V_{IJ}=V_{JI}^{\dagger}$ to relate the Hilbert spaces of nearby
patches, as depicted in Fig. 1; and (3) a collection of truncated operators to
act on the truncated Hilbert space at each patch.

Figure 1: An example of a chain of qubits (black dots) and spatial patches
(colored ovals) consisting of pairs of neighboring qubits. A local
wavefunction $\ket{\psi_{I}}$ is associated with each patch $I$, and the
Hilbert spaces of neighboring patches are related by connections $V_{IJ}$.

The local wavefunctions and connections are similar to the those within the
gauge picture Slagle (2022), except the Hilbert space is truncated. If the
full Hilbert space has dimension $N$, then $\ket{\Psi_{I}}$ and $\hat{U}_{IJ}$
in the gauge picture have dimensions $N$ and $N\times N$, respectively. Since
$N$ is exponentially large in system size, it is useful to truncate the full
Hilbert space dimension for approximate simulations. Therefore, we consider
truncated local wavefunctions $\ket{\psi_{I}}$ and connections $V_{IJ}$, which
have truncated dimensions $\chi_{I}$ and $\chi_{I}\times\chi_{J}$,
respectively, where typically $\chi_{I}\ll N$. We use capital and lower-case
Greek letters (e.g. $\ket{\Psi_{I}}$ vs $\ket{\psi_{I}}$) for wavefunctions in
the full and truncated Hilbert spaces, respectively. Similarly, we place hats
on operators that act within the full Hilbert space (e.g. $\hat{U}_{IJ}$),
while operators within the truncated Hilbert space (e.g. $V_{IJ}$) do not have
hats.

In order to calculate expectation values of local operators $\hat{A}_{I}$ in
the original Hilbert space, we must also define truncated operators, i.e.
$\chi_{I}\times\chi_{I}$ matrices $A_{I}$, that act on the truncated Hilbert
space. Throughout this work, $\hat{A}_{I}$ always denotes an operator that
acts within a patch $I$, and similar for $\hat{B}_{J}$, etc. The truncated
operators $A_{I}$ are notationally distinguished from the original operators
$\hat{A}_{I}$ by the lack of a hat. In Sec. I.1, we present a concrete mapping
to obtain a QGN and truncated operators. However, in many cases (e.g. Appendix
C.1.1 and E) the truncated operators can be taken to be a simple Kronecker
product, such as $\sigma_{I}^{\mu}=\mathbbl{1}\otimes\sigma^{\mu}$, where
$\mathbbl{1}$ is an identity matrix and $\sigma^{\mu}$ is a $2\times 2$ Pauli
matrix.

Ideally, we want the quantum gauge network to accurately encode approximate
expectation values. For example, if the QGN is an approximation for a
wavefunction $\ket{\Psi}$, then we would like the QGN to accurately encode
local expectation values; i.e. we want
$\braket{\psi_{I}}{A_{I}}{\psi_{I}}\approx\braket{\Psi}{\hat{A}_{I}}{\Psi}$.
Similarly, we typically also want expectation values of string operators to
also approximately match, e.g.

$\braket{\psi_{I}}{A_{I}V_{IJ}B_{J}V_{JK}C_{K}}{\psi_{K}}\approx\braket{\Psi}{\hat{A}_{I}\hat{B}_{J}\hat{C}_{K}}{\Psi}$
(5)

Note that in order to express a string operator that acts on multiple spatial
patches using a QGN, it is essential to insert connections $V_{IJ}$ between
operators and wavefunctions associated with different spatial patches.

Similar to the gauge picture of quantum dynamics, quantum gauge networks also
exhibit a local gauge symmetry:

$\displaystyle\ket{\psi_{I}}$ $\displaystyle\to\Lambda_{I}\ket{\Psi_{I}}$ (6)
$\displaystyle V_{IJ}$
$\displaystyle\to\Lambda_{I}V_{IJ}\Lambda_{J}^{\dagger}$ $\displaystyle A_{I}$
$\displaystyle\to\Lambda_{I}A_{I}\Lambda_{I}^{\dagger}$

where $\Lambda_{I}$ is a unitary matrix. Expectation values must be invariant
under this symmetry.

Since local expectation values
$\braket{\Psi}{\hat{A}_{I}}{\Psi}\approx\braket{\psi_{I}}{A_{I}}{\psi_{I}}$
are encoded in the local wavefunctions $\ket{\psi_{I}}$, we can think of the
local wavefunction as a purified reduced density matrix for a spatial patch.
The connections $V_{IJ}$ encode long-range correlations between the spatial
patches.

It is desirable for a QGN to at least approximately obey the following
consistency conditions

$\displaystyle V_{IJ}\ket{\psi_{J}}$ $\displaystyle\approx\ket{\psi_{I}}$ (7)
$\displaystyle V_{IJ}V_{JK}$ $\displaystyle\approx V_{IK}$

Although it is easy to make the first relation exact, the second will
typically only hold approximately. Typically, a QGN will only possess a
$V_{IJ}$ for nearby patches $I$ and $J$ (and not for far away patches). Thus,
the second relation only applies if all three connections ($V_{IJ}$, $V_{JK}$,
and $V_{IK}$) are contained in the QGN. The connections $V_{IJ}$ should have
singular values less than or equal to 1 (to ensure that expectation values are
never larger than the largest eigenvalue of the measured operator).

When $V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ holds exactly, QGN connected
correlation functions obey the usual identity (with a $V_{IJ}$ inserted):

$\displaystyle\big{\langle}\psi_{I}\big{|}\big{(}A_{I}-\braket{A_{I}}\\!\big{)}V_{IJ}\big{(}B_{J}-\braket{B_{J}}\\!\big{)}\big{|}\psi_{J}\big{\rangle}$
(8) $\displaystyle=$
$\displaystyle\braket{\psi_{I}}{A_{I}V_{IJ}B_{J}}{\psi_{J}}-\braket{A_{I}}\braket{B_{J}}$

where we abbreviate $\braket{A_{I}}=\braket{\psi_{I}}{A_{I}}{\psi_{I}}$ and
$\braket{B_{J}}=\braket{\psi_{J}}{B_{J}}{\psi_{J}}$. Therefore, if the QGN
connected correlation function [Eq. (8)] is small, we are guaranteed that
$\braket{\psi_{I}}{A_{I}V_{IJ}B_{J}}{\psi_{J}}\approx\braket{A_{I}}\braket{B_{J}}$,
as one should expect. The equations in this paragraph also hold for longer
chains of connections; e.g. they also hold if we replace $V_{IJ}$ with
$V_{IK}V_{KL}V_{LI}$.

###  I.1 QGN from Truncation Maps

In this subsection, we study a concrete construction to obtain a quantum gauge
network. The input for this QGN construction is a wavefunction $\ket{\Psi}$,
along with a _truncation map_ $Q_{I}$ for each patch of space. If we want to
construct a QGN from a density matrix instead, then $\ket{\Psi}$ should be
chosen to be a purification of the density matrix. The truncation maps are
$\chi_{I}\times N$ matrices that satisfy:

$\displaystyle Q_{I}Q_{I}^{\dagger}$ $\displaystyle=\hat{\mathbbl{1}}$ (9)
$\displaystyle Q_{I}^{\dagger}Q_{I}\ket{\Psi}$ $\displaystyle=\ket{\Psi}$

Therefore $Q_{I}^{\dagger}$ is an isometry matrix whose image includes the
wavefunction. (A matrix $M$ is isometric if $M^{\dagger}M=\mathbbl{1}$.)
Intuitively, each $Q_{I}$ maps a select subspace of states into a truncated
Hilbert space. In the next subsection, we will explain how one can obtain
useful truncation maps.

With this data, we can construct the following QGN:

$\displaystyle\ket{\psi_{I}}$ $\displaystyle=Q_{I}\ket{\Psi}$ (10)
$\displaystyle V_{IJ}$ $\displaystyle=Q_{I}Q_{J}^{\dagger}$

Note that $V_{II}=\mathbbl{1}$ due to Eq. (9). Local operators $\hat{A}_{I}$
with support on a spatial patch $I$ are also truncated:

$A_{I}=Q_{I}\hat{A}_{I}Q_{I}^{\dagger}$ (11)

Note that this equation can be used for both bosonic and fermionic operators.
In the limit that the bond dimension $\chi_{I}\to N$ approaches the full
Hilbert space dimension $N$, this truncation mapping will result in a QGN that
encodes correlation functions [e.g. Eq. (5)] exactly.

The operator norm of $V_{IJ}$ is bounded by
$||V_{IJ}||_{\text{op}}\leq||Q_{I}||_{\text{op}}||Q_{J}^{\dagger}||_{\text{op}}=1$.
Thus, the resulting connections $V_{IJ}$ have singular values that are less
than or equal to 1. One can verify that

$V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ (12)

from Eq. (7) holds exactly.

In practice, the truncation mapping can only be directly applied for
wavefunctions that are simple enough such that the truncation can be
efficiently computed. However, the construction is also useful for
theoretically understanding how a QGN can encode a wavefunction.

A trivial example of a QGN can be obtained from a product state wavefunction
$\ket{\Psi}=\ket{\psi_{1}}\otimes\ket{\psi_{2}}\otimes\cdots\ket{\psi_{n}}$.
The truncation maps can be chosen to be $Q_{I}=\ket{\psi_{I}}\bra{\Psi}$, with
$I=1,\ldots,n$. This results in a QGN with local wavefunctions
$\ket{\psi_{I}}$ and connections $V_{IJ}=\ket{\psi_{I}}\otimes\bra{\psi_{J}}$.
The truncated operators simply act within the on-site Hilbert space.
Additional examples of quantum gauge networks can be found in Appendix C.

Note that not all quantum gauge networks can be obtained from the truncation
mapping. For example, generating a QGN by sampling random numbers for
$\ket{\psi_{I}}$ and $V_{IJ}$ will result in a very unphysical QGN that is not
consistent with any wavefunction. This is in contrast to MPS or PEPS tensor
networks, for which a random number initialization still returns a physical
(although unnormalized) wavefunction.

In Appendix B, we show how to obtain truncation maps from the canonical form
of a matrix product state (MPS) Orús (2014a); Kin-Lic Chan  _et al._ ; Cirac
_et al._ . If the MPS has bond dimension $\chi$, then the QGN will have bond
dimensions equal to $d\chi^{2}$, where $d$ is the Hilbert space dimension at
each site. ($d=2$ for qubits.) The idea of the mapping is that the canonical
form of an MPS can consist of a center tensor at $I$ that is surrounded by
isometric tensors. The isometric tensors are then used to construct a
truncation map $Q_{I}$, while the center tensor is a local wavefunction
$\ket{\psi_{I}}$.

###  I.2 Truncation Map Construction

The accuracy of the truncation critically depends on a good choice of
truncation maps $Q_{I}$. Suppose we want to choose truncation maps such that
the quantum gauge network exactly encodes the expectation values of a chosen
collection of operator strings. For example, the expectation value of
$\hat{A}_{I}\hat{B}_{J}\hat{C}_{K}$ is encoded exactly if

$\braket{\psi_{I}}{A_{I}V_{IJ}B_{J}V_{JK}C_{K}}{\psi_{K}}=\braket{\Psi}{\hat{A}_{I}\hat{B}_{J}\hat{C}_{K}}{\Psi}$
(13)

Below, we show that this exact encoding can be achieved by choosing truncation
maps such that the image of each $Q_{I}^{\dagger}$ is the span of certain
strings of operators involved in the expectation values.

In general, we will find that the required bond dimension $\chi_{I}$ for each
patch $I$ is bounded by $\chi_{I}\leq 1+2p_{I}$ [Eq. (21)], where $p_{I}$ is
the number of chosen operator strings (which we want to encode) that act on
the patch $I$. We will also find that bond dimension $O(M^{k})$ [Eq. (30)] is
sufficient to exactly encode all $2k$-point correlation functions of $M$
different operators.

####  I.2.1 Warmup Example

Before presenting a generic algorithm for obtaining the truncation maps, let
us first discuss an instructive example. Suppose we want to ensure that the
QGN encodes the expectation value in Eq. (13) exactly for a particular choice
of local operators ($\hat{A}_{I}$, $\hat{B}_{J}$, and $\hat{C}_{K}$) and
spatial patches $I\neq J\neq K$. Below, we show that any choice of truncation
maps with the following images is sufficient:

$\displaystyle\text{im}(Q_{I}^{\dagger})$
$\displaystyle=\text{span}\big{\\{}\ket{\Psi},\;\hat{A}_{I}^{\dagger}\ket{\Psi}\big{\\}}$
(14) $\displaystyle\text{im}(Q_{J}^{\dagger})$
$\displaystyle=\text{span}\big{\\{}\ket{\Psi},\;\hat{A}_{I}^{\dagger}\ket{\Psi},\;\hat{C}_{K}\ket{\Psi}\big{\\}}$
$\displaystyle\text{im}(Q_{K}^{\dagger})$
$\displaystyle=\text{span}\big{\\{}\ket{\Psi},\;\hat{C}_{K}\ket{\Psi}\big{\\}}$

$\text{span}\big{\\{}\ket{\Psi_{1}},\ldots,\ket{\Psi_{m}}\big{\\}}$ denotes
the vector space spanned by the vectors
$\ket{\Psi_{1}},\ldots,\ket{\Psi_{m}}$. We set the bond dimensions to be equal
to the vector space dimension of the images:
$\chi_{I}=\text{dim}(\text{im}(Q_{I}^{\dagger}))$. Specifying these images
determines the truncation maps $Q_{I}$ up to a unitary gauge transformation
$Q_{I}\to\Lambda_{I}Q_{I}$, where $\Lambda_{I}$ is a unitary matrix. The
choice of gauge does not affect QGN expectation values. For examples of a
quantum gauge networks obtained in this way, see Appendix C.

On a computer, a $Q_{I}^{\dagger}$ with the desired image can be calculated
from the compact singular value decomposition
$M_{I}=Q_{I}^{\dagger}S_{I}R_{I}$. Here, $M_{I}$ is a matrix of column vector,
which each encode one of the wavefunctions contained in the span of
$\text{im}(Q_{I}^{\dagger})$. $S_{I}$ is a $\chi_{I}\times\chi_{I}$ diagonal
matrix of nonzero singular values, and $R_{I}^{\dagger}$ is an isometry
matrix.

Note that if $\ket{\Phi}\in\text{im}(Q_{I}^{\dagger})$, then
$Q_{I}^{\dagger}Q_{I}\ket{\Phi}=\ket{\Phi}$, which is a useful property. This
follows because $\ket{\Phi}\in\text{im}(Q_{I}^{\dagger})$ implies that there
exists $\ket{\phi}$ such that $Q_{I}^{\dagger}\ket{\phi}=\ket{\Phi}$, which
then implies that
$Q_{I}^{\dagger}Q_{I}\ket{\Phi}=Q_{I}^{\dagger}Q_{I}Q_{I}^{\dagger}\ket{\phi}=Q_{I}^{\dagger}\ket{\phi}=\ket{\Phi}$
where the middle equality follows from
$Q_{I}Q_{I}^{\dagger}=\hat{\mathbbl{1}}$ [Eq. (9)].

For each patch, the image must contain $\ket{\Psi}$ so that Eq. (9) is
satisfied. Next, we show that Eq. (14) is sufficient to exactly encode the
expectation value in Eq. (13):

$\displaystyle\braket{\psi_{I}}{A_{I}V_{IJ}B_{J}V_{JK}C_{K}}{\psi_{K}}$
$\displaystyle=$
$\displaystyle\braket{\Psi}{Q_{I}^{\dagger}\,(Q_{I}\hat{A}_{I}Q_{I}^{\dagger})\,V_{IJ}B_{J}V_{JK}\,(Q_{K}\hat{C}_{K}Q_{K}^{\dagger})\,Q_{K}}{\Psi}$
$\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{A}_{I}Q_{I}^{\dagger}V_{IJ}B_{J}V_{JK}Q_{K}\hat{C}_{K}}{\Psi}$
$\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{A}_{I}Q_{I}^{\dagger}\,(Q_{I}Q_{J}^{\dagger})\,B_{J}\,(Q_{J}Q_{K}^{\dagger})\,Q_{K}\hat{C}_{K}}{\Psi}$
$\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{A}_{I}Q_{J}^{\dagger}B_{J}Q_{J}\hat{C}_{K}}{\Psi}$
$\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{A}_{I}Q_{J}^{\dagger}(Q_{J}\hat{B}_{J}Q_{J}^{\dagger})Q_{J}\hat{C}_{K}}{\Psi}$
(15) $\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{A}_{I}\hat{B}_{J}\hat{C}_{K}}{\Psi}$

We used the identities $\ket{\psi_{I}}=Q_{I}\ket{\Psi}$ and
$V_{IJ}=Q_{I}Q_{J}^{\dagger}$ [Eq. (10)],
$A_{I}=Q_{I}\hat{A}_{I}Q_{I}^{\dagger}$ [Eq. (11)], and
$Q_{I}^{\dagger}Q_{I}\ket{\Phi}=\ket{\Phi}$ whenever
$\ket{\Phi}\in\text{im}(Q_{I}^{\dagger})$ in Eq. (14).

Note that even if all the local operators commute, the path of the string
matters. For example, Eq. (14) does not guarantee that
$\braket{\psi_{J}}{B_{J}V_{JI}A_{I}V_{IK}C_{K}}{\psi_{K}}=\braket{\Psi}{\hat{B}_{J}\hat{A}_{I}\hat{C}_{K}}{\Psi}$
even if $[\hat{A}_{I},\hat{B}_{J}]=0$.

If we want the QGN to encode expectation values for multiple operator strings,
then we must calculate the images for each operator string, and then take the
union. For example, if two operator strings respectively require
$\text{im}(Q_{I}^{\dagger})\supseteq S_{I}^{(1)}$ and
$\text{im}(Q_{I}^{\dagger})\supseteq S_{I}^{(2)}$, then we need
$\text{im}(Q_{I}^{\dagger})\supseteq\text{span}(S_{I}^{(1)},S_{I}^{(2)})$.

Equation (14) is not the unique choice for the images. For instance, the
images could be larger. That is, it is sufficient to replace the equalities
“=” in Eq. (14) with superset relations “$\supseteq$”. Alternatively, we could
replace $\text{im}(Q_{J}^{\dagger})$ with

$\text{im}(Q_{J}^{\dagger})=\text{span}\big{\\{}\ket{\Psi},\;\hat{B}_{J}\hat{C}_{K}\ket{\Psi},\;\hat{C}_{K}\ket{\Psi}\big{\\}}$
(16)

And with this $\text{im}(Q_{J}^{\dagger})$, we could then replace
$\text{im}(Q_{I}^{\dagger})$ with

$\text{im}(Q_{I}^{\dagger})=\text{span}\big{\\{}\ket{\Psi},\;\hat{B}_{J}\hat{C}_{K}\ket{\Psi}\big{\\}}$
(17)

The proof is similar to Eq. (15). By comparing these choices, we see that we
have some freedom in choosing a midpoint in the string of operator products,
which we make explicit in the generic case below.

####  I.2.2 Generic Case

More generally, suppose we want to ensure that a QGN exactly encodes the
expectation value of
$\hat{A}_{I_{1}}^{(1)}\hat{A}_{I_{2}}^{(2)}\cdots\hat{A}_{I_{M}}^{(M)}$, which
is a product of $M$ local operators. This expectation value is exactly encoded
in the QGN if

$\displaystyle\big{\langle}\psi_{I_{1}}\big{|}A_{I_{1}}^{(1)}V_{I_{1},I_{2}}A_{I_{2}}^{(2)}V_{I_{2},I_{3}}\cdots
A_{I_{M}}^{(M)}\big{|}\psi_{I_{M}}\big{\rangle}$ (18) $\displaystyle=$
$\displaystyle\big{\langle}\Psi\big{|}\hat{A}_{I_{1}}^{(1)}\hat{A}_{I_{2}}^{(2)}\cdots\hat{A}_{I_{M}}^{(M)}\big{|}\Psi\big{\rangle}$

where $I_{1},\ldots,I_{m}$ is a string of neighboring patches with $I_{m}\neq
I_{m+1}$. The exact encoding [Eq. (18)] is guaranteed by any choice of
truncation maps with images that contain:

$\text{im}(Q_{I_{m}}^{\dagger})\supseteq\begin{cases}\text{span}\big{\\{}\ket{\Psi},\ket{\Psi_{m-1}^{\text{L}}},\ket{\Psi_{m}^{\text{L}}}\big{\\}}&m<m_{0}\\\
\text{span}\big{\\{}\ket{\Psi},\ket{\Psi_{m-1}^{\text{L}}},\ket{\Psi_{m+1}^{\text{R}}}\big{\\}}&m=m_{0}\\\
\text{span}\big{\\{}\ket{\Psi},\ket{\Psi_{m}^{\text{R}}},\ket{\Psi_{m+1}^{\text{R}}}\big{\\}}&m>m_{0}\end{cases}$
(19)

where we recursively define

$\displaystyle\big{|}\Psi_{0}^{\text{L}}\big{\rangle}$
$\displaystyle=\big{|}\Psi_{M+1}^{\text{R}}\big{\rangle}=\big{|}\Psi\big{\rangle}$
(20) $\displaystyle\big{|}\Psi_{m}^{\text{L}}\big{\rangle}$
$\displaystyle=\hat{A}_{I_{m}}^{(m)\dagger}\big{|}\Psi_{m-1}^{\text{L}}\big{\rangle}$
$\displaystyle\big{|}\Psi_{m}^{\text{R}}\big{\rangle}$
$\displaystyle=\hat{A}_{I_{m}}^{(m)}\big{|}\Psi_{m+1}^{\text{R}}\big{\rangle}$

In Eq. (19), we are free to choose any half-integer
$m_{0}=1,\tfrac{3}{2},\ldots,M$ between 1 and $M$.

As noted previously, if we want the QGN to encode expectation values for $M$
different operator strings, then for each image we obtain a list of conditions
$\text{im}(Q_{I}^{\dagger})\supseteq S_{I}^{(m)}$ from Eq. (19) with
$m=1,2,\ldots,M$. We can then take the images to be the span of these vector
spaces: $\text{im}(Q_{I}^{\dagger})=\text{span}(S_{I}^{(1)},S_{I}^{(2)},\ldots
S_{I}^{(M)})$.

We can easily bound the necessary bond dimensions $\chi_{I}$ needed to exactly
encode the expectation values for many operator strings. Let $p_{I}$ be the
number of operator strings that involve the spatial patch $I$. Each
application of Eq. (19) increases the dimension of
$\text{im}(Q_{I}^{\dagger})$ by at most 2 [starting from $\chi_{I}=1$ since we
always have $\ket{\Psi}\in\text{im}(Q_{I_{m}}^{\dagger})$]. Therefore, this
procedure results in bond dimensions of at most

$\chi_{I}\leq 1+2p_{I}$ (21)

To prove that Eq. (19) implies Eq. (18), we first recall that
$Q_{I}^{\dagger}Q_{I}\ket{\Phi}=\ket{\Phi}$ whenever
$\ket{\Phi}\in\text{im}(Q_{I}^{\dagger})$. We thus obtain

$\displaystyle\big{\langle}\psi_{I_{1}}\big{|}A_{I_{1}}^{(1)}V_{I_{1},I_{2}}\cdots
A_{I_{{\lceil m_{0}\rceil}-1}}^{({\lceil m_{0}\rceil}-1)}Q_{{\lceil
m_{0}\rceil}-1}$ (22) $\displaystyle=$
$\displaystyle\big{\langle}\Psi_{I_{1}}\big{|}\hat{A}_{I_{1}}^{(1)}\cdots\hat{A}_{I_{{\lceil
m_{0}\rceil}-1}}^{({\lceil m_{0}\rceil}-1)}$

and

$\displaystyle Q_{{\lfloor m_{0}\rfloor}+1}^{\dagger}A_{I_{{\lfloor
m_{0}\rfloor}+1}}^{({\lfloor m_{0}\rfloor}+1)}\cdots
V_{I_{M-1},I_{M}}A_{I_{M}}^{(M)}\big{|}\psi_{I_{M}}\big{\rangle}$ (23)
$\displaystyle=$ $\displaystyle\hat{A}_{I_{{\lfloor
m_{0}\rfloor}+1}}^{({\lfloor
m_{0}\rfloor}+1)}\cdots\hat{A}_{I_{M}}^{(M)}\big{|}\Psi\big{\rangle}$

by again using the identities $\ket{\psi_{I}}=Q_{I}\ket{\Psi}$,
$V_{IJ}=Q_{I}Q_{J}^{\dagger}$, and $A_{I}=Q_{I}\hat{A}_{I}Q_{I}^{\dagger}$. If
$m_{0}$ is an integer, then inserting
$V_{m_{0}-1,m_{0}}A_{m_{0}}V_{m_{0},m_{0}+1}=Q_{m_{0}-1}Q_{m_{0}}^{\dagger}A_{m_{0}}Q_{m_{0}}Q_{m_{0}+1}^{\dagger}$
and Eqs. (22) and (23) into the first line of Eq. (18) yields

$\displaystyle\big{\langle}\psi_{I_{1}}\big{|}A_{I_{1}}^{(1)}V_{I_{1},I_{2}}\cdots
A_{I_{M}}^{(M)}\big{|}\psi_{I_{M}}\big{\rangle}$ (24) $\displaystyle=$
$\displaystyle\big{\langle}\Psi\big{|}\hat{A}_{I_{1}}^{(1)}\cdots\hat{A}_{I_{m_{0}-1}}^{({m_{0}-1})}\big{(}Q_{m_{0}}^{\dagger}Q_{m_{0}}\hat{A}_{m_{0}}Q_{m_{0}}^{\dagger}Q_{m_{0}}\big{)}$
$\displaystyle\hskip
28.45274pt\hat{A}_{I_{m_{0}+1}}^{({m_{0}+1})}\cdots\hat{A}_{I_{M}}^{(M)}\big{|}\Psi\big{\rangle}$
$\displaystyle=$
$\displaystyle\big{\langle}\Psi\big{|}\hat{A}_{I_{1}}^{(1)}\cdots\hat{A}_{I_{M}}^{(M)}\big{|}\Psi\big{\rangle}$

If $m_{0}$ is a half-integer, then Eq. (18) follows immediately from
$V_{{\lfloor m_{0}\rfloor},{\lceil m_{0}\rceil}}=Q_{{\lfloor
m_{0}\rfloor}}^{\dagger}Q_{{\lceil m_{0}\rceil}}$ and Eqs. (22) and (23). This
completes the proof.

####  I.2.3 Long-Range Correlation Functions

We can now show that quantum gauge networks can efficiently encode long-range
correlation functions. Suppose we want a QGN to exactly encode all two-point
correlation functions for some arbitrary collection of operators
$\hat{\tau}_{i}^{\mu}$ (indexed by $\mu$ and position $i$). That is, suppose
we want a QGN to satisfy

$\big{\langle}\psi_{I}\big{|}\tau_{i\in
I}^{\mu\dagger}V_{IJ}^{\text{string}}\tau_{j\in
J}^{\nu}\big{|}\psi_{J}\big{\rangle}=\big{\langle}\Psi\big{|}\hat{\tau}_{i}^{\mu\dagger}\hat{\tau}_{j}^{\nu}\big{|}\Psi\big{\rangle}$
(25)

Patches $I$ and $J$ can be any patches that respectively contain sites $i$ and
$j$. $\tau_{i\in I}^{\mu}=Q_{I}\hat{\tau}_{i}^{\mu}Q_{I}^{\dagger}$ denotes a
truncated operator on patch $I$. Since patches $I$ and $J$ could be far apart,
we insert a string of connections

$V_{IJ}^{\text{string}}=V_{IK_{1}}V_{K_{1}K_{2}}\cdots V_{K_{l}J}$ (26)

to connect patches $I$ and $J$, where $(I,K_{1},K_{2},\ldots,K_{l},J)$ is an
arbitrary string of neighboring patches. To exactly encode the above
correlation functions, Eq. (19) implies that it is sufficient for the images
$\text{im}(Q_{I}^{\dagger})$ to include the span of the actions of each
operator on the wavefunction:

$\text{im}(Q_{I}^{\dagger})\supseteq\text{span}\big{\\{}\ket{\Psi},\hat{\tau}_{j}^{\mu}\ket{\Psi}\text{
for each }\hat{\tau}_{j}^{\mu}\big{\\}}$ (27)

If there are $M$ many operators $\hat{\tau}_{i}^{\mu}$, then the above image
requires a QGN with bond dimension of at most $\chi_{I}=1+M$.

More generally, suppose we want to exactly encode all $2k$-point correlation
functions:

$\displaystyle\big{\langle}\Psi\big{|}\hat{\tau}_{i_{1}}^{\mu_{1}\dagger}\cdots\hat{\tau}_{i_{k}}^{\mu_{k}\dagger}\hat{\tau}_{j_{1}}^{\nu_{1}}\cdots\hat{\tau}_{j_{k}}^{\nu_{k}}\big{|}\Psi\big{\rangle}$
(28) $\displaystyle=\big{\langle}\psi_{I_{1}}\big{|}\tau_{i_{1}\in
I_{1}}^{\mu_{1}\dagger}V_{I_{1},I_{2}}^{\text{string}}\cdots\tau_{i_{k}\in
I_{k}}^{\mu_{k}\dagger}V_{I_{k},J_{1}}^{\text{string}}$ $\displaystyle\hskip
34.14322pt\tau_{j_{1}\in
J_{1}}^{\nu_{1}}V_{J_{1},J_{2}}^{\text{string}}\cdots\tau_{j_{k}\in
J_{k}}^{\nu_{k}}\big{|}\psi_{J_{k}}\big{\rangle}$

The left-hand-side is the expectation value of a product of $2k$ operators.
The right-hand-side can be any corresponding expectation value within a QGN
for any valid strings of connections. To exactly encode these correlation
functions, Eq. (19) implies that it is sufficient for the images to include
the span of the actions of all products of up to $k$ operators:

$\text{im}(Q_{I}^{\dagger})\supseteq\text{span}\left\\{\text{all
}\textstyle\prod_{r=1}^{R}\hat{\tau}_{j_{r}}^{\mu_{r}}\ket{\Psi}\text{ for
}R=0,1,\ldots,k\right\\}$ (29)

If there are $M$ many operators $\hat{\tau}_{i}^{\mu}$, then there are

$\chi=1+M+M^{2}+\cdots+M^{k}$ (30)

different operator products included in the span. Therefore, a QGN with bond
dimension of at most $\chi$ is sufficient to encode all $2k$-point correlation
functions of $M$ many operators.

Although the bond dimension increases exponentially with $k$, typically only
few-body operators can be measured in experiments. Therefore, encoding
$k$-point correlation functions with large $k$ may not be necessary. On the
other hand, for large systems, the number of operators $M$ is large, and it is
advantageous that $\chi$ only increases as a polynomial for large $M$ and
fixed $k$.

For a matrix product state (MPS), encoding all 2-point correlation functions
generically requires an exponentially large bond dimension
$\chi^{\text{MPS}}=2^{M/6}$ for $M$ many operators on a chain of qubits. For
example, this exponential scaling is required when the wavefunction is a
rainbow state. However, a matrix product operator (MPO) with bond dimension
$O(M)$ is sufficient to encode all 2-point correlation functions. See Appendix
C.5 for details. Thus, QGN bond dimension scaling can be similar to that of an
MPO.

See appendices C.3 and C.4 for analytical expressions of quantum gauge
networks that exactly encode normal-ordered correlation functions of coherent
boson wavefunctions and fermionic Slater determinant wavefunctions.

##  II Time Evolution Algorithm

Since the quantum gauge network is closely related to the gauge picture Slagle
(2022) of quantum dynamics, we can straight-forwardly modify the gauge picture
to obtain an approximate algorithm for simulating quantum dynamics using a
QGN. To do this, we simply replace $\ket{\Psi_{I}}\to\ket{\psi_{I}}$,
$\hat{U}_{IJ}\to V_{IJ}$, and $\hat{H}_{\langle I\rangle}\to H^{\prime}_{I}$
(defined below) in the gauge picture equations of motion [Eq. (2)] to obtain:

$\displaystyle\partial_{t}\ket{\psi_{I}}$
$\displaystyle=-iH^{\prime}_{I}\ket{\psi_{I}}$ (31)
$\displaystyle\partial_{t}V_{IJ}$
$\displaystyle=-iH^{\prime}_{I}V_{IJ}+iV_{IJ}H^{\prime}_{J}$

If the Hamiltonian $\hat{H}=\sum_{I}\hat{H}_{I}$ [Eq. (1)] is a sum of local
terms $\hat{H}_{I}$ each supported on a spatial patch $I$, then we define

$H^{\prime}_{I}=\sum_{J}^{J\cap I\neq\emptyset}V_{IJ}\,H_{J}\,V_{JI}$ (32)

$\sum_{J}^{J\cap I\neq\emptyset}$ sums over all patches $J$ that have
nontrivial overlap with the spatial patch $I$. $H^{\prime}_{I}$ is analogous
to $\hat{H}_{\langle I\rangle}$ in the gauge picture [Eq. (3)]. Recall that
$H_{I}$ is a $\chi_{I}\times\chi_{I}$ matrix in the truncated Hilbert space,
while $\hat{H}_{I}$ is $N\times N$ and acts on the full $N$-dimensional
Hilbert space. $H_{I}$ can be obtained using the truncation mapping (11),
which also holds when $\hat{H}_{I}$ is time-dependent.

For time-independent Hamiltonians, the QGN approximation

$E_{\text{QGN}}(t)=\sum_{I}\braket{\psi_{I}(t)}{H_{I}}{\psi_{I}(t)}$ (33)

of the energy expectation value $\braket{\Psi(t)}{\hat{H}}{\Psi(t)}$ is
conserved up to numerical integration errors for the time evolution in Eq.
(31). See Appendix D for a proof.

For a QGN with bond dimension $\chi$, the computation time of this algorithm
is dominated by multiplication of $\chi\times\chi$ matrices. Therefore, the
computation time for each time step scales as $O(n_{\text{V}}\chi^{3})$, where
$n_{\text{V}}$ is the number of connections $V_{IJ}$ used by the QGN. The
memory cost scales as $O(n_{\text{V}}\chi^{2})$ for storing the
connections.333Let $\widetilde{\chi}_{IJ}$ be the number of non-zero singular
values of $V_{IJ}$. If
$\widetilde{\chi}_{IJ}\ll\text{min}(\chi_{I},\chi_{J})$, then it is more
efficient to decompose $V_{IJ}=V_{I}^{(J)}V_{J}^{(I)\dagger}$ where
$V_{I}^{(J)}$ is a $\chi_{I}\times\widetilde{\chi}_{J}$ isometric matrix. If
all $\chi_{I}=\chi$ and $\widetilde{\chi}_{IJ}=\widetilde{\chi}$ are equal,
and if the $H_{I}$ can be encoded as sparse matrices with $O(\chi)$ nonzero
entries (which is typical), then the simulation and memory cost can
respectively scale as $O(n_{\text{V}}\chi\widetilde{\chi}^{2})$ and
$O(n_{\text{V}}\chi\widetilde{\chi})$. This is asymptotically less costly when
$\widetilde{\chi}\ll\chi$. However for the simulations in this work, this
decomposition is not useful since we chose connections for which
$\widetilde{\chi}/\chi$ only ranges from about 0.75 (in one spatial dimension)
to 0.5 (in three spatial dimensions). For local Hamiltonians, $n_{\text{V}}$
should be proportional to the number of lattice sites.

It is also possible to handle Hamiltonian terms that are not supported on a
single spatial patch. Consider the Hamiltonian

$\hat{H}=\sum_{J\cdots K}\sum_{\mu\cdots\nu}h_{J\cdots
K}^{\mu\cdots\nu}\;\hat{\tau}^{\mu}_{J}\cdots\hat{\tau}^{\nu}_{K}$ (34)

which consists of a sum $\sum_{J\cdots K}$ over spatial patches $J\cdots K$.
$h_{J\cdots K}^{\mu\cdots\nu}$ are (generically time-dependent) real
coefficients, and $\hat{\tau}^{\mu}_{J}$ denotes an operator indexed by $\mu$
with support on patch $J$. For example, we could take the patches to consist
of a single qubit, and the $\hat{\tau}^{\mu}_{J}$ could be Pauli operators.
For local Hamiltonians, $h_{J\cdots K}^{\mu\cdots\nu}$ will only be nonzero if
$J\cdots K$ are close together. Equation (32) could be generalized to

$\displaystyle H^{\prime}_{I}=$ $\displaystyle\sum_{J\cdots K}^{(J\cup\cdots
K)\cap I\neq\emptyset}\sum_{\mu\cdots\nu}$ (35) $\displaystyle\hskip
22.76228pt\tfrac{1}{2}h_{J\cdots
K}^{\mu\cdots\nu}\;(V_{IJ}\tau^{\mu}_{J}V_{JI})\cdots(V_{IK}\tau^{\nu}_{K}V_{KI})+h.c.$

where $\sum_{J\cdots K}^{(J\cup\cdots K)\cap I\neq\emptyset}$ sums over
spatial patches $J\cdots K$ such that the union $J\cup\cdots K$ has nontrivial
overlap with patch $I$. “$h.c.$” denotes the Hermitian conjugate of the
preceding terms and ensures that $H^{\prime}_{I}$ is Hermitian. When $J\cdots
K$ involves more than two patches, we must decide on an ordering of the
$J\cdots K$ for each $I$. Eq. (35) follows from a similar generalized
expression for $\hat{H}_{\langle I\rangle}$ in the gauge picture Slagle (2022)
after projecting onto the Hermitian part (via the $h.c.$) and replacing
$\hat{U}_{IJ}\to V_{IJ}$ and $\hat{\tau}^{\mu}_{J}\to\tau^{\mu}_{J}$. However,
unlike for Eq. (32), the energy will not be conserved exactly for this choice
of $H^{\prime}_{I}$.

###  II.1 Fermion Quench

To benchmark this QGN algorithm, we simulate the dynamics of a quantum quench
and compare to exact methods. First, we study the quench dynamics for a model
of spinless fermions in one, two, and three spatial dimensions. In Appendix E,
we also study the quench to a near-critical transverse field Ising model on a
square lattice. In all cases, we find that increasing the bond dimension
increases the simulation accuracy. This is as expected, since exact
simulations of the gauge picture are reproduced once the bond dimension
reaches the full Hilbert space dimension (or less when there are conserved
quantities).

(a) 22 (b)

(c)

(d) $4\times 4$ (e)

(f)

Figure 2: Simulation data for the time dynamics of the fermionic Hamiltonian
(36) with $V=1$ following a quench from a checkerboard initial state in 1D and
2D periodic lattices. (a) An initial state of 11 spinless fermions (black
disks) on a periodic chain of 22 sites (black and white disks). (b) Starting
from this initial state, we show the fermion number expectation value
$\braket{\hat{n}_{i_{0}}(t)}$ vs time $t$ at a site $i_{0}$ with zero fermions
at time $t=0$. Simulation results are shown for the exact value (black line)
and the quantum gauge network (QGN) with different bond dimensions $\chi$
(colored lines). The legend also shows the number of CPU core hours used for
each simulation. (c) The error of the QGN approximation
$\braket{n_{i_{0}}}_{\text{QGN}}$ [Eq. (37)] to the exact expression
$\braket{\hat{n}_{i_{0}}}$. Panels (d-f) are similar, but for a periodic
$4\times 4$ square lattice. In both dimensions, we see that increasing the
bond dimension increases the accuracy of the QGN simulations.

(a) 22 (b)

(c)

(d) $4\times 4$ (e)

(f)

(g) $4\times 4\times 4$

(h)  (i)

Figure 3: Panels (a-f) are the same as Fig. 2, except we simulate non-
interacting fermions with $V=0$. By comparing to Fig. 2, we see that the QGN
simulates $V=0$ with roughly the same accuracy as $V=1$. Panels (g-i) are
similar, but for a periodic $4\times 4\times 4$ cubic lattice.

We initialize the system with a checkerboard pattern of spinless fermions, for
which $\braket{\hat{n}_{i}(0)}=1-\braket{\hat{n}_{j}(0)}$ at time $t=0$ for
nearest-neighbor sites $i$ and $j$; as shown in Figs. 2a and 2d. We then time-
evolve the system using the following Hamiltonian with nearest-neighbor
hoppings and nearest-neighbor repulsive interactions:

$\displaystyle\hat{H}^{\text{Fermi}}$ $\displaystyle=\sum_{\langle
ij\rangle}\hat{H}^{\text{Fermi}}_{\langle ij\rangle}$ (36)
$\displaystyle\hat{H}^{\text{Fermi}}_{I=\langle ij\rangle}$
$\displaystyle=-\hat{c}_{i}^{\dagger}\hat{c}_{j}-\hat{c}_{j}^{\dagger}\hat{c}_{i}+V\hat{n}_{i}\hat{n}_{j}$

Each patch $I$ is composed of a pair of nearest-neighbor sites, which are
summed over by $\sum_{\langle ij\rangle}$ in the first line. $\hat{c}_{i}$ is
a fermion annihilation operator, and
$\hat{n}_{i}=\hat{c}_{i}^{\dagger}\hat{c}_{i}$ is the fermion number operator.

We initialize the QGN using truncation maps, as described in Sec. I.1.
Choosing the truncation maps $Q_{I}$ requires choosing a subspace of states to
keep for each patch (i.e. a choice for the image of $Q_{I}^{\dagger}$, as in
Sec. I.2). We use the following method to select subspaces of states that we
expect will acquire the highest weight after a short time evolution. (1) We
begin with images of $Q_{I}^{\dagger}$ that only contain the initial state.
(2) For each patch $I$ of nearest-neighbor sites, we add states to the image
of $Q_{I}^{\dagger}$ that can be obtained from the current image by swapping
the two sites within the patch. (3) For each patch $I$, we add states that are
included in patches $J$ that overlap with patch $I$. (4) We repeat steps 2 and
3 until we reach a bond dimension that is as large as desired.

With this initialization algorithm, we do not choose the bond dimension $\chi$
precisely (only roughly). This has the advantage that $\chi$ is chosen such
that the QGN remains symmetric under lattice symmetries. Each image consists
of a span of certain eigenstates of the number operators $\hat{n}_{i}$. Note
that since the initial state has a definite fermion number, this procedure
only adds states with the same fermion number, which is desirable since the
Hamiltonian also conserves the fermion number. If we were to repeat steps 2
and 3 until no additional states could be added, then we would add all states
with the initial fermion number and the QGN simulation would be exact.

We compare the quantum gauge network results to the exact values.444Exact
expressions for small system sizes can be obtained by calculating the full
wavefunction $\ket{\Psi(t)}=e^{-i\hat{H}t}\ket{\Psi(0)}$ using a sparse matrix
representation for the Hamiltonian. Figure 2 shows one-dimensional (1D) and
two-dimensional (2D) simulation data with interaction strength $V=1$. In both
spatial dimensions, we find that the QGN can accurately simulate the dynamics
for a short time, and the time over which the simulations are accurate
increases with increasing bond dimension. The most time-consuming simulation
consumed 90 CPU hours, which took less than half a day on an 8-core laptop.

We calculate the number expectation value $\braket{\hat{n}_{i_{0}}(t)}$ for
sites $i_{0}$ that initially have zero fermions. In the QGN, these are
estimated as follows:

$\braket{n_{i}}_{\text{QGN}}=\sum_{I\ni
i}^{\text{mean}}\braket{\psi_{I}}{n_{i\in I}}{\psi_{I}}$ (37)

$\sum_{I\ni i}^{\text{mean}}$ averages over all patches $I$ that contain the
site $i$, and $n_{i\in I}=Q_{I}\hat{n}_{i}Q_{I}^{\dagger}$ is the truncated
[Eq. (11)] number operator at site $i$ for patch $I$. In this example,
$\braket{\psi_{I}}{n_{i\in I}}{\psi_{I}}$ is equal for all patches $I$ that
contain site $i$ due to spatial symmetries. However in other models with less
symmetry, simulation errors can make these expectation values differ for
different patches.

If we were to integrate the equations of motion exactly, then the energy
expectation value [Eq. (33)] would be conserved exactly. Since exact
integration is not practical, we use a modified RK4 Runge-Kutta method for
numerical integration with time step $\delta_{\text{t}}=0.05$. Due to this
approximation, the energy per site changed by at most $10^{-3}$ for all data
shown. See Appendix F for more details.

In order to compare to exact methods in three dimensions (3D) with many sites,
we repeat the comparison in Fig. 3 with no interactions ($V=0$). Free (i.e.
non-interacting $V=0$) fermion systems are efficient to simulate exactly Afo
and are typically about as challenging for tensor network methods as
interacting fermionic systems. By comparing the 1D and 2D data in Figs. 2 and
3, we indeed see that $V=1$ and $V=0$ appear to be roughly equally challenging
for the QGN. Therefore, we expect that the comparison we make for free
fermions in 3D is representative of the interacting $V=1$ model (for which we
could not perform exact simulations). We find that accurate simulations can be
maintained for a short period of time, which increases with the bond dimension
and decreases with increasing spatial dimension.

For the free fermion simulations, the energy expectation value and total
fermion number within the QGN appears to be conserved exactly (up to floating
point precision). We have not yet investigated why this conservation occurs in
the free fermion system. When interactions ($V\neq 0$) are included, the QGN
conserves the energy expectation value up to numerical integration errors, and
the QGN fails to conserve the total charge expectation value.

##  III Outlook

Quantum gauge networks offer several advantages for simulating quantum
dynamics, especially in comparison to many tensor network methods. (1) The
computation time for simulating time dynamics is only $\chi^{3}$, where $\chi$
is the bond dimension. Notably, this computation time does not increase with
the spatial dimension, which makes quantum gauge networks a promising tool in
two or more dimensions. (2) Fermionic models are simple to handle. (Unlike
MERA or PEPS, fermionic swap gate Orús (2014b); Corboz and Vidal (2009) are
not needed.) (3) The code is simple since the only tensors involved are
vectors $\psi_{I}$ and matrices $V_{IJ}$, and the code does not get more
complicated in larger spatial dimensions. (4) The energy expectation value can
be conserved (up to integration error). (5) Lattice symmetries can be
maintained exactly. (6) Time discretization errors are small in practice since
the dynamics can be integrated using very accurate Runge-Kutta methods.
(Trotter-Suzuki expansions Childs  _et al._ (2021); Vanhecke  _et al._ (2019)
are not needed.)

There are numerous important future directions for the study of quantum gauge
networks: (1) In comparison to other tensor networks, understanding quantum
gauge networks is conceptually more demanding since the wavefunction is not
directly encoded. The truncation mapping (10) is an example for which we can
understand how the QGN relates to a wavefunction. But we do not know how to
tell if a given QGN is consistent with any choice of wavefunction and
truncation maps.555This problem may be related to the QMA complete Liu (2006)
quantum marginal problem Klyachko (2006); Chen  _et al._ (2016); Mazziotti
(2012). We also do not know to what extent the truncation mapping can produce
all quantum gauge networks that are useful approximations for quantum
wavefunctions. (2) Can we optimize a QGN to find approximate ground states or
excited eigenstates? This is relatively challenging for quantum gauge networks
because a QGN can encode unphysical states, which means that some sort of
(possibly approximate) constraint must be imposed on the QGN during energy
minimization. (3) How well can a QGN encode topological states Wen (2017); Gu
_et al._ (2009); Buerschaper  _et al._ (2009); Williamson  _et al._ (2017);
Soejima  _et al._ (2020)?

There are also many opportunities to significantly improve our QGN time
evolution algorithm: (4) For the TEBD algorithm Vidal (2004), it is
straightforward to increase or optimally truncate the bond dimension during
the simulation. In this work, we initialized the quantum gauge network using a
simple basis of states in the number basis. We expect that this simple
initialization is far from optimal, and that dynamically adding and removing
more optimally chosen states throughout the time evolution could greatly
improve simulation accuracy. (5) Can we obtain exact energy conservation when
the Hamiltonian terms act on multiple spatial patches [as in Eq. (34)]? (6)
Can we obtain charge conservation for charge-conserving Hamiltonians? (7) The
TEBD algorithm allows one to upper-bound the simulation error in terms of the
truncation error. Can we also estimate the error of a QGN time evolution
without comparing to other algorithms? (8) Could we simulate infinite system
sizes when the Hamiltonian is translation-invariant. (9) Finally, after
improving QGN algorithms, benchmarking QGN methods against other methods
Paeckel  _et al._ (2019); Vidal (2004); White and Feiguin (2004); Haegeman
_et al._ (2016); Leviatan  _et al._ (2017); Mendl (2018); Czarnik  _et al._
(2019); Bauernfeind and Aichhorn (2020); Rakovszky  _et al._ (2022); Yang and
White (2020); Kloss  _et al._ (2020); Alhambra and Cirac (2021); Lin  _et al._
(2021); Schmitt and Heyl (2020); Gutiérrez and Mendl (2022); Lin and Pollmann
(2022); Yehorova and Kretchmer (2022) will be useful.

We find it intriguing that the consistency conditions
($V_{IJ}\ket{\psi_{J}}\approx\ket{\psi_{I}}$ and $V_{IJ}V_{JK}\approx V_{IK}$)
in Eq. (7) are precisely the equations for a classical lattice gauge theory
Münster and Walzl (2000); Kogut (1979) to be in its ground state when coupled
to a Higgs field, where $V_{IJ}$ plays the role of the gauge connection and
$\ket{\psi_{I}}$ is the Higgs field. (See Appendix A for more details.) The
primary difference is that in lattice gauge theory, the gauge connections are
typically chosen to be unitary matrices. But $V_{IJ}$ is not a unitary matrix;
its singular values are only constrained to be less than or equal to 1. With
unitary gauge connections, no information is encoded locally in the classical
ground state. (Information is only encoded in non-contractible Wilson loops.)
It is remarkable that by simply relaxing the unitary constraint on the gauge
connections (as in a QGN), grounds states of classical lattice gauge theory
coupled to a Higgs field are capable of locally encoding approximate quantum
wavefunctions. Gauge theory plays a foundational role within the standard
model of particle physics. As such, it may be interesting to study the
emergent physics of quantum gauge networks, viewed not as a computational
tool, but instead as a new kind of classical lattice gauge theory that
exhibits aspects of emergent quantum mechanics Slagle and Preskill (2022);
Aaronson (2004); ’t Hooft (2020); Adler (2012); Vanchurin (2019); Nelson
(2012); Hall  _et al._ (2014).

###### Acknowledgements.

We thank Gunhee Park, Lesik Motrunich, and Garnet Chan for helpful
conversations. K.S. was partially supported by the Walter Burke Institute for
Theoretical Physics at Caltech; and the U.S. Department of Energy, Office of
Science, National Quantum Information Science Research Centers, Quantum
Science Center.

## References

  * Slagle (2022) Kevin Slagle, “The Gauge Picture of Quantum Dynamics,” (2022), arXiv:arXiv:2210.09314 .
  * Orús (2019) Román Orús, “Tensor networks for complex quantum systems,” Nature Reviews Physics 1, 538–550 (2019), arXiv:1812.04011 .
  * Orús (2014a) Román Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Annals of Physics 349, 117–158 (2014a), arXiv:1306.2164 .
  * (4) Garnet Kin-Lic Chan, Anna Keselman, Naoki Nakatani, Zhendong Li, and Steven R. White, “Matrix Product Operators, Matrix Product States, and ab initio Density Matrix Renormalization Group algorithms,” arXiv:1605.02611 .
  * (5) Ignacio Cirac, David Perez-Garcia, Norbert Schuch, and Frank Verstraete, “Matrix Product States and Projected Entangled Pair States: Concepts, Symmetries, and Theorems,” arXiv:2011.12127 .
  * Ran  _et al._ (2020) Shi-Ju Ran, Emanuele Tirrito, Cheng Peng, Xi Chen, Luca Tagliacozzo, Gang Su, and Maciej Lewenstein, in  _Tensor Network Contractions_ (Springer, 2020) arXiv:1708.09213 .
  * Bridgeman and Chubb (2017) Jacob C. Bridgeman and Christopher T. Chubb, “Hand-waving and interpretive dance: an introductory course on tensor networks,” Journal of Physics A Mathematical General 50, 223001 (2017), arXiv:1603.03039 .
  * Zaletel and Pollmann (2020) Michael P. Zaletel and Frank Pollmann, “Isometric Tensor Network States in Two Dimensions,” Phys. Rev. Lett. 124, 037201 (2020), arXiv:1902.05100 .
  * (9) Katharine Hyatt and E. M. Stoudenmire, “DMRG Approach to Optimizing Two-Dimensional Tensor Networks,” arXiv:1908.08833 .
  * Haghshenas  _et al._ (2019) Reza Haghshenas, Matthew J. O’Rourke, and Garnet Kin-Lic Chan, “Conversion of projected entangled pair states into a canonical form,” Phys. Rev. B 100, 054404 (2019), arXiv:1903.03843 .
  * Tepaske and Luitz (2021) Maurits S. J. Tepaske and David J. Luitz, “Three-dimensional isometric tensor networks,” Physical Review Research 3, 023236 (2021), arXiv:2005.13592 .
  * Vidal (2008) G. Vidal, “Class of Quantum Many-Body States That Can Be Efficiently Simulated,” Phys. Rev. Lett. 101, 110501 (2008), arXiv:quant-ph/0610099 .
  * Evenbly and Vidal (2014) G. Evenbly and G. Vidal, “Class of Highly Entangled Many-Body States that can be Efficiently Simulated,” Phys. Rev. Lett. 112, 240502 (2014), arXiv:1210.1895 .
  * Evenbly and Vidal (2009) G. Evenbly and G. Vidal, “Algorithms for entanglement renormalization,” Phys. Rev. B 79, 144108 (2009), arXiv:0707.1454 .
  * Acuaviva  _et al._ (2022) Arturo Acuaviva, Visu Makam, Harold Nieuwboer, David Pérez-García, Friedrich Sittner, Michael Walter, and Freek Witteveen, “The minimal canonical form of a tensor network,” (2022), arXiv:2209.14358 .
  * Ferrari  _et al._ (2022) Giovanni Ferrari, Giuseppe Magnifico, and Simone Montangero, “Adaptive-weighted tree tensor networks for disordered quantum many-body systems,” Phys. Rev. B 105, 214201 (2022), arXiv:2111.12398 .
  * (17) The time dynamics of a free fermion Hamiltonian $\hat{H}=\sum_{ij}h_{ij}\,\hat{c}_{i}^{\dagger}\hat{c}_{j}$ can be simulated exactly by calculating the time-evolved filled single-fermion wavefunctions $\ket{\phi_{\alpha}(t)}=e^{-iht}\ket{\phi_{\alpha}(0)}$. The wavefunction $\ket{\Psi}=\prod_{\alpha}^{\text{filled}}\big{(}\sum_{i}\braket{i}{\phi_{\alpha}}\hat{c}_{i}^{\dagger}\big{)}\ket{0}$ is never explicitly calculated. $\prod_{\alpha}^{\text{filled}}$ denotes the product over the filled single-fermion wavefunctions, and $\ket{0}$ is the empty state with no fermions. Then $\braket{\hat{n}_{i}(t)}=\sum_{\alpha}^{\text{filled}}|\braket{i}{\phi_{\alpha}(t)}|^{2}$, where $\ket{i}$ is the single-fermion wavefunction for a fermion at site $i$.
  * Orús (2014b) Román Orús, “Advances on tensor network theory: symmetries, fermions, entanglement, and holography,” European Physical Journal B 87, 280 (2014b), arXiv:1407.6552 .
  * Corboz and Vidal (2009) Philippe Corboz and Guifré Vidal, “Fermionic multiscale entanglement renormalization ansatz,” Phys. Rev. B 80, 165129 (2009), arXiv:0907.3184 .
  * Childs  _et al._ (2021) Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu, “Theory of trotter error with commutator scaling,” Phys. Rev. X 11, 011020 (2021), arXiv:1912.08854 .
  * Vanhecke  _et al._ (2019) Bram Vanhecke, Laurens Vanderstraeten, and Frank Verstraete, “Symmetric cluster expansions with tensor networks,” (2019), arXiv:1912.10512 .
  * Liu (2006) Yi-Kai Liu, “Consistency of local density matrices is qma-complete,” in  _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques_ , edited by Josep Díaz, Klaus Jansen, José D. P. Rolim, and Uri Zwick (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 438–449, arXiv:quant-ph/0604166 .
  * Klyachko (2006) Alexander A. Klyachko, “Quantum marginal problem and N-representability,” in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 36 (2006) pp. 72–86, arXiv:quant-ph/0511102 .
  * Chen  _et al._ (2016) Jianxin Chen, Zhengfeng Ji, Nengkun Yu, and Bei Zeng, “Detecting consistency of overlapping quantum marginals by separability,” Phys. Rev. A 93, 032105 (2016), arXiv:1509.06591 .
  * Mazziotti (2012) David A. Mazziotti, “Structure of fermionic density matrices: Complete $n$-representability conditions,” Phys. Rev. Lett. 108, 263002 (2012), arXiv:1112.5866 .
  * Wen (2017) Xiao-Gang Wen, “Colloquium: Zoo of quantum-topological phases of matter,” Reviews of Modern Physics 89, 041004 (2017), arXiv:1610.03911 .
  * Gu  _et al._ (2009) Zheng-Cheng Gu, Michael Levin, Brian Swingle, and Xiao-Gang Wen, “Tensor-product representations for string-net condensed states,” Phys. Rev. B 79, 085118 (2009), arXiv:0809.2821 .
  * Buerschaper  _et al._ (2009) Oliver Buerschaper, Miguel Aguado, and Guifré Vidal, “Explicit tensor network representation for the ground states of string-net models,” Phys. Rev. B 79, 085119 (2009), arXiv:0809.2393 .
  * Williamson  _et al._ (2017) Dominic J. Williamson, Nick Bultinck, and Frank Verstraete, “Symmetry-enriched topological order in tensor networks: Defects, gauging and anyon condensation,” (2017), arXiv:1711.07982 .
  * Soejima  _et al._ (2020) Tomohiro Soejima, Karthik Siva, Nick Bultinck, Shubhayu Chatterjee, Frank Pollmann, and Michael P. Zaletel, “Isometric tensor network representation of string-net liquids,” Phys. Rev. B 101, 085117 (2020), arXiv:1908.07545 .
  * Vidal (2004) Guifré Vidal, “Efficient Simulation of One-Dimensional Quantum Many-Body Systems,” Phys. Rev. Lett. 93, 040502 (2004), arXiv:quant-ph/0310089 .
  * Paeckel  _et al._ (2019) Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, and Claudius Hubig, “Time-evolution methods for matrix-product states,” Annals of Physics 411, 167998 (2019), arXiv:1901.05824 .
  * White and Feiguin (2004) Steven R. White and Adrian E. Feiguin, “Real-Time Evolution Using the Density Matrix Renormalization Group,” Phys. Rev. Lett. 93, 076401 (2004), arXiv:cond-mat/0403310 .
  * Haegeman  _et al._ (2016) Jutho Haegeman, Christian Lubich, Ivan Oseledets, Bart Vandereycken, and Frank Verstraete, “Unifying time evolution and optimization with matrix product states,” Phys. Rev. B 94, 165116 (2016), arXiv:1408.5056 .
  * Leviatan  _et al._ (2017) Eyal Leviatan, Frank Pollmann, Jens H. Bardarson, David A. Huse, and Ehud Altman, “Quantum thermalization dynamics with Matrix-Product States,” (2017), arXiv:1702.08894 .
  * Mendl (2018) Christian B. Mendl, “Time evolution of matrix product operators with energy conservation,” (2018), arXiv:1812.11876 .
  * Czarnik  _et al._ (2019) Piotr Czarnik, Jacek Dziarmaga, and Philippe Corboz, “Time evolution of an infinite projected entangled pair state: An efficient algorithm,” Phys. Rev. B 99, 035115 (2019), arXiv:1811.05497 .
  * Bauernfeind and Aichhorn (2020) Daniel Bauernfeind and Markus Aichhorn, “Time dependent variational principle for tree Tensor Networks,” SciPost Physics 8, 024 (2020), arXiv:1908.03090 .
  * Rakovszky  _et al._ (2022) Tibor Rakovszky, C. W. von Keyserlingk, and Frank Pollmann, “Dissipation-assisted operator evolution method for capturing hydrodynamic transport,” Phys. Rev. B 105, 075131 (2022), arXiv:2004.05177 .
  * Yang and White (2020) Mingru Yang and Steven R. White, “Time-dependent variational principle with ancillary Krylov subspace,” Phys. Rev. B 102, 094315 (2020), arXiv:2005.06104 .
  * Kloss  _et al._ (2020) Benedikt Kloss, David Reichman, and Yevgeny Bar Lev, “Studying dynamics in two-dimensional quantum lattices using tree tensor network states,” SciPost Physics 9, 070 (2020), arXiv:2003.08944 .
  * Alhambra and Cirac (2021) Álvaro M. Alhambra and J. Ignacio Cirac, “Locally Accurate Tensor Networks for Thermal States and Time Evolution,” PRX Quantum 2, 040331 (2021), arXiv:2106.00710 .
  * Lin  _et al._ (2021) Sheng-Hsuan Lin, Michael Zaletel, and Frank Pollmann, “Efficient Simulation of Dynamics in Two-Dimensional Quantum Spin Systems with Isometric Tensor Networks,” (2021), arXiv:2112.08394 .
  * Schmitt and Heyl (2020) Markus Schmitt and Markus Heyl, “Quantum Many-Body Dynamics in Two Dimensions with Artificial Neural Networks,” Phys. Rev. Lett. 125, 100503 (2020), arXiv:1912.08828 .
  * Gutiérrez and Mendl (2022) Irene López Gutiérrez and Christian B. Mendl, “Real time evolution with neural-network quantum states,” Quantum 6, 627 (2022), arXiv:1912.08831 .
  * Lin and Pollmann (2022) Sheng-Hsuan Lin and Frank Pollmann, “Scaling of Neural‑Network Quantum States for Time Evolution,” Physica Status Solidi B Basic Research 259, 2100172 (2022), arXiv:2104.10696 .
  * Yehorova and Kretchmer (2022) Dariia Yehorova and Joshua S. Kretchmer, “A multi-fragment real-time extension of projected density matrix embedding theory: Non-equilibrium electron dynamics in extended systems,” (2022), arXiv:2209.06368 .
  * Münster and Walzl (2000) G. Münster and M. Walzl, “Lattice Gauge Theory - A short Primer,” (2000), arXiv:hep-lat/0012005 .
  * Kogut (1979) John B. Kogut, “An introduction to lattice gauge theory and spin systems,” Rev. Mod. Phys. 51, 659–713 (1979).
  * Slagle and Preskill (2022) Kevin Slagle and John Preskill, “Emergent Quantum Mechanics at the Boundary of a Local Classical Lattice Model,” (2022), arXiv:2207.09465 .
  * Aaronson (2004) Scott Aaronson, “Multilinear formulas and skepticism of quantum computing,” in  _Proceedings of the Thirty-Sixth Annual ACM Symposium on Theory of Computing_, STOC ’04 (Association for Computing Machinery, New York, NY, USA, 2004) p. 118–127, arXiv:quant-ph/0311039 .
  * ’t Hooft (2020) Gerard ’t Hooft, “Deterministic Quantum Mechanics: the Mathematical Equations,” (2020), arXiv:2005.06374 .
  * Adler (2012) Stephen L Adler, “Quantum theory as an emergent phenomenon: Foundations and phenomenology,” Journal of Physics: Conference Series 361, 012002 (2012).
  * Vanchurin (2019) Vitaly Vanchurin, “Entropic Mechanics: Towards a Stochastic Description of Quantum Mechanics,” Foundations of Physics 50, 40–53 (2019), arXiv:1901.07369 .
  * Nelson (2012) Edward Nelson, “Review of stochastic mechanics,” Journal of Physics: Conference Series 361, 012011 (2012).
  * Hall  _et al._ (2014) Michael J. W. Hall, Dirk-André Deckert, and Howard M. Wiseman, “Quantum Phenomena Modeled by Interactions between Many Classical Worlds,” Physical Review X 4, 041013 (2014), arXiv:1402.6144 .
  * Vidal (2007) G. Vidal, “Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension,” Phys. Rev. Lett. 98, 070201 (2007), arXiv:cond-mat/0605597 .
  * Garcia  _et al._ (2019) Stephan Ramon Garcia, Matthew Okubo Patterson, and William T. Ross, “Partially isometric matrices: a brief and selective survey,” (2019), arXiv:1903.11648 .
  * Hamer (2000) C. J. Hamer, “Finite-size scaling in the transverse Ising model on a square lattice,” Journal of Physics A Mathematical General 33, 6683–6698 (2000), arXiv:cond-mat/0007063 .

##  Appendix A Higgsed Lattice Gauge Theory

In this appendix, we review more details regarding the connection between
Higgsed lattice gauge theory and quangum gauge networks. Consider a lattice of
vertices connected by edges, e.g. a triangular lattice. Each vertex of the
lattice hosts a Higgs field $\psi_{i}\in\mathbb{C}^{N}$ with $|\psi_{i}|=1$,
i.e. a normalized complex vector with $N$ components. Each edge of the lattice
hosts a $\text{U}(N)$ gauge connection $U_{ij}$, i.e. an $N\times N$ unitary
matrix. The energy of a classical $\text{U}(N)$ lattice gauge theory coupled
to a Higgs field (on a lattice with triangular plaquettes) is

$E=-\sum_{ijk}^{\text{plaquettes}}\operatorname{tr}U_{ij}U_{jk}U_{ki}-\sum_{\langle
ij\rangle}\text{Re}(\psi_{i}^{*}\cdot U_{ij}\cdot\psi_{j})$ (38)

where $U_{ji}=U_{ij}^{\dagger}$ and $\sum_{ijk}^{\text{plaquettes}}$ sums over
all plaquettes of the lattice (in two or more dimensions).
$U_{ij}U_{jk}U_{ki}$ is a product of gauge fields around the edges of a
plaquette (which we assumed to be a triangle only for notational simplicity).
$\sum_{\langle ij\rangle}$ sums over vertices $i$ and $j$ connected by an
edge. $\text{Re}(\psi_{i}^{*}\cdot U_{ij}\cdot\psi_{j})$ denotes the real part
of $\psi_{i}^{*}\cdot U_{ij}\cdot\psi_{j}$.

The energy of this classical lattice gauge theory is minimized when
$U_{ij}U_{jk}=U_{ik}$ and $U_{ij}\cdot\psi_{j}=\psi_{i}$, which is analogous
to the consistency conditions in Eqs. (4) and (7). However, there is an
important difference in each case: (1) In the gauge picture, $\hat{U}_{IJ}$ is
an $N\times N$ unitary matrix where $N$ is the full Hilbert space dimension,
which increases exponentially with system size. However in lattice gauge
theory, $N$ is typically taken to be a fixed integer. (2) In a quantum gauge
network, $V_{IJ}$ is not a unitary matrix; its singular values are only
constrained to be less than or equal to 1.

##  Appendix B Matrix Product State Mapping

In this appendix, we show that any matrix product state (MPS) with bond
dimension $\chi$ can be mapped to a quantum gauge network with bond dimension
$d\chi^{2}$, where $d$ is the Hilbert space dimension at each site. ($d=2$ for
qubits). Before explaining the mapping, we first briefly review MPS canonical
forms.

###  B.1 MPS Review

An matrix product state is a compact representation of a wavefunction, where
the wavefunction amplitudes are given by matrix products. Orús (2014a); Kin-
Lic Chan  _et al._ ; Cirac  _et al._ An MPS is specified by a
$\chi_{i}^{\text{MPS}}\times\chi_{i+1}^{\text{MPS}}$ rectangular matrix
$M_{i}^{(s_{i})}$ for each site $i$ and local state $\ket{s_{i}}$.
Equivalently, if $s_{i}=1,2,\ldots,d_{i}$ can take on $d_{i}$ different
states, then each $M_{i}$ can be viewed as a $\chi_{i}^{\text{MPS}}\times
d_{i}\times\chi_{i+1}^{\text{MPS}}$ tensor. We restrict
$\chi_{1}^{\text{MPS}}=\chi_{n+1}^{\text{MPS}}=1$. The MPS wavefunction for a
chain of $n$ sites is

$\displaystyle\big{|}\Psi_{\text{MPS}}\big{\rangle}$
$\displaystyle=\sum_{s_{1}s_{2}\cdots
s_{n}}\operatorname{tr}\\!\big{(}M_{1}^{(s_{1})}M_{2}^{(s_{2})}\cdots
M_{n}^{(s_{n})}\big{)}\big{|}s_{1}s_{2}\cdots s_{n}\big{\rangle}$
$\displaystyle\begin{array}[]{c}\\\ {\scalebox{1.5}{=}}\\\
{\scriptstyle(n=5)}\end{array}\;\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$
(42)

The second line shows a tensor network diagram for $n=5$ sites. The blue
circles represent the tensors $M_{i}$; lines between tensors denote contracted
indices; dangling lines denote uncontracted indices (the states $\ket{s_{i}}$
in this case); and we suppress the bond dimension
$\chi_{1}^{\text{MPS}}=\chi_{n+1}^{\text{MPS}}=1$ lines that are traced out.

There is a gauge redundancy
$M_{i}^{(s_{i})}\to\Lambda_{i}M_{i}^{(s_{i})}\Lambda_{i+1}^{\dagger}$ (with
unitary $\Lambda_{i}$) between neighboring matrices that does not affect the
encoded wavefunction. This redundancy is often used to compute a transformed
MPS in a canonical form centered at a specific site $i$:

$\displaystyle\big{|}\Psi_{\text{MPS}}\big{\rangle}$
$\displaystyle=\sum_{s_{1}\cdots
s_{n}}\operatorname{tr}\\!\Big{(}L_{1}^{(s_{1})}\cdots
L_{i-1}^{(s_{i-1})}C_{i}^{(s_{i})}$ (43) $\displaystyle\hskip
52.63777ptR_{i+1}^{(s_{i+1})}\cdots R_{n}^{(s_{n})}\Big{)}\big{|}s_{1}\cdots
s_{n}\big{\rangle}$ $\displaystyle\begin{array}[]{c}{\scriptstyle(i=3)}\\\
{\scalebox{1.5}{=}}\\\
{\scriptstyle(n=5)}\end{array}\;\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$
(47)

$L_{i}$, $C_{i}$, and $R_{i}$ are obtained from $M_{i}$ using gauge
transformations such that the following identities are obeyed:

$\displaystyle=\sum_{s_{i}}L_{i}^{(s_{i})\dagger}L_{i}^{(s_{i})}=\mathbbl{1}$
(48)
$\displaystyle=\sum_{s_{i}}\operatorname{tr}\\!\big{(}C_{i}^{(s_{i})\dagger}C_{i}^{(s_{i})}\big{)}=1$
$\displaystyle=\sum_{s_{i}}R_{i}^{(s_{i})}R_{i}^{(s_{i})\dagger}=\mathbbl{1}$

Thus, the $L_{i}$ and $R_{i}$ tensors are isometries, which pick out an
orthonormal basis of states for the orthogonality center $C_{i}$, which is
normalized like a wavefunction (but in a truncated Hilbert space). One utility
of the canonical form is that local expectation values are easy to compute:

$\braket{\Psi_{\text{MPS}}}{\hat{A}_{i}}{\Psi_{\text{MPS}}}=\sum_{s_{i}s^{\prime}_{i}}\operatorname{tr}\\!\Big{(}C_{i}^{(s_{i})\dagger}C_{i}^{(s^{\prime}_{i})}\Big{)}\braket{s_{i}}{\hat{A}_{i}}{s^{\prime}_{i}}$
(49)

$\hat{A}_{i}$ is an operator that only acts on site $i$.

It is also possible to obtain multiple simultaneous canonical forms, one for
each $i=1,2,\ldots,n$, while sharing the same isometries ($L_{i}$ and
$R_{i}$). These shared isometries obey

$\displaystyle L_{i}^{(s_{i})}C_{i+1}^{(s_{i+1})}$
$\displaystyle=C_{i}^{(s_{i})}R_{i+1}^{(s_{i+1})}$ (50)
$\displaystyle=\;\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$

This simultaneous canonical form can be obtained by sweeping across the MPS
multiple times using SVD decompositions or from the Vidal gauge Vidal (2007).
The resulting bond dimensions of the $\chi_{i}^{\text{MPS}}\times
d_{i}\times\chi_{i+1}^{\text{MPS}}$ tensors obey $\chi_{i}^{\text{MPS}}\leq
d_{i}\chi_{i+1}^{\text{MPS}}$ and
$\chi_{i}^{\text{MPS}}d_{i}\geq\chi_{i+1}^{\text{MPS}}$.

###  B.2 Quantum Gauge Network from MPS

To construct a quantum gauge network from an MPS, we choose $Q_{i}^{\dagger}$
to map the truncated Hilbert space of the MPS orthogonality center $C_{i}$ to
the full Hilbert space:

$\displaystyle Q_{i}^{\dagger}$ $\displaystyle=\sum_{s_{1}\cdots
s_{n}}\sum_{\alpha,\beta}\Big{(}L_{1}^{(s_{1})}\cdots
L_{i-1}^{(s_{i-1})}\Big{)}_{1,\alpha}$ (51) $\displaystyle\hskip
62.59596pt\Big{(}R_{i+1}^{(s_{i+1})}\cdots
R_{n}^{(s_{n})}\Big{)}_{\beta,1}\ket{s_{1}\cdots s_{n}}\bra{\alpha
s_{i}\beta}$ $\displaystyle\begin{array}[]{c}{\scriptstyle(i=3)}\\\
{\scalebox{1.5}{=}}\\\
{\scriptstyle(n=5)}\end{array}\;\;\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$
(55)

Thus, $Q_{i}^{\dagger}$ is just the canonical MPS centered at $i$, but with
the orthogonality center $C_{i}$ removed. Here, we choose the spatial patches
to consists of just a single site, and we make no notational distinction
between capital $(I,J,K)$ and lower case $(i,j,k)$ spatial indices letters.

Using Eq. (10) and the MPS canonical form identities [Eq. (48)], we find that
the local wavefunctions are equal to MPS orthogonality centers, while the
connections are equal to tensor products of two MPS isometries:

$\displaystyle\ket{\psi_{i}}$
$\displaystyle=\sum_{\alpha,s_{i},\beta}\big{(}C_{i}^{(s_{i})}\big{)}_{\alpha,\beta}\,\ket{\alpha
s_{i}\beta}$
$\displaystyle=\;\,\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$
$\displaystyle V_{i,i+1}$
$\displaystyle=\sum_{\alpha,s_{i},\alpha^{\prime}}\sum_{\beta,s_{i+1},\beta^{\prime}}\big{(}L_{i}^{(s_{i})}\big{)}_{\alpha\alpha^{\prime}}\big{(}R_{i+1}^{(s_{i+1})}\big{)}^{*}_{\beta\beta^{\prime}}$
$\displaystyle\hskip 75.39963pt\ket{\alpha
s_{i}\beta}\bra{\alpha^{\prime}s_{i+1}\beta^{\prime}}$ (56)
$\displaystyle=\quad\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$

We use the simultaneous canonical forms, which obey Eq. (50) and guarantee
that $V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ [Eq. (12)] is satisfied.

The bond dimensions of the quantum gauge network are therefore
$\chi_{i}=\chi_{i}^{\text{MPS}}d_{i}\chi_{i+1}^{\text{MPS}}$. The singular
values of the $\chi_{i}\times\chi_{i+1}$ matrix $V_{i,i+1}$ consist of
$(\chi_{i+1}^{\text{MPS}})^{2}$ ones, while the rest are zero. Thus, these
$V_{i,i+1}$ are partial isometry matrices Garcia  _et al._ (2019), which are
matrices that obey $VV^{\dagger}V=V$ (or equivalently matrices whose singular
values are either zero or one).

##  Appendix C QGN Examples

In this appendix, we discuss several examples of quantum gauge networks.

###  C.1 Mixed State Example

As a simple example of a quantum gauge network, below we construct a QGN for
the following mixed state of $n$ qubits:

$\displaystyle\hat{\rho}^{\text{mix}}$
$\displaystyle=\tfrac{1}{2}\ket{\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}\bra{\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}$
(57)
$\displaystyle\,+\tfrac{1}{2}\ket{\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}\bra{\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}$

In order to apply the truncation mapping [Eq. (10)], we first find a
purification of the density matrix:

$\displaystyle\ket{\Psi^{\text{mix}}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\ket{\uparrow_{0}\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}$
(58)
$\displaystyle\,+\tfrac{1}{\sqrt{2}}\ket{\downarrow_{0}\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}$

The wavefunction $\ket{\Psi^{\text{mix}}}$ hosts an additional auxiliary qubit
with index $I=0$. In this example, we take the spatial patches to consist of a
single qubit. $\ket{\Psi^{\text{mix}}}$ is a purification of
$\hat{\rho}^{\text{mix}}$ because tracing out the $I=0$ qubit yields the
density matrix:
$\hat{\rho}^{\text{mix}}=\operatorname{tr}_{0}\ket{\Psi^{\text{mix}}}\bra{\Psi^{\text{mix}}}$.

We follow Sec. I.2 to derive truncation maps $Q_{I}$ by choosing the images
$\text{im}(Q_{I}^{\dagger})$ to be carefully chosen subspaces of states. To
this end, we first note that any operator is a linear combination of Pauli
strings, and the only Pauli strings with a nonzero expectation value for the
state $\hat{\rho}^{\text{mix}}$ are products of
$\hat{\sigma}_{I}^{\text{z}}\hat{\sigma}_{J}^{\text{z}}$ operators. We will
therefore ensure that expectation values of products of
$\hat{\sigma}_{I}^{\text{z}}\hat{\sigma}_{J}^{\text{z}}$ operators are
retained by the truncation. It turns out that this is sufficient to exactly
encode all correlation functions for this example. Equation (27) thus implies
that we only need to include the action of a single
$\hat{\sigma}_{I}^{\text{z}}$ in the image of $Q_{I}^{\dagger}$:

$\displaystyle\text{im}(Q_{I}^{\dagger})$
$\displaystyle=\text{span}\big{\\{}\ket{\Psi^{\text{mix}}},\;\hat{\sigma}_{I}^{\text{z}}\ket{\Psi^{\text{mix}}}\big{\\}}$
(59)
$\displaystyle=\text{span}\big{\\{}\ket{\uparrow_{0}\uparrow_{1}\cdots\uparrow_{n}},\;\ket{\downarrow_{0}\downarrow_{1}\cdots\downarrow_{n}}\big{\\}}$

This image has dimension $\chi_{I}=2$. A natural gauge choice for $Q_{I}$
consistent with the above is:

$\displaystyle Q_{I}$
$\displaystyle=\ket{\uparrow}\bra{\uparrow_{0}\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}$
(60)
$\displaystyle\,+\ket{\downarrow}\bra{\downarrow_{0}\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}$

Equation (10) then yields the following QGN:

$\displaystyle\ket{\psi_{I}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\ket{\uparrow}+\tfrac{1}{\sqrt{2}}\ket{\downarrow}$
(61) $\displaystyle V_{IJ}$
$\displaystyle=\ket{\uparrow}\bra{\uparrow}+\ket{\downarrow}\bra{\downarrow}$

with truncated operators [Eq. (11)]

$\displaystyle\sigma_{I}^{\text{x}}$ $\displaystyle=\sigma_{I}^{\text{y}}=0$
(62) $\displaystyle\sigma_{I}^{\text{z}}$
$\displaystyle=\ket{\uparrow}\bra{\uparrow}-\ket{\downarrow}\bra{\downarrow}$

This quantum gauge network exactly encodes all expectation values of the
original reduced density matrix. For example,

$\displaystyle\braket{\psi_{I}}{\sigma_{I}^{\mu}}{\psi_{I}}$
$\displaystyle=\operatorname{tr}\hat{\rho}^{\text{mix}}\,\hat{\sigma}_{I}^{\mu}=0$
(63)
$\displaystyle\braket{\psi_{I}}{\sigma_{I}^{\mu}\,V_{IJ}\,\sigma_{J}^{\nu}}{\psi_{J}}$
$\displaystyle=\operatorname{tr}\hat{\rho}^{\text{mix}}\,\hat{\sigma}_{I}^{\mu}\hat{\sigma}_{J}^{\nu}=\begin{cases}1&\mu=\nu=\text{z}\\\
0&\text{otherwise}\end{cases}$

####  C.1.1 Kronecker Product Operators

If we want to preserve the algebra of more of the truncated operators, then we
should include their action in the images. For example, if we want to preserve
the on-site algebra of the truncated Pauli operators, then we should instead
choose:

$\displaystyle\text{im}(Q_{I}^{\dagger})$
$\displaystyle=\text{span}\big{\\{}\ket{\Psi^{\text{mix}}},\;\hat{\sigma}_{I}^{\text{x}}\ket{\Psi^{\text{mix}}},\;\hat{\sigma}_{I}^{\text{y}}\ket{\Psi^{\text{mix}}},\;\hat{\sigma}_{I}^{\text{z}}\ket{\Psi^{\text{mix}}}\big{\\}}$
(64)
$\displaystyle=\text{span}\big{\\{}\ket{\uparrow_{0}\uparrow_{1}\cdots\uparrow_{n}},\;\hat{\sigma}_{I}^{\text{x}}\ket{\uparrow_{0}\uparrow_{1}\cdots\uparrow_{n}},$
$\displaystyle\hskip
39.26494pt\ket{\downarrow_{0}\downarrow_{1}\cdots\downarrow_{n}},\;\hat{\sigma}_{I}^{\text{x}}\ket{\downarrow_{0}\downarrow_{1}\cdots\downarrow_{n}}\big{\\}}$

We can then pick the following truncation map:

$\displaystyle Q_{I}$
$\displaystyle=\ket{\uparrow_{0}\uparrow_{I}}\bra{\uparrow_{0}\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}$
(65)
$\displaystyle\,+\ket{\uparrow_{0}\downarrow_{I}}\bra{\uparrow_{0}\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}\hat{\sigma}_{I}^{\text{x}}$
$\displaystyle\,+\ket{\downarrow_{0}\uparrow_{I}}\bra{\downarrow_{0}\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}\hat{\sigma}_{I}^{\text{x}}$
$\displaystyle\,+\ket{\downarrow_{0}\downarrow_{I}}\bra{\downarrow_{0}\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}$

The resulting QGN follows from Eq. (10):

$\displaystyle\ket{\psi_{I}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\ket{\uparrow_{0}\uparrow_{I}}+\tfrac{1}{\sqrt{2}}\ket{\downarrow_{0}\downarrow_{I}}$
(66) $\displaystyle V_{IJ}$
$\displaystyle=\ket{\uparrow_{0}\uparrow_{I}}\bra{\uparrow_{0}\uparrow_{J}}+\ket{\downarrow_{0}\downarrow_{I}}\bra{\downarrow_{0}\downarrow_{J}}$

Now the truncated Pauli operators are their natural Kronecker products:

$\displaystyle\sigma_{I}^{\text{x}}$
$\displaystyle=\mathbbl{1}_{0}\otimes\big{(}\quad\;\;\ket{\uparrow_{I}}\bra{\downarrow_{I}}+\;\;\ket{\uparrow_{I}}\bra{\downarrow_{I}}\big{)}$
(67) $\displaystyle\sigma_{I}^{\text{y}}$
$\displaystyle=\mathbbl{1}_{0}\otimes\big{(}-i\ket{\uparrow_{I}}\bra{\downarrow_{I}}+i\ket{\uparrow_{I}}\bra{\downarrow_{I}}\big{)}$
$\displaystyle\sigma_{I}^{\text{z}}$
$\displaystyle=\mathbbl{1}_{0}\otimes\big{(}\quad\;\;\ket{\uparrow_{I}}\bra{\uparrow_{I}}+\;\;\ket{\downarrow_{I}}\bra{\downarrow_{I}}\big{)}$

where
$\mathbbl{1}_{0}=\ket{\uparrow_{0}}\bra{\uparrow_{0}}+\ket{\downarrow_{0}}\bra{\downarrow_{0}}$.
These truncated operators obey their usual on-site algebra, e.g.
$\sigma_{I}^{\text{x}}\sigma_{I}^{\text{y}}=i\sigma_{I}^{\text{z}}$.

###  C.2 Cat State Example

Encoding expectation values that act on many qubits in more than one spatial
dimension is less straight-forward. For example, consider the following cat
state

$\displaystyle\ket{\Psi^{\text{cat}}}$
$\displaystyle=\tfrac{1}{\sqrt{2}}\ket{\uparrow_{1}\uparrow_{2}\cdots\uparrow_{n}}$
(68)
$\displaystyle\,+\tfrac{1}{\sqrt{2}}\ket{\downarrow_{1}\downarrow_{2}\cdots\downarrow_{n}}$

The only Pauli string expectation value that can distinguish
$\ket{\Psi^{\text{cat}}}$ from the mixed state $\hat{\rho}^{\text{mix}}$ in
Eq. (57) is the highly-nonlocal product of Pauli $\hat{\sigma}^{\text{x}}$
operators on every qubit:

$\displaystyle\operatorname{tr}\hat{\rho}^{\text{mix}}\prod_{I=1}^{n}\hat{\sigma}^{\text{x}}_{I}$
$\displaystyle=0$ (69)
$\displaystyle\bigg{\langle}\Psi^{\text{cat}}\bigg{|}\prod_{I=1}^{n}\hat{\sigma}^{\text{x}}_{I}\bigg{|}\Psi^{\text{cat}}\bigg{\rangle}$
$\displaystyle=1$

Any Pauli string that does not act on all qubits will have an equal
expectation value for $\hat{\rho}^{\text{mix}}$ and $\ket{\Psi^{\text{cat}}}$.

(a)  (b)

(c)

Figure 4: Three different paths that cover all lattice sites.

Now suppose that $\ket{\Psi^{\text{cat}}}$ is a wavefunction for a square
lattice of qubits. A quantum gauge network for $\ket{\Psi^{\text{cat}}}$
should reproduce the same nonlocal expectation value:

$\braket{\psi_{1}}{\sigma^{\text{x}}_{1}V_{12}\sigma^{\text{x}}_{2}\cdots
V_{n-1,n}\sigma^{\text{x}}_{n}}{\psi_{n}}=1$ (70)

where the sites $1,2,\ldots,n$ snake across the square lattice, as depicted in
Fig. 4a. However, one may want other choices of paths [e.g. Fig. 4b or 4c] for
this string operator to also lead to the same expectation value. This can be
achieved by adding these additional string operators to the procedure in Sec.
I.2 at the cost of increasing the bond dimension. But if these additional
string operators are not included in the QGN construction, then the
expectation value of these excluded strings will not be encoded correctly.
This example demonstrates the issue that a quantum gauge network can seem to
encode different values for the same nonlocal expectation value depending on
the path chosen.

###  C.3 Bosonic Coherent States

The normal ordered expectation values of bosonic coherent states can be
encoded within a quantum gauge network in a rather trivial way. A bosonic
coherent state is specified by complex numbers $\Theta_{i}$ and takes the
following form:

$\ket{\Theta}=\exp\\!\bigg{(}\\!\sum_{i}\Theta_{i}\hat{b}_{i}^{\dagger}\bigg{)}\ket{0}$
(71)

$\hat{b}_{i}$ is a boson annihilation operator, which satisfies the
commutation relations $[\hat{b}_{i},\hat{b}_{j}]=0$ and
$[\hat{b}_{i},\hat{b}_{j}^{\dagger}]=\delta_{ij}$, and $\ket{0}$ is the vacuum
state with no bosons: $\hat{b}_{i}\ket{0}=0$.

The coherent state is an eigenstate of the annihilation operators:
$\hat{b}_{i}\ket{\Theta}=\Theta_{i}\ket{\Theta}$. Therefore, if we only want
the QGN to encode normal ordered expectation values, then Eq. (29) implies
that the images of $Q_{I}^{\dagger}$ only need to contain one state:
$\ket{\Theta}$. We thus obtain a QGN with trivial bond dimensions $\chi_{I}=1$
via the truncation map $Q_{I}=\ket{\widetilde{0}}\bra{\Theta}$. Here,
$\ket{\widetilde{0}}$ labels a state in a Hilbert space of dimension 1. With
this truncation mapping, the local wavefunctions are
$\ket{\psi_{I}}=\ket{\widetilde{0}}$, and the connections are
$V_{IJ}=\ket{\widetilde{0}}\bra{\widetilde{0}}=\mathbbl{1}$. The truncated
[Eq. (11)] annihilation operator at a site $i$ in patch $I$ is simply

$b_{i\in I}=Q_{I}\hat{b}_{i}Q_{I}^{\dagger}=\Theta_{i}$ (72)

This QGN encodes all normal ordered expectation values exactly, e.g.

$\displaystyle\big{\langle}\Theta\big{|}\hat{b}_{i}^{\dagger}\hat{b}_{j}\big{|}\Theta\big{\rangle}$
$\displaystyle=\big{\langle}\psi_{I}\big{|}b_{i\in I}^{\dagger}V_{IJ}b_{j\in
J}\big{|}\psi_{J}\big{\rangle}$ (73) $\displaystyle=\Theta_{i}^{*}\Theta_{j}$
$\displaystyle\big{\langle}\Theta\big{|}\hat{b}_{i}^{\dagger}\hat{b}_{j}^{\dagger}\hat{b}_{k}\hat{b}_{l}\big{|}\Theta\big{\rangle}$
$\displaystyle=\big{\langle}\psi_{I}\big{|}b_{i\in I}^{\dagger}V_{IJ}b_{j\in
J}^{\dagger}V_{JK}b_{k\in K}V_{KL}b_{l\in L}\big{|}\psi_{L}\big{\rangle}$
$\displaystyle=\Theta_{i}^{*}\Theta_{j}^{*}\Theta_{k}\Theta_{l}$

However, expectation values of operators that are not normal ordered are not
encoded correctly by this QGN. For example, $\braket{\psi_{I}}{b_{i\in
I}V_{IJ}b_{j\in J}^{\dagger}}{\psi_{J}}=\Theta_{i}^{*}\Theta_{j}$ while
$\braket{\Theta}{\hat{b}_{i}\hat{b}_{j}^{\dagger}}{\Theta}=\delta_{ij}+\Theta_{i}^{*}\Theta_{j}$.
These additional expectation values could be encoded exactly by adding
additional states to the images $\text{im}(Q_{I}^{\dagger})$, as outlined in
Sec. I.2.3.

###  C.4 Fermion Slater Determinants

We can analytically construct a quantum gauge network that exactly encodes all
normal-ordered two-fermion correlation functions
$\braket{\hat{c}_{i}^{\dagger}\hat{c}_{j}}$ for a fermionic Slater
wavefunction. If there are $n_{\text{f}}$ filled states, we can construct a
QGN with bond dimension $\chi_{I}=1+n_{\text{f}}$. This is more efficient than
the $1+M$ upper bound in Eq. (30), where $M=n$ is the number of operators
whose correlation functions we wish to encode; here we consider all fermion
annihilation operators $\hat{c}_{i}$ with $i=1,\ldots,n$.

A Slater determinant wavefunction can be expressed as

$\ket{\Phi}=\prod_{\alpha=1}^{n_{\text{f}}}\sum_{i=1}^{n}\Phi_{\alpha
i}\,\hat{c}_{i}^{\dagger}\ket{0}$ (74)

using second-quantized Fock states. $i=1,\ldots,n$ indexes the $n$ different
single-particle states. $\hat{c}_{i}$ is a fermion annihilation operator,
which satisfies the anticommutation relations
$\\{\hat{c}_{i},\hat{c}_{j}\\}=0$ and
$\\{\hat{c}_{i},\hat{c}_{j}^{\dagger}\\}=\delta_{ij}$. We fill $K$ many
fermion orbitals, which are index by $\alpha$ and encoded by the matrix
elements $\Phi_{\alpha i}$. The orbitals are assumed to be orthonormalized:
$\Phi\cdot\Phi^{\dagger}=\mathbbl{1}$. The inner product of two Slater
determinant wavefunctions is
$\braket{\Phi}{\Phi^{\prime}}=\det(\Phi^{\prime}\cdot\Phi^{\dagger})$.

The action of an annihilation operator on a Slater determinant wavefunctions
is

$\hat{c}_{i}\ket{\Phi}=\sum_{\alpha=1}^{K}(-1)^{\alpha-1}\Phi_{\alpha
i}\ket{\Phi^{-}_{\alpha}}$ (75)

We define

$\ket{\Phi^{-}_{\alpha}}=\prod_{\alpha^{\prime}\neq\alpha}\sum_{i}\Phi_{\alpha^{\prime}i}\hat{c}_{i}^{\dagger}\ket{0}$
(76)

to be the Slater determinant wavefunction where we do not fill orbital
$\alpha$, but the other $n_{\text{f}}-1$ orbitals are still filled.

Equation (27) implies that the images of $Q_{I}^{\dagger}$ only need to
contain the states $\ket{\Phi}$ and $\ket{\Phi^{-}_{\alpha}}$. Let
$\ket{\phi}$ and $\ket{\phi^{-}_{\alpha}}$ label a basis of $1+n_{\text{f}}$
states for the local QGN Hilbert spaces. Then we can choose truncation maps

$Q_{I}=\ket{\phi}\bra{\Phi}+\sum_{\alpha=1}^{n_{\text{f}}}\ket{\phi^{-}_{\alpha}}\bra{\Phi^{-}_{\alpha}}$
(77)

With this choice, the local wavefunctions are
$\ket{\psi_{I}}=Q_{I}\ket{\Phi}=\ket{\phi}$, and the connections are
$V_{IJ}=Q_{I}Q_{J}^{\dagger}=\mathbbl{1}$. The truncated fermion operators
follow from Eq. (75):

$c_{i\in
I}=Q_{I}\hat{c}_{i}Q_{I}^{\dagger}=\sum_{\alpha=1}^{K}(-1)^{\alpha-1}\Phi_{\alpha
i}\ket{\phi^{-}_{\alpha}}\bra{\phi}$ (78)

This QGN exactly encodes all normal-ordered two-fermion correlation functions:

$\displaystyle\braket{\Phi}{\hat{c}_{i}^{\dagger}\hat{c}_{j}}{\Phi}$
$\displaystyle=\braket{\psi_{I}}{c_{i\in
I}^{\dagger}V_{IK_{1}}V_{K_{1}K_{2}}\cdots V_{K_{l}J}c_{j\in J}}{\psi_{J}}$
(79) $\displaystyle=(\Phi^{\dagger}\cdot\Phi)_{ij}$

where $V_{IK_{1}}V_{K_{1}K_{2}}\cdots V_{K_{l}J}$ is any string of connections
that connect patches $I$ and $J$.

Above, we only worked out analytical expressions for two-fermion correlation
functions. However, analytical expressions for many-fermion correlation
functions should also be possible.

The above QGN is rather trivial in the sense that the connections $V_{IJ}$ are
all identity matrices. This is because we did not take advantage of spatial
locality. If many of the fermion orbitals are spatially local, we expect that
an approximate QGN encoding can be achieved with significantly smaller bond
dimensions and non-identity $V_{IJ}$.

###  C.5 Rainbow State

In Sec. I.2.3, we showed that all $2k$-point correlation functions of $M$ many
operators can be encoded exactly by a QGN with bond dimension $O(M^{k})$. This
is significantly more efficient than a matrix product state (MPS), which can
require bond dimension $\chi^{\text{MPS}}=2^{n/2}$ to encode all two-point
correlation functions of certain states with $n$ qubits, e.g. the rainbow
state. In the $n$-qubit rainbow state, pairs of qubits $i$ and $n+1-i$ are
maximally entangled in a Bell state. The two-point correlation functions of
the rainbow state are
$\braket{\hat{\sigma}_{i}^{\mu}\hat{\sigma}_{j}^{\nu}}=-\delta_{\mu\nu}\delta_{n+1-i,j}$.
$\delta_{\mu\nu}$ denotes the Kronecker delta function. The rainbow state is
the unique state with these correlation functions. Therefore, in order for an
MPS to encode these 2-point correlation functions, the MPS must encode the
rainbow state. Encoding the rainbow state requires MPS bond dimension
$\chi^{\text{MPS}}=2^{n/2}=2^{M/6}$, where $M=3n$ is the number of Pauli
operators $\hat{\sigma}_{i}^{\mu}$ for $n$ qubits.

However, a matrix product operator (MPO) with bond dimension
$\chi^{\text{MPO}}=1+M/2$ is sufficient to encode the 2-point correlation
functions of the rainbow state by encoding the (unphysical) density matrix
$\hat{\rho}=2^{-n}\big{[}\hat{\mathbbl{1}}-\sum_{i=1}^{n/2}\sum_{\mu=\text{x},\text{y},\text{z}}\hat{\sigma}_{i}^{\mu}\hat{\sigma}_{n+1-i}^{\mu}\big{]}$.
This density matrix is unphysical because it has negative eigenvalues. (A QGN
bond dimension of $\chi=1+M/2$ would also be sufficient for this example if we
restrict the allowed operator strings to never change direction.)

##  Appendix D Energy Conservation

Below, we prove that the QGN equations of motion [Eq. (31)] preserve the
energy expectation value [Eq. (33)] exactly when the local Hamiltonian terms
$\hat{H}_{I}$ are time-independent and each supported on a single spatial
patch [as in Eqs. (1) and (32)].

$\displaystyle\partial_{t}\sum_{I}\braket{\psi_{I}}{H_{I}}{\psi_{I}}$
$\displaystyle=$
$\displaystyle\sum_{I}i\braket{\psi_{I}}{[H^{\prime}_{I},H_{I}]}{\psi_{I}}$
$\displaystyle=$ $\displaystyle\sum_{IJ}^{I\cap
J\neq\emptyset}i\braket{\psi_{I}}{[V_{IJ}H_{J}V_{JI},H_{I}]}{\psi_{I}}$ (80)
$\displaystyle=$ $\displaystyle\sum_{IJ}^{I\cap
J\neq\emptyset}\frac{i}{2}\braket{\psi_{I}}{[V_{IJ}H_{J}V_{JI},H_{I}]}{\psi_{I}}$
$\displaystyle\quad\;\;+\frac{i}{2}\braket{\psi_{J}}{[V_{JI}H_{I}V_{IJ},H_{J}]}{\psi_{J}}$
$\displaystyle=$ $\displaystyle\sum_{IJ}^{I\cap
J\neq\emptyset}\frac{i}{2}\braket{\psi_{I}}{[V_{IJ}H_{J}V_{JI},H_{I}]+[H_{I},V_{IJ}H_{J}V_{JI}]}{\psi_{I}}$
$\displaystyle=$ $\displaystyle\,0$

The first three equalities respectively follow from Eq. (31) for
$\partial_{t}\ket{\psi_{I}}$; Eq. (32) for $H^{\prime}_{I}$; and symmetrizing
the sum over $I\leftrightarrow J$. $\sum_{IJ}^{I\cap J\neq\emptyset}$ denotes
the sum over all patches $I$ and $J$ that have nonzero overlap. The final
equality follows from the antisymmetry of the commutator. The second to last
equality follows from:

$\displaystyle\braket{\psi_{J}}{[V_{JI}H_{I}V_{IJ},H_{J}]}{\psi_{J}}$ (81)
$\displaystyle=$
$\displaystyle\braket{\psi_{J}}{V_{JI}H_{I}V_{IJ}H_{J}-H_{J}V_{JI}H_{I}V_{IJ}}{\psi_{J}}$
$\displaystyle=$
$\displaystyle\braket{\psi_{I}}{H_{I}V_{IJ}H_{J}V_{JI}-V_{IJ}H_{J}V_{JI}H_{I}}{\psi_{I}}$
$\displaystyle=$
$\displaystyle\braket{\psi_{I}}{[H_{I},V_{IJ}H_{J}V_{JI}]}{\psi_{I}}$

which follows from $V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ [Eq. (12)].

##  Appendix E Ising Model Quench

In this appendix, we benchmark the quantum gauge network by studying the
dynamics following a quench to a near-critical Ising model. We start form the
initial state $\ket{\Psi(0)}=\otimes_{i}\ket{\rightarrow_{i}}$ where all
$\braket{\sigma_{i}^{\text{x}}}=1$, and then we time evolve with a near-
critical transverse field Ising Hamiltonian

$\hat{H}^{\text{Ising}}=-\sum_{\langle
ij\rangle}\hat{\sigma}_{i}^{\text{z}}\hat{\sigma}_{j}^{\text{z}}-h\sum_{i}\hat{\sigma}_{i}^{\text{x}}$
(82)

with $h=3$ on a two-dimensional $4\times 4$ square lattice with periodic
boundary conditions. (The critical point is at $h_{c}\approx 3.045$ Hamer
(2000).) This system size is chosen so that we can compare to exact methods
that calculate the full wavefunction
$\ket{\Psi(t)}=e^{-i\hat{H}^{\text{Ising}}t}\ket{\Psi(0)}$.

In order to make use of Eqs. (1) and (32), we define the Hamiltonian
$\hat{H}_{I}$ on each spatial patch to be

$\hat{H}_{I=\langle
ij\rangle}^{\text{Ising}}=-\hat{\sigma}_{i}^{\text{z}}\hat{\sigma}_{j}^{\text{z}}-\tfrac{1}{4}h\,(\hat{\sigma}_{i}^{\text{x}}+\hat{\sigma}_{j}^{\text{x}})$
(83)

We take each spatial patch $I$ to be a pair of nearest-neighbor sites $\langle
ij\rangle$. Note that in the sum $\hat{H}=\sum_{I}\hat{H}_{I}$ from Eq. (1),
each site is summed over four times on a square lattice; thus we require the
above $\tfrac{1}{4}$ factor in front of $h$.

We initialize the QGN using truncation maps (as described in Sec. I.1), which
are chosen using a method similar to the one described in Sec. II.1. However
in this spin model, we do not have a conserved charge. Therefore, we modify
step 2 of the method in Sec. I.1 [paragraph below Eq. (36)] to the following:
(2) For each patch $I$, we add states to the image of $Q_{I}^{\dagger}$ that
can be obtained from the current image of $Q_{I}^{\dagger}$ by acting with
Pauli operators within the patch $I$.

The truncation at each patch $I$ only retains states consisting of a span of
eigenstates of the $\hat{\sigma}_{i}^{\text{x}}$ operators. With a natural
gauge choice for the truncation maps, the truncated Pauli operators take the
form of a Kronecker product:

$\sigma_{i\in
I}^{\mu}=Q_{I}\hat{\sigma}_{i}^{\mu}Q_{I}^{\dagger}=\mathbbl{1}\otimes\sigma^{\mu}$
(84)

$\sigma_{i\in I}^{\mu}$ is the truncated [Eq. (11)] Pauli operator at site $i$
for patch $I$, and $\sigma^{\mu}$ is a $2\times 2$ Pauli matrix.

In Fig. 5, we show QGN simulation data for Pauli expectation values
$\braket{\hat{\sigma}^{\mu}_{i}(t)}$ and compare to the exact values. We see
that the simulation errors expectation value decrease as we increase the bond
dimension.

In the QGN, the expectation values are estimated as

$\braket{\hat{\sigma}_{i}^{\mu}}_{\text{QGN}}=\sum_{I\ni
i}^{\text{mean}}\braket{\psi_{I}}{\sigma_{i\in I}^{\mu}}{\psi_{I}}$ (85)

where $\sum_{I\ni i}^{\text{mean}}$ averages over all patches $I$ that contain
the site $i$. In this example, $\braket{\psi_{I}}{\sigma_{i\in
I}^{\mu}}{\psi_{I}}$ is equal for all patches $I$ that contain site $i$ due to
spatial symmetries. However in other models with less symmetry, simulation
errors can make these expectation values differ for different patches.

(a)

(b)

Figure 5: Simulation data for the time dynamics of the near-critical
transverse field Ising Hamiltonian (82) on a periodic $4\times 4$ square
lattice following a quench from the state
$\ket{\Psi(0)}=\otimes_{i}\ket{\rightarrow_{i}}$. (a) The quantum gauge
network (QGN) approximation for $\braket{\sigma_{i}^{\text{x}}}$ vs time $t$
for different bond dimensions $\chi$ (colored lines) vs the exact value (black
line). The legend also shows the number of CPU core hours used for each
simulation. Due to symmetry, the $\braket{\sigma_{i}^{\text{y}}}$ and
$\braket{\sigma_{i}^{\text{z}}}$ expectation values (not shown) are exactly
zero for all time for both the QGN and the exact value. (b) The error
$\braket{\sigma_{i}^{\text{x}}}_{\text{QGN}}-\braket{\sigma_{i}^{\text{x}}}_{\text{exact}}$
of the QGN approximation to $\braket{\sigma_{i}^{\text{x}}}$. The error stays
small for longer times as we increase the bond dimension $\chi$.

If we were to integrate the equations of motion exactly, then the energy
expectation value [Eq. (33)] would be conserved exactly. Since exact
integration is not practical, we use a modified RK4 Runge-Kutta method for
integration with time step $\delta_{\text{t}}=0.02$. Due to this
approximation, the energy per site changed by $5\times 10^{-4}$ and $1\times
10^{-4}$ for the $\chi=88$ and $2028$ simulations, respectively. See Appendix
F for more details.

##  Appendix F Modified Runge-Kutta

We use a modified RK4 Runge-Kutta method to integrate the differential
equations. RK4 is a forth-order Runge-Kutta method that results in an
$O(\delta_{\text{t}}^{4})$ error at time $t\sim 1$, where $\delta_{\text{t}}$
is the time step size. However, the straight-forward application of Runge-
Kutta will not preserve $V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ [Eq. (12)]
exactly; it will only be preserved up to $O(\delta_{\text{t}}^{4})$ error. In
this work, we chose to modify the Runge-Kutta method slightly such that
$V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$ is preserved exactly (i.e. up to
floating point precision). We end up making an additional approximation that
increases the simulation error to $O(\delta_{\text{t}}^{3})$ (which we were
satisfied with). It would be useful to improve the approximation such that
$O(\delta_{\text{t}}^{4})$ and smaller errors can be achieved while
maintaining $V_{IJ}\ket{\psi_{J}}=\ket{\psi_{I}}$.

Instead of integrating $\partial_{t}\ket{\psi_{I}}$ and $\partial_{t}V_{IJ}$
directly at each time step, we use a modified Runge-Kutta method to obtain
estimates for the unitary evolution

$U_{I}^{\text{RK}}(t+\delta_{\text{t}},t)=e^{-i\delta_{\text{t}}H_{I}^{\text{RK}}(t)}$
(86)

We then update the QGN from time $t$ to $t+\delta_{\text{t}}$ as follows:

$\displaystyle\ket{\psi_{I}(t+\delta_{\text{t}})}$
$\displaystyle=U_{I}^{\text{RK}}(t+\delta_{\text{t}},t)\ket{\psi_{I}(t)}$ (87)
$\displaystyle V_{IJ}(t+\delta_{\text{t}})$
$\displaystyle=U_{I}^{\text{RK}}(t+\delta_{\text{t}},t)V_{IJ}(t)U_{J}^{\text{RK}}(t+\delta_{\text{t}},t)^{\dagger}$

We obtain

$H_{I}^{\text{RK}}(t)=\sum_{k=1}^{s}b_{k}\widetilde{H}_{I}^{(k)}(t)$ (88)

using the Runge-Kutta coefficients $b_{k}$, where $s$ is the number of Runge-
Kutta stages. $s=4$ for RK4.

In order to calculate $\widetilde{H}_{I}^{(k)}(t)$, we first recursively
define $H_{I}^{(k)}(t+c_{k}\delta_{\text{t}})$ (without the tilde), which we
set equal to $H^{\prime}_{I}$ [Eq. (32)] at time $t+c_{k}\delta_{\text{t}}$
evaluated using the QGN that is updated from time $t$ to
$t+c_{k}\delta_{\text{t}}$ by the unitary

$\widetilde{U}^{(k)}_{I}(t+c_{k}\delta_{\text{t}},t)=e^{-i\sum_{l=1}^{k-1}a_{kl}\delta_{\text{t}}\widetilde{H}_{I}^{(l)}(t)}$
(89)

via Eq. (87) (with $\widetilde{U}^{(k)}_{I}$ and $t+c_{k}\delta_{\text{t}}$ in
place of $U_{I}^{\text{RK}}$ and $t+\delta_{\text{t}}$). $a_{kl}$ and $c_{k}$
are additional Runge-Kutta coefficients. For $k=1$, $c_{1}=0$ and
$\widetilde{H}_{I}^{(1)}=H_{I}^{(1)}=H^{\prime}_{I}$ at time $t$. We find that
choosing $\widetilde{H}_{I}^{(l)}(t)=H_{I}^{(l)}(t+c_{k}\delta_{\text{t}})$
results in $O(\delta_{\text{t}}^{2})$ simulation errors after time $t\sim 1$
using the RK4 coefficients. We instead use

$\displaystyle\widetilde{H}_{I}^{(l)}(t)=\tfrac{1}{2}\widetilde{U}^{(l)}_{I}(t+c_{l}\delta_{\text{t}},t)^{\dagger}$
$\displaystyle
H_{I}^{(l)}(t+c_{k}\delta_{\text{t}})\widetilde{U}^{(l)}_{I}(t+c_{l}\delta_{\text{t}},t)$
$\displaystyle+\;\tfrac{1}{2}$ $\displaystyle
H_{I}^{(l)}(t+c_{k}\delta_{\text{t}})$ (90)

for which we observe an $O(\delta_{\text{t}}^{3})$ simulation error after time
$t\sim 1$. For RK4, the tableau of coefficients is

$\begin{array}[]{c|cccc}c_{1}&&&&\\\ c_{2}&a_{21}&&&\\\
c_{3}&a_{31}&a_{32}&&\\\ c_{4}&a_{41}&a_{42}&a_{43}&\\\
\hline\cr&b_{1}&b_{2}&b_{3}&b_{4}\end{array}\;\;\scalebox{1.2}{=}\;\;\begin{array}[]{c|cccc}0&&&&\\\
1/2&1/2&&&\\\ 1/2&0&1/2&&\\\ 1&0&0&1&\\\ \hline\cr&1/6&1/3&1/3&1/6\end{array}$
(91)