Quantum Gauge Networks: A New Kind of Tensor Network

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


Introduction
Tensor network algorithms [2,3,4,5,6,7] are very useful for simulating strongly-correlated quantum physics.In one spatial dimension, matrix product state (MPS) algorithms [3,4,5] are often the best available tool for this task.Although still useful, tensor network algorithms in higher dimensions [3,5,8,9,10,11,12,13,14,15,16] 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 [1], which adds "gauge fields" to Schrödinger's picture in order to make spatial locality explicit in the equations of motion.In the gauge picture, one first chooses a collection of possiblyoverlapping 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.Each patch is assigned a local wavefunction |Ψ I ⟩ (with the same Hilbert space dimension as the usual wavefunction), and the Hilbert spaces of neighboring patches are related by unitary connections ÛIJ (which act on the entire Hilbert space), as depicted in Fig. 1.In the simplest setting, the Hamiltonian is written as a sum over terms ĤI that act within a single patch I: The local wavefunctions and connections time-evolve according to where Ĥ⟨I⟩ = is the sum of local Hamiltonian terms supported on patches that overlap with patch I. Typically, we initialize ÛIJ (0) = 1 and |Ψ I (0)⟩ = |Ψ(0)⟩ at time t = 0, where 1 is the identity operator and |Ψ(t)⟩ is the usual wavefunction in the Schrödinger picture.
The expectation value ⟨Ψ| ÂI |Ψ⟩ of a local operator ÂI that only acts within the patch I can be evaluated in the gauge picture as ⟨Ψ I | ÂI |Ψ I ⟩.To calculate an expectation value ⟨Ψ| ÂI BJ |Ψ⟩ for a product of operators acting on different patches, a connection ÛIJ must be inserted in the gauge picture, as in ⟨Ψ I | ÂI ÛIJ BJ |Ψ J ⟩.Such correlation functions can be used to calculate the density matrix from the gauge picture local wavefunctions and connections.For example, for a system of n qubits, the density matrix is where σµ i are Pauli operators, and here we take each patch to consist of just a single qubit (for simplicity) so that I = 1, . . ., n indexes the qubits/patches along some path.The local wavefunction |Ψ I ⟩ is local in the sense that its dynamics are local and connections are required to extract information about operators outside the patch I. Time evolution preserves the following network of relations: ÛIJ |Ψ J ⟩ = |Ψ I ⟩ ÛIJ ÛJK = ÛIK (5) along with Û † IJ = ÛJI and ÛII = 1.

V I J
ψ I  ψ J  Figure 1: An example of a chain of qubits (black dots) and spatial patches (colored ovals) consisting of pairs of neighboring qubits.
In the gauge picture, a local wavefunction |ΨI ⟩ is associated with each patch I, and the Hilbert spaces of neighboring patches are related by unitary connections ÛIJ .A quantum gauge network analogously consists of local wavefunctions |ψI ⟩ in truncated Hilbert spaces that are related by non-unitary connections VIJ .
In this work, we truncate the Hilbert spaces of the local wavefunctions and connections in the gauage picture so that approximate quantum dynamics 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.Furthermore, the connections allow us to utilize different truncated Hilbert spaces for different patches of space.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 [3]), a QGN only involves matrices and vectors (rather than tensors with many indices) regardless of the spatial dimension.As such, it is natural for a QGN algorithm to only require a computation time (e.g.CPU time) that scales as O(χ 3 ), regardless of the number of spatial dimensions, where χ is the dimension of the truncated Hilbert spaces.χ can be viewed as the bond dimension of the QGN.Thus, the natural O(χ 3 ) computation time for a QGN algorithm is the same as for 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 [8, 9, 10], for which the computation time of the TEBD 3 algorithm is O(χ 12 ) [11]. 1 Thus, QGN 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.For a QGN, the number of parameters scales as χ 2 (due to the matrix-valued connections).Therefore, the computation time per parameter scales as 1 More generally, the TEBD 3 computation time is O(D 10 χ 2 ) + O(Dχ 6 ) when the bond dimension χ of the central bonds is different from the bond dimension D of the other bonds.[11] In 2D, the TEBD 2 computation time scales as χ 7 .[8] χ 3 /χ 2 = χ 1.5 for a QGN.This is the same ratio as MPS algorithms, which are very efficient in 1D.For the isometric tensor network TEBD 3 algorithm in 3D, this ratio scales as χ 12 /χ 6 = χ 2 [11], 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, quantum gauge networks depart from this habit.However, a density matrix can be computed from a QGN similar to Eq. (4) for the gauge picture.Unlike a generic PEPS [3] but similar to a matrix product state (MPS) or isometric tensor network [8, 9, 10], 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. 2, 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. (32)], while an MPS of qubits can require an exponentially larger bond dimension 2 M/6 for k = 1.In Sec. 3, 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.

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 consist of just a single lattice site.Another natural choice is to take the patches to have the same support as the Hamiltonian terms.That is, we might choose patches that are pairs of nearest-neighbor sites if the Hamiltonian terms act on nearest-neighbor sites.A QGN then consists of (1) a local wavefunction |ψ I ⟩ for each spatial patch; (2) non-unitary connections V IJ = V † JI 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.
The local wavefunctions and connections are similar to the those within the gauge picture [1], except the Hilbert space is truncated.If the full Hilbert space has dimension N , then |Ψ I ⟩ and ÛIJ in the gauge picture have dimensions N and N × 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 |ψ I ⟩ and connections V IJ , which have truncated dimensions χ I and χ I ×χ J , respectively, where typically χ I ≪ N .We use capital and lower-case Greek letters (e.g.|Ψ I ⟩ vs |ψ 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.Û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 ÂI in the original Hilbert space, we must also define truncated operators, i.e. χ I × χ I matrices A I , that act on the truncated Hilbert space.Throughout this work, ÂI always denotes an operator that acts within a patch I, and similar for BJ , etc.The truncated operators A I are notationally distinguished from the original operators ÂI by the lack of a hat.In Sec.2.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 σ µ I = 1 ⊗ σ µ , where 1 is an identity matrix and σ µ is a 2 × 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 |Ψ⟩, then we would like the QGN to accurately encode local expectation values; i.e. we want ⟨ψ I |A I |ψ I ⟩ ≈ ⟨Ψ| ÂI |Ψ⟩.Similarly, we typically also want expectation values of string operators to also approximately match, e.g.
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, a density matrix can be extracted from a QGN.Thus, a QGN most generally encodes a mixed state rather than a pure state.For example, analogous to Eq. (4) for a system of n qubits, a density matrix can be approximately extracted from a QGN via where σµ i are Pauli operators.Here, we take each patch to consist of just a single qubit (for simplicity), and we use I = 1, . . ., n to index the qubits/patches along a string of nearest-neighbors, e.g. as in Fig. 2.However, due to the approximations induced by the QGN, different paths (e.g.those in Fig. 2) can yield different density matrices.Also similar to the gauge picture, quantum gauge networks exhibit a local gauge symmetry: where Λ I is a unitary matrix.Expectation values must be invariant under this symmetry.
Since local expectation values ⟨Ψ| ÂI |Ψ⟩ ≈ ⟨ψ I |A I |ψ I ⟩ are encoded in the local wavefunctions |ψ 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 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 |ψ J ⟩ = |ψ I ⟩ holds exactly, QGN connected correlation functions obey the usual identity (with a V IJ inserted): where we abbreviate ⟨A I ⟩ = ⟨ψ I |A I |ψ I ⟩ and Therefore, if the QGN connected correlation function [Eq.(10)] is small, we are guaranteed that ⟨ψ I |A I V IJ B J |ψ J ⟩ ≈ ⟨A I ⟩ ⟨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 .original Hilbert space

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 |Ψ⟩, 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 |Ψ⟩ should be chosen to be a purification of the density matrix.The truncation maps are χ I ×N matrices (where N is the dimension of the full Hilbert space) that satisfy: Therefore Q † I is an isometry matrix whose image includes the wavefunction.(A matrix M is isometric if M † M = 1.)Intuitively, each Q I maps a select subspace of states into a truncated Hilbert space, as depicted in Fig. 3.In the next subsection, we will explain how one can obtain useful truncation maps.
With this data, we can construct the following QGN: Note that V II = 1 due to Eq. (11).Local operators ÂI with support on a spatial patch I are also truncated: Note that this equation can be used for both bosonic and fermionic operators.In the limit that the bond dimension χ I → N approaches the full Hilbert space dimension N , this truncation mapping will result in a QGN that encodes correlation functions [e.g.Eq. (6)] exactly.
The operator norm of V IJ is bounded by Thus, the resulting connections V IJ have singular values that are less than or equal to 1.One can verify that from Eq. (9) 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 Note that not all quantum gauge networks can be obtained from the truncation mapping.For example, generating a QGN by sampling random numbers for |ψ 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) [3,4,5].If the MPS has bond dimension χ, then the QGN will have bond dimensions equal to dχ 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 |ψ I ⟩.

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 ÂI BJ ĈK is encoded exactly if Below, we show that this exact encoding can be achieved by choosing truncation maps such that the image of each Q † I is the span of certain strings of operators involved in the expectation values.
In general, we will find that the required bond dimension χ I for each patch I is bounded by χ I ≤ 1 + 2p I [Eq.(23)], 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. (32)] is sufficient to exactly encode all 2k-point correlation functions of M different operators.

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. (15) exactly for a particular choice of local operators ( ÂI , BJ , and ĈK ) and spatial patches I ̸ = J ̸ = K.Below, we show that any choice of truncation maps with the following images is sufficient: span |Ψ 1 ⟩ , . . ., |Ψ m ⟩ denotes the vector space spanned by the vectors |Ψ 1 ⟩ , . . ., |Ψ m ⟩.We set the bond dimensions to be equal to the vector space dimension of the images: . Specifying these images determines the truncation maps Q I up to a unitary gauge transformation On a computer, a Q † I with the desired image can be calculated from the compact singular value decomposition Here, M I is a matrix of column vector, which each encode one of the wavefunctions contained in the span of im(Q † I ).S I is a χ I × χ I diagonal matrix of nonzero singular values, and For each patch, the image must contain |Ψ⟩ so that Eq. (11) is satisfied.Next, we show that Eq. ( 16) is sufficient to exactly encode the expectation value in Eq. (15): )], and ).Note that even if all the local operators commute, the path of the string matters.For example, Eq. ( 16) 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 im(Q I ).Equation ( 16) 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. ( 16) with superset relations "⊇".Alternatively, we could replace im(Q And with this im(Q † J ), we could then replace im( The proof is similar to Eq. (17).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.

Generic Case
More generally, suppose we want to ensure that a QGN exactly encodes the expectation value of Â(1) where I 1 , . . ., I m is a string of neighboring patches with I m ̸ = I m+1 .The exact encoding [Eq.(20)] is guaranteed by any choice of truncation maps with images that contain: where we recursively define In Eq.We can then take the images to be the span of these vector spaces: im(Q
We can easily bound the necessary bond dimensions χ 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. ( 21) increases the dimension of im(Q † I ) by at most 2 [starting from χ I = 1 since we always have |Ψ⟩ ∈ im(Q † Im )].Therefore, this procedure results in bond dimensions of at most To prove that Eq. ( 21) implies Eq. (20), we first recall that We thus obtain and by again using the identities and Eqs. ( 24) and (25) into the first line of Eq. (20) yields If m 0 is a half-integer, then Eq. (20) follows immediately from V ⌊m0⌋,⌈m0⌉ = Q † ⌊m0⌋ Q ⌈m0⌉ and Eqs.(24) and (25).This completes the proof.

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 τ µ i (indexed by µ and position i).That is, suppose we want a QGN to satisfy Patches I and J can be any patches that respectively contain sites i and j. τ µ i∈I = Q I τ µ i Q † I denotes a truncated operator on patch I. Since patches I and J could be far apart, we insert a string of connections to connect patches I and J, where (I, K 1 , K 2 , . . ., K l , J) is an arbitrary string of neighboring patches.To exactly encode the above correlation functions, Eq. ( 21) implies that it is sufficient for the images im(Q † I ) to include the span of the actions of each operator on the wavefunction: If there are M many operators τ µ i , then the above image requires a QGN with bond dimension of at most More generally, suppose we want to exactly encode all 2k-point correlation functions: 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. ( 21) implies that it is sufficient for the images to include the span of the actions of all products of up to k operators: If there are M many operators τ µ i , then there are different operator products included in the span.Therefore, a QGN with bond dimension of at most χ 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 kpoint 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 χ 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 χ 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.

Time Evolution Algorithm
Since the quantum gauge network is closely related to the gauge picture [1] 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 |Ψ I ⟩ → |ψ I ⟩, ÛIJ → V IJ , and Ĥ⟨I⟩ → H ′ I (defined below) in the gauge picture equations of motion [Eq.(2)] to obtain: If the Hamiltonian Ĥ = I ĤI [Eq.(1)] is a sum of local terms ĤI each supported on a spatial patch I, then we define sums over all patches J that have nontrivial overlap with the spatial patch I. H ′ I is analogous to Ĥ⟨I⟩ in the gauge picture [Eq.(3)].Recall that H I is a χ I × χ I matrix in the truncated Hilbert space, while ĤI is N × N and acts on the full Ndimensional Hilbert space.H I can be obtained using the truncation mapping (13), which also holds when ĤI is time-dependent.
For time-independent Hamiltonians, the QGN approximation of the energy expectation value ⟨Ψ(t)| Ĥ|Ψ(t)⟩ is conserved up to numerical integration errors for the time evolution in Eq. (33).See Appendix D for a proof.
For a QGN with bond dimension χ, the computation time of this algorithm is dominated by multiplication of χ × χ matrices.Therefore, the computation time for each time step scales as O(n V χ 3 ), where n V is the number of connections V IJ used by the QGN.The memory cost scales as O(n V χ 2 ) for storing the connections. 2For local Hamiltonians, n V should be proportional to the number of lattice sites.
2 Let χ IJ be the number of non-zero singular values of V IJ .If χ IJ ≪ min(χ I , χ J ), then it is more efficient to decompose χ and χ IJ = χ are equal, and if the H I can be encoded as sparse matrices with O(χ) nonzero entries (which is typical), then the simulation and memory cost can respectively scale as O(n V χ χ 2 ) and O(n V χ χ).This is asymptotically less costly when χ ≪ χ.However for the simulations in this work, this decomposition is not useful since we chose connections for which χ/χ only ranges from about 0.75 (in one spatial dimension) to 0.5 (in three spatial dimensions).
It is also possible to handle Hamiltonian terms that are not supported on a single spatial patch.Consider the Hamiltonian Ĥ = and τ µ J denotes an operator indexed by µ with support on patch J.For example, we could take the patches to consist of a single qubit, and the τ µ J could be Pauli operators.For local Hamiltonians, (34) could be generalized to where ." denotes the Hermitian conjugate of the preceding terms and ensures that H ′ I is Hermitian.When J • • • K involves more than two patches, we must decide on an ordering of the J • • • K for each I. Eq. (37) follows from a similar generalized expression for Ĥ⟨I⟩ in the gauge picture [1] after projecting onto the Hermitian part (via the h.c.) and replacing ÛIJ → V IJ and τ µ J → τ µ J .However, unlike for Eq.(34), the energy will not be conserved exactly for this choice of H ′ I .

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).We initialize the system with a checkerboard pattern of spinless fermions, for which ⟨n i (0)⟩ = 1 − ⟨n j (0)⟩ at time t = 0 for nearest-neighbor sites i and j; as shown in Figs.5a and 5d.We then timeevolve the system using the following Hamiltonian with nearest-neighbor hoppings and nearest-neighbor repulsive interactions: Each patch I is composed of a pair of nearest-neighbor sites, which are summed over by ⟨ij⟩ in the first line.
ĉi is a fermion annihilation operator, and ni = ĉ † i ĉi is the fermion number operator.
We initialize the QGN using truncation maps, as described in Sec.2.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 , as in Sec.2.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 that only contain the initial state.(2) For each patch I of nearest-neighbor sites, we add states to the image of Q † I that can be obtained from the current image by swapping the two sites within the patch.(3) Stop if the bond dimension is sufficiently large.(4) For each patch I, we add states that are included in patches J that overlap with patch I. (5) Go back to the second step.See Fig. 4 for states that get included for the 1D lattice.
With this initialization algorithm, we do not choose the bond dimension χ precisely.For example, Fig. 4 shows that only χ = 2, 4, 8, 10, 16, . . .are allowed for the 1D chain.This restriction has the advantage that χ is chosen such that the QGN remains symmetric under lattice symmetries.Each image consists of a span of certain eigenstates of the number operators ni .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 continue the procedure 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. 3Figure 5 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 timeconsuming simulation consumed 90 CPU hours, which took less than half a day on an 8-core laptop.
We calculate the number expectation value ⟨n i0 (t)⟩ for sites i 0 that initially have zero fermions.In the QGN, these are estimated as follows: mean I∋i averages over all patches I that contain the site i, and )] number operator at site i for patch I.In this example, ⟨ψ I |n i∈I |ψ 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.(35)] would be conserved exactly.Since exact integration is not practical, we use a modified RK4 Runge-Kutta method for numerical integration with time step δ 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. 6 with no interactions (V = 0).Free (i.e.non-interacting V = 0) fermion systems are efficient to simulate exactly [17] and are typically about as challenging for tensor network methods as interacting fermionic systems.By comparing the 1D and 2D data in Figs. 5 and 6, 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).For short time evolutions, we find that the QGN simulation errors decrease as the bond dimension χ is increased.
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 ̸ = 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.

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 scales as χ 3 , where χ 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.Furthermore, the computational time per variational parameter is χ 3 /χ 2 = χ 1.5 , which is the same remarkably efficient ratio as MPS algorithms.
(2) Fermionic models are simple to handle.(Unlike MERA or PEPS, fermionic swap gate [18, 19] are not needed.)(3) The code is simple since the only tensors involved are vectors ψ 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.
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 (12) 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. 4We 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) Imaginary time evolution in the gauge picture does not yield ground state physics, as it does in Schrödinger's picture.[The same is true for Heisenberg's picture where ÂH (t = −iτ ) = e + Ĥτ Â(0)e − Ĥτ , which isn't even Hermitian.]Thus, imaginary time evolution may not appear to be a useful tool for obtaining approximate ground states using a QGN.However, the imaginary time evolution of a Hamiltonian H can equivalently be expressed as the real time evolution of the non-local Hamiltonian Ĥim (t) = −i[ Ĥ, |Ψ(t)⟩ ⟨Ψ(t)|].Although the nonlocality is non-ideal, this kind of time evolution could be implemented in a QGN with all-to-all connectivity of the connections to yield a QGN ground state algorithm.(4) How well can a QGN encode topological states [26,27,28,29,30]?
There are also many opportunities to significantly improve our QGN time evolution algorithm: (5) For the TEBD algorithm [31], 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.
We find it intriguing that the consistency conditions (V IJ |ψ J ⟩ ≈ |ψ I ⟩ and V IJ V JK ≈ V IK ) in Eq. (9) are precisely the equations for a classical lattice gauge theory [49,50] to be in its ground state when coupled to a Higgs field, where V IJ plays the role of the gauge connection and |ψ 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 noncontractible 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 [51,52,53,54,55,56,57]. [

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 ψ i ∈ C N with |ψ i | = 1, i.e. a normalized complex vector with N components.Each edge of the lattice hosts a U(N ) gauge connection U ij , i.e. an N × N unitary matrix.The energy of a classical U(N ) lattice gauge theory coupled to a Higgs field (on a lattice with triangular plaquettes) is where 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).⟨ij⟩ sums over vertices i and j connected by an edge.

Re(ψ
The energy of this classical lattice gauge theory is minimized when U ij U jk = U ik and U ij •ψ j = ψ i , which is analogous to the consistency conditions in Eqs.(5) and (9).However, there is an important difference in each case: (1) In the gauge picture, ÛIJ is an N ×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.

B Matrix Product State Mapping
In this appendix, we show that any matrix product state (MPS) with bond dimension χ can be mapped to a quantum gauge network with bond dimension dχ 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
A matrix product state is an efficient representation of a wavefunction, where the wavefunction amplitudes are given by matrix products.[3,4,5] The MPS wavefunction for a chain of n sites is 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 |s i ⟩ in this case); and we suppress the bond dimension χ MPS 1 = χ MPS n+1 = 1 lines that are traced out.There is a gauge redundancy (with unitary Λ 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 [58,3] centered at a specific site i: L i , C i , and R i are obtained from M i using gauge transformations such that the following identities are obeyed: 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: Â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, . . ., n, while sharing the same isometries (L i and R i ).These shared isometries obey This simultaneous canonical form can be obtained by sweeping across the MPS multiple times using SVD decompositions or from the Vidal gauge [59].The resulting bond dimensions of the

B.2 Quantum Gauge Network from MPS
We construct a quantum gauge network from an MPS using the truncation maps Q i defined in Sec.2.2, which specify the truncated Hilbert space used by the QGN for each patch.We choose Q † i to map the truncated Hilbert space of the MPS orthogonality center C i to the full Hilbert space: Thus, Q i 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. (12) and the MPS canonical form identities [Eq.( 43)], we find that the local wavefunctions are equal to MPS orthogonality centers, while the connections are equal to tensor products of two MPS isometries: We use the simultaneous canonical forms, which obey Eq. (45) and guarantee that The bond dimensions of the quantum gauge network are therefore The singular values of the χ i × χ i+1 matrix V i,i+1 consist of (χ MPS i+1 ) 2 ones, while the rest are zero.Thus, these V i,i+1 are partial isometry matrices [60], which are matrices that obey V V † V = V (or equivalently matrices whose singular values are either zero or one).

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: In order to apply the truncation mapping [Eq.(12)], we first find a purification of the density matrix: The wavefunction |Ψ mix ⟩ hosts an additional auxiliary qubit with index I = 0.In this example, we take the spatial patches to consist of a single qubit.|Ψ mix ⟩ is a purification of ρmix because tracing out the I = 0 qubit yields the density matrix: We follow Sec.2.2 to derive truncation maps Q I by choosing the images im(Q † I ) 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 ρmix are products of σz I σz J operators.We will therefore ensure that expectation values of products of σz I σz J operators are retained by the truncation.It turns out that this is sufficient to exactly encode all correlation functions for this example.Equation (29) thus implies that we only need to include the action of a single σz I in the image of This image has dimension χ I = 2.A natural gauge choice for Q I consistent with the above is: Equation (12) then yields the following QGN: with truncated operators [Eq.(13)] This quantum gauge network exactly encodes all expectation values of the original reduced density matrix.For example,

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: We can then pick the following truncation map: The resulting QGN follows from Eq. (12): Now the truncated Pauli operators are their natural Kronecker products: where 1 0 = |↑ 0 ⟩ ⟨↑ 0 | + |↓ 0 ⟩ ⟨↓ 0 |.These truncated operators obey their usual on-site algebra, e.g.

C.2 Cat State Example
Encoding expectation values that act on many qubits in more than one spatial dimension can be less straight-forward.For example, consider the following cat state The only Pauli string expectation value that can distinguish |Ψ cat ⟩ from the mixed state ρmix in Eq. ( 48) is the highly-nonlocal product of Pauli σx operators on every qubit: Any Pauli string that does not act on all qubits will have an equal expectation value for ρmix and |Ψ cat ⟩.Now suppose that |Ψ cat ⟩ is a wavefunction for a square lattice of qubits.
A quantum gauge network for |Ψ cat ⟩ should reproduce the same nonlocal expectation value: where the sites 1, 2, . . ., n snake across the square lattice, as depicted in Fig. 2a.However, one may want other choices of paths [e.g.Fig. 2b or 2c] 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.2.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 Θ i and takes the following form: This QGN encodes all normal ordered expectation values exactly, e.g.

Θ b †
However, expectation values of operators that are not normal ordered are not encoded correctly by this QGN.For example, ⟨ψ These additional expectation values could be encoded exactly by adding additional states to the images im(Q † I ), as outlined in Sec.2.2.3.

C.4 Fermion Slater Determinants
We can analytically construct a quantum gauge network that exactly encodes all normal-ordered twofermion correlation functions ⟨ĉ † i ĉj ⟩ for a fermionic Slater wavefunction.If there are n f filled states, we can construct a QGN with bond dimension χ I = 1 + n f .This is more efficient than the 1 + M upper bound in Eq. (32), where M = n is the number of operators whose correlation functions we wish to encode; here we consider all fermion annihilation operators ĉi with i = 1, . . ., n.
A Slater determinant wavefunction can be expressed as using second-quantized Fock states.i = 1, . . ., n indexes the n different single-particle states.ĉi is a fermion annihilation operator, which satisfies the anticommutation relations {ĉ i , ĉj } = 0 and {ĉ i , ĉ † j } = δ ij .We fill K many fermion orbitals, which are index by α and encoded by the matrix elements Φ αi .The orbitals are assumed to be orthonormalized: The action of an annihilation operator on a Slater determinant wavefunctions is We define With this choice, the local wavefunctions are |ψ I ⟩ = Q I |Φ⟩ = |ϕ⟩, and the connections are The truncated fermion operators follow from Eq. (66): This QGN exactly encodes all normal-ordered twofermion correlation functions: 0 since states with two fermions annihilated from |Ψ⟩ are not included in the truncation.However, analytical expressions for a QNG that encodes 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.2.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 χ 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 ⟨σ µ i σν j ⟩ = −δ µν δ n+1−i,j .δ µν 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 χ MPS = 2 n/2 = 2 M/6 , where M = 3n is the number of Pauli operators σµ i for n qubits.
However, a matrix product operator (MPO) with bond dimension χ MPO = 1 + M/2 is sufficient to encode the 2-point correlation functions of the rainbow state by encoding the (unphysical) density matrix ρ = 2 −n 1 − n/2 i=1 µ=x,y,z σµ i σµ n+1−i .This density matrix is unphysical because it has negative eigenvalues.(A QGN bond dimension of χ = 1+M/2 would also be sufficient for this example if we restrict the allowed operator strings to never change direction.)

D Energy Conservation
Below, we prove that the QGN equations of motion [Eq.(33)] preserve the energy expectation value [Eq.(35)] exactly when the local Hamiltonian terms ĤI are time-independent and each supported on a single spatial patch [as in Eqs.(1) and (34)].
The first three equalities respectively follow from Eq. (33) for ∂ t |ψ I ⟩; Eq. ( 34) for H ′ I ; and symmetrizing the sum over I ↔ J.

I∩J̸ =∅ IJ
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:

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 |Ψ(0)⟩ = ⊗ i |→ i ⟩ where all ⟨σ x i ⟩ = 1, and then we time evolve with a near-critical transverse field Ising Hamiltonian with h = 3 on a two-dimensional 4 × 4 square lattice with periodic boundary conditions.(The critical point is at h c ≈ 3.045 [61].)This system size is chosen so that we can compare to exact methods that calculate the full wavefunction |Ψ(t)⟩ = e −i ĤIsing t |Ψ(0)⟩.
In order to make use of Eqs.(1) and (34), we define the Hamiltonian ĤI on each spatial patch to be ĤIsing I=⟨ij⟩ = −σ z i σz j − 1 4 h (σ x i + σx j ) (74) We take each spatial patch I to be a pair of nearestneighbor sites ⟨ij⟩.Note that in the sum Ĥ = I ĤI from Eq. (1), each site is summed over four times on a square lattice; thus we require the above 1  4 factor in front of h.
We initialize the QGN using truncation maps (as described in Sec.2.1), which are chosen using a method similar to the one described in Sec.3.1.However in this spin model, we do not have a conserved charge.Therefore, we modify step 2 of the method in Sec.2.1 [paragraph below Eq. (38)] to the following: (2) For each patch I, we add states to the image of Q † I that can be obtained from the current image of Q † I 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 σx i operators.With a natural gauge choice for the truncation maps, the truncated Pauli operators take the form of a Kronecker product: σ µ i∈I is the truncated [Eq.(13)] Pauli operator at site i for patch I, and σ µ is a 2 × 2 Pauli matrix.
In Fig. 7, we show QGN simulation data for Pauli expectation values ⟨σ µ 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 averages over all patches I that contain the site i.In this example, ⟨ψ I |σ µ i∈I |ψ 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.(35)] would be conserved exactly.Since exact integration is not practical, we use a modified RK4 Runge-Kutta method for integration with time step δ t = 0.02.Due to this approximation, the energy per site changed by 5 × 10 −4 and 1 × 10 −4 for the χ = 88 and 2028 simulations, respectively.See Appendix F for more details.

F Modified Runge-Kutta
We use a modified RK4 Runge-Kutta method to integrate the differential equations.RK4 is a forthorder Runge-Kutta method that results in an O(δ  error at time t ∼ 1, where δ t is the time step size.However, the straight-forward application of Runge-Kutta will not preserve V IJ |ψ J ⟩ = |ψ I ⟩ [Eq.( 14)] exactly; it will only be preserved up to O(δ 4 t ) error.In this work, we chose to modify the Runge-Kutta method slightly such that V IJ |ψ J ⟩ = |ψ 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(δ 3  t ) (which we were satisfied with).It would be useful to improve the approximation such that O(δ 4 t ) and smaller errors can be achieved while maintaining V IJ |ψ J ⟩ = |ψ I ⟩.
Instead of integrating ∂ t |ψ I ⟩ and ∂ t V IJ directly at each time step, we use a modified Runge-Kutta method to obtain estimates for the unitary evolution We obtain using the Runge-Kutta coefficients b k , where s is the number of Runge-Kutta stages.s = 4 for RK4.
In order to calculate G where t k = t + c k δ t and a kl δt G (l) a kl and c k are additional Runge-Kutta coefficients.
For k = 1, c 1 = 0 so that t = t, and V for which we observe an O(δ 3 t ) simulation error after 5 Equivalently, the right-hand-side of Eq. ( 80) is H ′ I from Eq. ( 34) at time t k evaluated using the QGN that is updated from time t to t k by U

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

Figure 3 :
Figure 3: Each truncation map QI maps a subspace of the original Hilbert space on to the truncated Hilbert space associated with patch I.
The truncation maps can be chosen to be Q I = |ψ I ⟩ ⟨Ψ|, with I = 1, . . ., n.This results in a QGN with local wavefunctions |ψ I ⟩ and connections V IJ = |ψ I ⟩ ⊗ ⟨ψ J |.The truncated operators simply act within the on-site Hilbert space.Additional examples of quantum gauge networks can be found in Appendix C.
(21), we are free to choose any half-integer m 0 = 1,3  2 , . . ., 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 im(Q † I ) ⊇ S (m) I from Eq. (21) with m = 1, 2, . . ., M .

4 χ = 2 Figure 4 :
Figure4: States included in the truncated Hilbert space for a particular patch (labelled in red) of the 1D lattice.These are the states generated from five iterations of the algorithm described in the paragraph below Eq.(38).The initial state is shown in the top row, where white and black squares denote empty and filled fermions along the 1D chain (where we only show a subset of the chain near the patch).Below the initial state are additional states generated by the algorithm, where a red dot denotes a site with a different fermion number than the initial state.Each iteration ends on step (3), where we check if the number of states is sufficiently large.Each iteration is separated by a green line in the figure.

Figure 5 :
Figure 5: Simulation data for the time dynamics of the fermionic Hamiltonian(38) 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 ⟨ni 0 (t)⟩ vs time t at a site i0 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 χ (colored lines).The legend also shows the number of CPU core hours used for each simulation.(c) The error of the QGN approximation ⟨ni 0 ⟩ QGN [Eq.(39)] to the exact expression ⟨ni 0 ⟩.Panels (d-f) are similar, but for a periodic 4 × 4 square lattice.In both dimensions, we see that increasing the bond dimension increases the accuracy of the QGN simulations (for sufficiently small times).

Figure 6 :
Figure 6: Panels (a-f) are the same as Fig. 5, except we simulate non-interacting fermions with V = 0.By comparing to Fig. 5, 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 × 4 × 4 cubic lattice.

bi
is a boson annihilation operator, which satisfies the commutation relations [ bi , bj ] = 0 and [ bi , b † j ] = δ ij , and |0⟩ is the vacuum state with no bosons: bi |0⟩ = 0.The coherent state is an eigenstate of the annihilation operators: bi |Θ⟩ = Θ i |Θ⟩.Therefore, if we only want the QGN to encode normal ordered expectation values, then Eq. (31) implies that the images of Q † I only need to contain one state: |Θ⟩.We thus obtain a QGN with trivial bond dimensions χ I = 1 via the truncation map Q I = | 0⟩ ⟨Θ|.Here, | 0⟩ labels a state in a Hilbert space of dimension 1.With this truncation mapping, the local wavefunctions are |ψ I ⟩ = | 0⟩, and the connections are V IJ = | 0⟩ ⟨ 0| = 1.The truncated [Eq.(13)] annihilation operator at a site i in patch I is simply

Figure 7 :
Figure 7: Simulation data for the time dynamics of the near-critical transverse field Ising Hamiltonian (73) on a periodic 4×4 square lattice following a quench from the state |Ψ(0)⟩ = ⊗i |→i⟩.(a) The quantum gauge network (QGN) approximation for ⟨σ xi ⟩ vs time t for different bond dimensions χ (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 ⟨σ y i ⟩ and ⟨σ z i ⟩ expectation values (not shown) are exactly zero for all time for both the QGN and the exact value.(b) The error ⟨σ x i ⟩ QGN − ⟨σ x i ⟩ exact of the QGN approximation to ⟨σ x i ⟩.The error stays small for longer times as we increase the bond dimension χ.

I 5 G
(t k ) (without the tilde) as

( 1 )
IJ (t 1 ) = V IJ (t)and G I from Eq. (34) at time t.We find that choosing G

I
(t l ) results in O(δ 2 t )simulation errors after time t ∼ 1 using the RK4 coefficients.We instead use G

I
(t k , t).

time t ∼ 1 . 1 c 2 a 21 c 3 a 31 a 32 c 4 a 41 a 42 a 43 b 1 b 2 b 3 b 4
For RK4, the tableau of coefficients is c where Λ 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.