Variational Tensor Network Operator

We propose a simple and generic construction of the variational tensor network operators to study the quantum spin systems by the synergy of ideas from the imaginary-time evolution and variational optimization of trial wave functions. By applying these operators to simple initial states, accurate variational ground state wave functions with extremely few parameters can be obtained. Furthermore, the framework can be applied to study spontaneously symmetry breaking, symmetry protected topological, and intrinsic topologically ordered phases, and we show that symmetries of the local tensors associated with these phases can emerge directly after the optimization without any gauge fixing. This provides a universal way to identify quantum phase transitions without prior knowledge of the system.


I. INTRODUCTION
Tensor networks (TN) are preeminent theoretical and computational approaches to study quantum many-body systems [1]. Among them, matrix product states (MPS) and projected entangled pair states (PEPS) receive particular attention to describe the ground states of many-body systems in the thermodynamic limit in one and higher dimensions [2,3]. The success of MPS and PEPS lies in the fact that they efficiently parameterize the area-law states and encapsulate all the information in a single tensor. The former guarantees that they are successful ansatz to represent the ground states of gapped local Hamiltonians, and the latter reduces the problem of classifying phases of matter into studying the symmetries of the local tensors. However, in practical simulations, significant efforts are still devoted to address the following two questions: Given the translationally invariant Hamiltonian, how to optimize the MPS and PEPS ansatz? And after the optimization, how to classify different phases in terms of the local tensors?
One natural way to search for the optimal ground states is to variationally optimize all the parameters of the local tensor with proper gauge-fixing conditions. This includes the well-known density matrix renormalization group (DMRG) and the variational uniform matrix product state (VUMPS) for MPS [4][5][6][7], as well as the gradient-based optimization for PEPS, where the gradients of the parameters are computed using either the sum of several tensor diagrams or the differentiable programming techniques [8][9][10]. While variational optimization (VarOpt) makes full use of the TN ansatz and yields the most reliable ground states, it is computationally demanding, and the optimization process may potentially be trapped to local minima due to a large number of free parameters. This downside is particularly severe for the 2D PEPS and the higher-dimensional tensor networks. Another method to search for the optimal tensor is to perform the imaginary time evolution (ITE). The basic idea is that the ground state can be obtained by evolving a random initial state with the operator e −τĤ for a sufficiently long time. While the resulting state can be represented as a TN with infinite bond dimensions using the Trotter decomposition, the area law guarantees that one can reasonably truncate the bond dimensions. This method contains the infinite time-evolving block decimation (iTEBD) in 1D [11][12][13] and the simple and full updates method in 2D [14][15][16][17]. ITE method is computationally cheap as it contains no variational parameters. However, the Trotter error due to the discretization in time and the necessary truncation leads to loss of accuracy. Therefore, the TN wave functions obtained from ITE are usually biased by the choice of initial state and may not capture the long-range correlation and entanglement accurately.
After the optimization, a recurring theme is to determine the phase of the system from the optimized tensor. While the classification of several quantum phases in terms of MPS and PEPS are well developed [18][19][20][21][22], implementing those theories in practice requires prior knowledge of the system. Furthermore, the techniques developed for a particular phase of matter cannot be easily generalized to others. This comes from the fact that the MPS and PEPS representations are not unique, and the characteristics of their building block tensors in different phases become clearer in different gauge choices. For example, since the one-dimensional symmetry-protected topologically ordered (SPT) phases are characterized by the projective representations of the symmetry group, the specific numerical approach was developed to extract the corresponding representations from the MPS [23]. On the other hand, the intrinsic topologically ordered phases (TO) are characterized by the global symmetries on the virtual bonds, and one should impose the virtual symmetry during the optimization [24,25]. We note that some unbiased techniques to study the topologically ordered phases without imposing symmetries during the optimization were proposed recently [26,27]. However, one should perform tedious gauge fixing procedures after the optimization, and this still requires the knowledge of what kind of TO phase we are studying.
In this paper, we address the aforementioned questions by proposing a simple and generic construction of the variational tensor network operator (TNO) which can be derived merely from the knowledge of the model Hamiltonian. We refer to this operator as Generic-TNO or GTNO in the following. The generic tensor network state (GTNS) containing variational parameters is then obtained by applying GTNO to some simple initial states. GTNO combines the idea of ITE and VarOpt by constructing a variational operator which projects the simple initial state to the ground state. Specifically, consider the translationally invariant Hamiltonian of the formĤ = h (l) i1,i2,...,i l , where h (l) i1,i2,...,i l is the local interaction with an l-site support. The method of ITE, i.e., applying e −τĤ for a sufficiently long time, implies that the ground state can be obtained by applying some linear combinations of the operators generated from the Hamiltonian. While in general one cannot construct the exact ground-state projector efficiently fromĤ, one can use the local interaction h (l) i1,i2,...,i l to build a variational operator satisfying all the symmetries ofĤ. The key insight of GTNO is that all the powers of h (l) i1,i2,...,i l can be decomposed into the smaller tensors. It is then natural to construct a variational TNO using those tensors as the building blocks. On the side of obtaining ground states, GTNS needs far fewer parameters but yields comparable energy with VarOpt. On the side of classifying phases of matter, symmetries of the local GTNS emerge themselves after the optimization. Different from most of the other numerical approaches that actively look for the symmetries of the local tensor, GTNS let those symmetries come looking for it after the optimization, which is extraordinary. As we will demonstrate below, without strictly imposing any symmetry before or performing gauge fixing after the optimization, GTNS can be used to distinguish several quantum phases, including the spontaneously symmetry breaking (SSB), SPT, and the intrinsic TO phases. Our study provides a unified protocol to study the general quantum spin systems and paves the way to explore new exotic phases.
The paper is organized as follows. In Sec.II, we motivate the construction of GTNO by noting that one can exponetiate the local interaction h (l) i1,i2,...,i l to build an efficient operator to project the initial state to the ground state. We then show that any power of the two-site interaction can be expressed as the contraction of two rank-3 tensors, which form the building blocks of the GTNO. In Sec.III, we explicitly construct the 1D GTNO, the generic matrix product operator (GMPO) from those building block tensors and discuss the symmetry properties as well as the physical meaning. After establishing the GMPO form, we demonstrate its power by studying the 1D Heisenberg, transverse field Ising (TFI), and transverse field cluster (TFC) models. In Sec.IV, we generalize the GTNO to the 2D systems, the generic projected entangled product operator (GPEPO). The GPEPO is then applied to study the 2D Heisenberg, TFI, and the toric code model in a magnetic field. There, we also compare the results obtained by the GPEPO framework with other numerical approaches. Finally, we conclude in Sec.V and discuss several possible extensions of the present work.
The exact ground state projector lim τ →∞ e −τĤ can then be decomposed into the tensor product of all the local projectors , which can be further expanded as Therefore, the exact projector to the ground state subspace can be written as the sum of tensor products of h (l) i1,i2,...,i l with proper coefficients. On the other hand, if the Hamiltonian is not frustration-free, the exact ground state projector should be written as the power series expansion of the global Hamilto-nianĤ: ..,i l ] = 0 in general,Ĥ m produces some operators that cannot be written as a tensor product of local terms, inhibiting an efficient TN representation. However, the local nature of h (l) i1,i2,...,i l implies that most of the terms expanded in Eq. (2) do have a tensor-product structure. Therefore, it is still possible to establish an efficient TNO using merely the power of h (l) i1,i2,...,i l rather thanĤ. As we will show in Sec.III and IV, the non-commuting property of the local interaction can be remedied using a symmetrizing procedure in the GTNO framework.
The key observation of our work is that all the powers of h (l) i1,i2,...,i l can be decomposed into the smaller tensors. The variational operator can then be naturally written as a TNO composed of those tensors. In the following, we demonstrate how to obtain those building block tensors by considering the l = 2 case, i.e., the nearest-neighbor interactions, and show that they naturally inherit the on-site symmetries from the Hamiltonian.
For any two-site interaction h (2) i,j , one can express the αth power of it as a sum of the tensor product of two local operators: Fig. 1(a)]. In the second equality, we regard the collection of (A α 0 , A α 1 , · · · , A α Dα−1 ) and (B α 0 , B α 1 , · · · , B α Dα−1 ) as the vector operators (A α | and |B α ), i.e., the vector in the virtual space v with its element as the physical operator. Here, we use |·) to denote a vector in the virtual Hilbert space as opposed to the physical one |· . It then follows that (A α |B α ) denotes the inner product in the virtual space.
If the Hamiltonian respects the global on-site symmetry G: U (g)ĤÛ † (g) =Ĥ, ∀g in G withÛ (g) = i u i (g), where u i (g) is the local representation of g, then (h (2) i,j ) α is invariant under the local on-site transformation: [see Fig. 1(b)]. This implies that where V α (g) is a unitary representation for the symmetry group G [ Fig. 1(c)]. In other words, under the physical transformation, (A α | and |B α ) transform like the bra and ket vectors in the virtual Hilbert space, respectively. Interestingly, for several physical models, (h (2) i,j ) α , α = 0, 1, ..., ∞ are not linearly independent of one another. In other words, one can find an integer n such that for any m ≥ n, (A m |B m ) is not linearly independent of {(A α |B α )} for α = 0, 1, ..., n − 1. This closed condition indicates that any polynomial function of h (2) i,j , i.e., P h m , can be written as linear combination of (A α |B α ) for α = 0, 1, ..., n − 1 with coefficients c α : For example, one can show that all the Hamiltonians composed of the Pauli matrices satisfy Eq. (6) using the relation σ i σ j = δ ij I + i ijk σ k . The benefit of the closed condition is that Eq. (6) implies that the local projector in Eq. (1) can be written as the linear combination of (A α |B α ), α = 0, ..., n − 1. Therefore, if the Hamiltonian is frustration-free, one can always obtain the exact ground state projector by variationally optimizing n parameters c α . In Sec.III and Sec.IV, we will show that GTNO can always represent the exact projector to the ground-state subspace of a frustration-free Hamiltonian. On top of that, GTNO allows accurate description for the non frustration-free systems and can be used to distinguish different quantum phases. Before explicitly constructing the GTNO using (A α | and |B α ) in the next section, we provide two concrete examples of the above decompositions.
Case I: Ising model. The local interaction of the Ising model is of the form which is invariant under the global Z 2 transformation: Given the knowledge of |A α ) and |B α ), we can identify the virtual unitary representation V α as V 0 = 1 and Case II: Heisenberg model. The local interaction of spinhalf Heisenberg model is of the form which is invariant under the global SU (2) transformation: Since σ i σ j = δ ij I + i ijk σ k , we find that n = 2 with By direct computation of |u γ (θ)A 0 u γ (θ) † ) and |u γ (θ)A 1 u γ (θ) † ), we identify V 0 = 1 and V 1 = e −iJ γ 3 θ , where J γ 3 is the generator of the three-dimensional irreducible SU (2) representation.

A. General Construction
We now demonstrate how to construct the 1D GTNO, i.e., the GMPO, given the knowledge of (A α | and |B α ), α = 0, ..., n − 1. Let (A| = ⊕ n−1 α=0 (A α | and |B) = ⊕ n−1 α=0 |B α ), the GMPO with bond dimension D = n−1 α=0 D α can be constructed as The constructions of GMPO given the pairs of (A α | and |B α ). Note that the outgoing and incoming bonds are always attributed to the three-leg tensors (A| and |B), respectively. (b) The representation power of GMPO can be improved by applying it several times and assigning independent parameters for all the blocks of where A v1 B v2 denotes the fully symmetrized product. This can be schematically represented in Fig. 2(a). Here we use the symbol "∼" to emphasize the fact that GMPO is not exactly equal to A v1 B v2 /2! but contains n 2 variational parameters for each block of A α v1 B α v2 /2 with α, α = 0, ..., n − 1 in general. Note that the use of the fully symmetrized product is to take the non-commuting property of different local interactions h (2) i,j into account, which acts trivially if the Hamiltonian is frustration-free. In other words, if [h (2) i,j , h This implies that with proper choice of the free parameters, GMPO can represent the projector to the ground state subspace of the frustration-free Hamiltonian.
Besides, since (A α | and |B α ) transform like the bra and ket vectors in the virtual Hilbert space under the on-site physical transformation [see Eq. (5)], it is straightforward to show that s1,s2 where V = ⊕ n−1 α=0 V α . In terms of the TN representation, .
The implication of Eq. (13) can be understood by considering the global physical operatorĜ formed by contracting all the virtual indices of G: G = s ,s · · · G s j sj G s j+1 sj+1 · · · |s s|, where s = (· · · , s j , s j+1 , · · · ). Since the virtual matrices V (g) cancel out when contracting all the virtual bonds, one can easily see that the global GMPOĜ commutes with the global on-site symmetry operation [Ĝ,Û ] = 0. In other words, GMPO inherits all the on-site symmetries of the Hamiltonian.
Aside from the on-site symmetries, we would like to require GTNO to respect all the symmetries of the Hamiltonian, as the spirit of GTNO is to construct a variational operator generated from the Hamiltonian. This also allows us to further reduced the number of parameters. For example, from Eq. (8), we find that (A| = (B| = I σ z for the Ising model. Here we use the vertical line to distinguish operators arising from different α = 0, ..., n − 1. Enforcing the inversion symmetry and the normalization condition, the GMPO with two variational parameters c 1 , c 2 can be written as or equivalently, G 00 = I, To get more intuitive understanding of Eq. (14), we represent the corresponding global operatorĜ schematically in Fig. 3(b) by denoting the virtual index 0(1) as the black(red) string. One can observe that different types of string correspond to different physical actions generated from the power of h i ,j ] = 0 if i = i or j = j in general, GMPO symmetrizes the physical action locally through A v1 B v2 /2! when different strings are connected. The resulting global operatorĜ then generates several physical operators, both local and non-local, arising from the power of global HamiltonianĤ m for some m. The parameters then correspond to the weights of the different connections of strings.
Nonetheless, the current GMPO constructed using Eq. (11), which we term the simple GMPO, cannot generate all the terms existing in the power ofĤ. This is manifested by considering the case of the Heisenberg model. From Eq. (10), we get |A) = |B) = I σ x σ y σ z . The GMPO with bond dimension D = 4 and two parameters can then be written as (imposing the inversion symmetry and the normalization condition) where we allow different signs for the G and G T to capture both the ferromagnetic and anti-ferromagnetic interaction. The vanishing of the off-diagonal components in the second diagonal block upon the fully symmetrized product is due to the anti-commutation relation for the Pauli matrices. This reflects the fact that when considering the power ofĤ, h j,k always appear equal number of times as h i,j , and thus the terms like σ x i (σ x j σ y j )σ y k will never exist. However, the simple GMPO cannot generate some long-range physical actions like σ x i (σ x j σ y j )(σ y k σ z k )σ z l which do exist from the power ofĤ. To include those terms in our variational ansatz, one can directly apply the simple GMPO several times (say, s times) as shown in Fig. 2(b). The resulting GMPO has bond dimension D = n−1 α=0 D α s with n 2s variational parameters.
Schematically, this will introduce more string configurations generated from the power ofĤ. Another way to enhance the representation power of GMPO is to increase the number of times each A α and B α appears (say, t α times): The resulting GMPO constructed using Eq.
Physically, this introduces more independent weights for different string configurations. For example, by introducing |A 1 ) and |B 1 ) two times for the Ising model, the GMPO becomes Regarding both the virtual index 1 and 2 as the red string (since they represent the same physical action), it is then obvious that the additional parameters allow us to assign more independent coefficients for the string configurations in Fig. 3(b). In the following, we choose to apply the simple GMPO several times and assign independent parameters to each GMPO for simplicity. In other words, for the s-th order variational GMPS |ψ (s) =Ĝ s |ψ 0 with some initial state |ψ 0 , the total GMPO has bond dimension D = n−1 α=0 D α s with s × n 2 free parameters. We will show that even with such a simple construction, it is already adequate in reflecting the ground-state properties. Note that ifÛ |ψ 0 = |ψ 0 , the resulting GMPS |ψ (s) will also respect the on-site symmetry since [Ĝ,Û ] = 0. This turns out to be a very efficient way to construct the symmetric TN wave functions.
Since the spirit of GTNO lies in superpoisng several operators generating from the Hamiltonian, what really matters is the fact that the building-block tensors (A α | and |B α ) satisfy Eq. (5) and Eq. (6). Hence, we generalize the meaning of the superscript α and it no longer needs to represent the power of h (2) i,j . For example, given the XXZ model h which has exactly the same form as the Heisenberg model but with more variational parameters. This is reasonable since XXZ model has less symmetries than the Heisenberg model, and hence its GMPO is less restrictive.

B. Adding One-Site Terms
Now, we consider the more general Hamiltonian by adding the one-site terms: as the special case of h (2) i,j , here we treat them as different objects, which allows us to assign more variational parameters to the ansatz. The generalization is straightforward: consider the β-th power of h (1) i , which corresponds to a two-leg tensor, i.e., a physical one-site operator, C β . If the Hamiltonian is invariant under the on-site symmetry transformation: In other words, under the physical transformation, C β transforms like a scalar in the virtual Hilbert space. If h (1) i is made out of the Pauli matrices, any polynomial function of h For example, the one-site transverse magnetic field h GMPO can then be constructed as In terms of the TN representation, the left hand side of the Eq. (21) corresponds to The n 2 ×n o variational parameters can then be assigned for each block of A α B α C β , where α, α = 0, ..., n − 1 and β = 0, ..., n o − 1. Note that since C β transforms trivially under the physical transformation, the virtual gauge transformation matrix V (g) is fully determined by |A) and |B). Another useful point in the practical simulations is that one can easily write down a program to construct the GMPO using Eq. (21) without even knowing the explicit form of it.
From Sec.III C to Sec.III E, we will use GMPO to study the 1D Heisenberg, transverse field Ising (TFI), and transverse field cluster (TFC) models. Specifically, we construct the symmetric GMPS by applying the GMPO to some simple initial states respecting the symmetries of the Hamiltonian, and then variationally optimize the parameters of GMPS to search for the ground states. Here we emphasize that GMPS is not merely a class of symmetric TNS [28]. The idea of the latter is to sort out all the possible classes of the TNS given the symmetry requirements. After the classification, one should variationally optimize all the classes to nail down the phases of the system. In contrast, the main objective of the GMPS is not the symmetries but the Hamiltonian of the given system. Constructing the GMPO directly from the Hamiltonian and then applying it to some symmetric initial state allows us to obtain the symmetric ansatz effortlessly, bypassing the tedious classification of the symmetric TNS. Furthermore, as we will demonstrate below, there exist extra symmetries of the local tensor that only emerge for some parameter ranges of the GMPS. These emergent symmetries serve as the probes of several quantum phase transitions, including the spontaneously symmetry breaking and the symmetry protected topological phase transitions.

C. Case I: Heisenberg Model
We now use the GMPO of the form in Eq. (15) to study the S = 1/2 Heisenberg model. Since the system is SU (2)-invariant and critical, the ground state cannot be represented exactly as an injective MPS with a finite bond dimension [29]. We choose the initial state as a Majumdar-Gosh state [30] respecting both the translational and SU (2) symmetries, which can be written as a non-injective MPS |ψ 0 = s (· · · M sn−1 M sn M sn+1 )|s with local tensor M = Applying the GMPO several times then generates the SU (2) symmetric GMPS |ψ (s) =Ĝ s |ψ 0 with 2s free parameters. To find the optimal parameters for the lowest-energy state, we employ the techniques of the differentiable programming tensor network [31,32], which uses the automatic differentiation to compute the gradient of all the parameters with respect to the energy. This allows us to update all the parameters each time we evaluate the energy. Our code implementation for the 1D calculations is publicly available at Ref. [33] Tab. I shows the energy deviation δE = (E calc − E exact )/|E exact | comparing with the exact solution [34] and the correlation length for s = 1, 2, 3. While the Hamiltonian is not frustration-free, GMPO can get reliable SU (2)-symmetric ground state with δE = 1.31 ≈ 10 −5 and the extremely large correlation length ξ ≈ 205 unit cells by optimizing merely six free parameters. This demonstrates that GMPO accurately captures the long-range correlation of the system.

D. Case II: Transverse Field Ising Model
The Hamiltonian of the TFI model is written aŝ which enjoys the global Z 2 symmetryÛ = i σ x i . This model is exactly solvable by using the Jordan-Wigner transformation and exhibits the phase transition between the polarized and spontaneously symmetry-broken (SSB) phases at h x = 1. Using Eq. (21), the GMPO with bond dimension D = 2 composed of 4 variational parameters can be written as To gain some intuition, we first consider the extreme limits h x = 0, ∞ such that the Hamiltonian is frustration free. When h x → ∞, the ground state is polarized in the x direction and can be obtained by applying the projector to any random state |ψ 0 having non-zero overlap with the ground state. Up to a normalization constant, this corresponds to the GMPO with parameters c 3 = 1 and c 1 = c 2 = c 4 = 0.
On the other hand, when h x = 0, the ground state is two-fold degenerate and its subspace can be obtained by applying the projector This corresponds to the GMPO with parameters c 1 = c 2 = 1 and c 3 = c 4 = 0. Therefore, GMPO reproduces the projectors to the ground state subspace in the two extreme limits. Now, we choose the initial state as the product state polarized in the x-direction: |ψ 0 = ⊗ i |+ i , which is the only product state respecting the Z 2 spin-flip symmetry and having non-zero overlap with the true ground state for all h x ≥ 0. The resulting s-th order GMPS |ψ (s) =Ĝ s |ψ 0 then also respects the Z 2 symmetry. For the first-order GMPS |ψ (1) , since σ x |+ = |+ , c 3 and c 4 can be discarded. Fig. 4(a) shows the evolution of c 1 and c 2 with respect to the transverse field for |ψ (1) . One can clearly observe that both of them exhibit a sharp change at h x ≈ 0.828, indicating the occurrence of the phase transition. To understand the physical meaning of the significant difference between the local tensors above and below h x ≈ 0.828, we first consider the limit h x = 0. The GMPO yields the Greenberger-Horne-Zeilinger (GHZ) state |ψ (1) = 1 i=0 |i, i, · · · , i when c 1 = c 2 = 1, and the local tensor is of the form which is non-injective and exhibits a global virtual Z 2 symmetry [see Fig. 4(b)], which detects whether the local tensor respects the extra virtual Z 2 symmetry. Here |·| denotes the Frobenius norm of the tensor. Fig. 4(c) shows the VOP as the function of the transverse field. We found that the VOP of |ψ (1) remains 0 until h x increases to around 0.828, coinciding with the sharp change signaling in the the evolution of the parameters. Assuming that the virtual Z 2 symmetry is encoded in the first layer of GMPO, the VOP in Eq. (29) can be generalized to the higher-order GMPS by choosing W (s) = I ⊗(s−1) ⊗σ x . While this may not necessary be the case as the GMPO at different layers commute, we found that it works extremely well if we always adopt the optimized state |ψ (s) as the initial parameters for |ψ (s+1) . From Fig. 4(c), one can clearly see that the predicted critical points for the higher-order GMPS are getting closer to the exact point h x = 1. We note that while the VOP remains zero in the SSB phase, its value in the polarized phase does not converge to a fixed value when increasing the bond dimensions. This shows that the VOP can only be served as a qualitative probe to identify phases and cannot be used to extract the scaling behavior near the phase transition point. To study the critical behavior, one should adopt a more formal approach by identifying the entanglement order parameters (EOP) proposed in Ref. [25]. Remarkably, GTNO can also be used to compute the EOP even though we do not strongly impose the symmetry constraint of the local tensor as implemented in Ref. [25] (see App. A for details). This demonstrates GTNO's flexibility to combine with other established frameworks. Finally, we compare the variational energy with the exact energy via the relative energy deviation in Fig. 4(d). As expected, the relative energy deviation is efficiently lowered by considering the higher-order ansatz.

E. Case III: Transverse Field Cluster Model
The Hamiltonian of the TFC model is written as [35,36] which enjoys the Z 2 × Z 2 = a, b|a 2 = b 2 = e, ab = ba symmetry characterized by the two generatorŝ This model is exactly solvable and exhibits a phase transition between the symmetryprotected topological and polarized phases at h x = 1 [37]. Importantly, both phases respect the Z 2 × Z 2 symmetry and cannot be distinguished by any local order parameters. Instead, they should be distinguished by the string order parameters [38] and the degeneracy of the entanglement spectra [39,40]. Phrasing in the TN language, the corresponding virtual gauge transformation matrix V (g) of the symmetric MPS should form a non-trivial projective representations of the symmetry groupĜ [18,19,40].
While the interaction is three-site, one can treat them as the two-site interactions by blocking two spins into a single unit cell as follows: The generators of Z 2 × Z 2 symmetry then becomesÛ a = k (σ x ⊗ I) k , andÛ b = (I ⊗ σ x ) k . From Eq. (32), the sets of |A α ) and |B α ) with n = 2 can be easily identified as with the corresponding virtual transformation matrix Similarly, C β with n o = 4 can be identified as C 0 = I ⊗ I, C 1 = I ⊗ σ x , C 2 = σ x ⊗ I, and C 3 = σ x ⊗ σ x . For simplicity, we let C β with β = 1, 2, 3 on each A α B α share the same parameter. The resulting GMPO after imposing the normalization condition and inversion symmetry will then have 11 free parameters. Now, we choose the initial state as the product state polarized in the x-direction: |ψ 0 = ⊗ i |+ i , which is the only product state respecting the Z 2 × Z 2 spin-flip symmetry and having non-zero overlap with the true ground state for all h x ≥ 0. The resulting s-th order GMPS |ψ (s) =Ĝ s |ψ 0 then also respects the Z 2 × Z 2 symmetry. Fig. 5(a) shows the evolution of the variational parameters with respect to the transverse field for |ψ (1) . Similar to the TFIM, one can observe that all the parameters exhibit a sharp change at h x ≈ 0.828, signaling the occurrence of the phase transition. To identify the VOP, we first note that from Eq. (34), the virtual gauge transformation of GMPO V (i) = ⊕ α [V (i) ] α is the linear representation of the Z 2 × Z 2 symmetry. Therefore, the corresponding symmetric GMPS describes the topological trivial phases in general. To see how the non-trivial projective representation emerges in the SPT phase, we consider the zero field limit h x = 0. The projector to the ground state subspace can be represented by GMPO with all the parameters for the twosite interaction, i.e., A α B α C 0 , becoming equal, while all the parameters for the one-site interactions, i.e., A α B α C β =0 , vanish. In this extreme limit, one can show that the GMPO satisfies the following physical-virtual transformation: , with W a = σ x ⊗ σ z and W b = I ⊗ σ x . Since W a W b = −W b W a , it forms the projective representation of the Z 2 ×Z 2 symmetry. The VOP can then be defined as the difference between the original and transformed tensors as shown in [ Fig. 5(b)]. Using Eq. (13), one can derive that W a,(s) = I ⊗(2s−2) ⊗ W a for the s-th ordered GMPS. Fig. 5(c) shows the VOP as the function of h x for |ψ (s) , s = 1, 2, 3. The non-zero value of VOP signals the occurrence of the phase transition and the estimated transition point is moving closer to the exact critical point h x = 1 when considering the higherorder GMPS. We also compute the relative energy deviation, showing that higher order GMPS yields more accurate energy of the ground state [ Fig. 5(d)].

IV. GPEPO
In this section, we generalize the GMPO construction to higher dimensions. While we focus on the 2D square lattice, generalization to other geometries is straightforward. Given (A| = ⊕ n−1 α=0 (A α | and |B) = ⊕ n−1 α=0 |B α ), the GPEPO with bond dimension D = n−1 α=0 D α and n 4 free parameters can be constructed as or schematically, .
Similar to the GMPO, the outgoing and incoming bonds are always attributed to (A| and |B), respectively. This immediately guarantees that the physical transformation can be pulled through the virtual bonds (using Eq. (5)): .
Therefore, the global operatorĜ formed by contracting all the virtual bonds of G v1v2v3v4 inherits all the on-site symmetries from the Hamiltonian. By applying the GPEPO to some simple symmetric PEPS tensor, one can then construct the symmetric GPEPS ansatz to variationally study 2D systems. The incorporation of the onsite terms is straightforward as well: given the set of two-leg tensors C β = sis i (C sis i ) β |s i s i |} (β = 0, 1, · · · , n o − 1), we consider the GPEPO We note that while Eq. (39) seems to be complicate, one can easily write down a program to construct the GPEPO without knowing its explicit form. In Sec.IV A-IV C, we will use the GPEPO to study the 2D Heisenberg, TFI, and toric code [41] in a magnetic field (TCX) models. Note that one cannot apply Eq. (39) to the TCX model, as it consists of the four-site interactions. While we currently have no simple recipe to write down the general four-site GPEPO, we provide explicit construction of the GPEPO for the TXC model in Sec.IV C. Besides, we remark that the most unbiased initial state for GTNO is the simplest state respecting all the symmetries of the Hamiltonian. As demonstrated in Sec.III, the symmetric GMPS can still describe the SSB phases, where the degeneracy of the ground states is encoded in the virtual symmetries. However, to compare with other numerical approaches for the 2D cases, we will use the prior knowledge about the system to choose the initial state for 2D Heisenberg and TFI models, which is also assumed in Ref. [42] and [43] Different from MPS, PEPS has no exact canonical form and cannot even be contracted exactly. This renders an efficient optimization of the parameters in a PEPS to be much more challenging than the MPS. Recent progress has been made to effectively evaluate the gradient of the energy functional either diagrammatically [8,9] or using the automatic differentiation approach developed in the machine learning community [10]. We will adopt the latter approach using peps-torch  [44] to optimize the ansatz once the GPEPS has been established for the desired models. Peps-torch has been demonstrated with a very good capability for various spin systems such as frustrated Heisenberg antiferromagnet [42,45], chiral spin liquid [46], and many other quantum magnetism [47,48], as well as bosonic systems [49].

A. Case I: Heisenberg Model
Now we study the 2D anti-ferromagnetic Heisenberg model using the GPEPO evaluated through Eq. (36) with |A) = |B) = I σ x σ y σ z . Since the ground state develops Neel order, we choose the initial state as the product state polarized in the z-direction |ψ (0) = ⊗ i | ↑ i , and then rotate the physical σ z basis by unitaries −iσ y at one of the sublattices. The resulting GPEPS |ψ (n) =Ĝ|ψ (0) then acquires the U (1) symmetry, the subgroup of SU (2). To yield more accurate ground states of GPEPO for the given bond dimension D = 4, we note that since the ground states is not SU (2) symmetric, one can relieve the SU (2)-symmetric constraint on the GPEPO. Therefore, we consider the n = 2 with (A 0 | = (B 0 | = I , (A 1 | = (B 1 | = σ x σ y , and (A 2 | = (B 2 | = σ z . The resulting GPEPO constructed using Eq. (36) is then only invariant under the U (1) on-site transformationÛ (θ) = i e −iθσ z i /2 . Such GPEPO is also feasible for studying the XXZ model.
Tab. II shows energies and magnetizations obtained by two different GPEPOs and VarOpt. Here, the PEPS of VarOpt is also optimized using peps-torch with the only condition of imposing the C 4v symmetry. The GPEPSs, with extremely few parameters (5 and 16 for SU (2) and U (1) GPEPO, respectively), yield very close energies and magnetizations compared to the results by VarOpt, which contains 110 parameters. These outcomes provide a clear evidence that GPEPS faithfully capture the physics of the Heisenberg model using much less parameters.

B. Case II: Transverse Field Ising Model
We now turn to studying the 2D TFI model, whose Hamiltonian readsĤ where V = [1] ⊕ [−1] = σ z . This implies that the GPEPS will respect the global Z 2 symmetry only when the initial state |ψ 0 respects the symmetry, i.e., θ/π = 0.25. Eq. (41) motivates us to define the VOP as the norm difference between the original and symmetry-transformed tensor [ Fig. 6(a)]. Fig. 6(b) shows the evolution of θ/π, VOP, and S z with the change of h x for |ψ (1) =Ĝ|ψ 0 . The physical order parameter S z = σ z /2 estimates the phase transition point at h x ≈ 3.09. Besides, the parameter of the initial state θ/π sticks to a fixed value 0.25 for h x > ∼ 0.309, consistent with our previous claim that the GPEPS will respect the spin-flip symmetry only when |ψ 0 is invariant under the symmetry. The phase transition is also captured by the vanishing of VOP in the symmetry-preserved phase. We also consider the 2nd order GPEPS, and the results are similar to |ψ (1) with more accurate estimated critical point at h x ≈ 0.3051.
Tab. III shows the comparison between GPEPS with previous results using VarOpt [9] and the perturbation method [43]. For D = 2, GPEPS estimates the same critical point as the VarOpt did, which is far better than the results obtained using the perturbation method. However, we note that the number of parameters for GPEPS is only slightly lesser than that of  [9] and the perturbation method [43]. Numbers in the parentheses represent the number of variational parameters. The critical point is located at hc = 3.04438(2) as estimated by the Monte Carlo simulations [51].
VarOpt. To demonstrate the power of GPEPO, we also consider the 2nd order GPEPS. With merely 19 free parameters, GPEPS predicts a more accurate critical point than the VarOpt did for D = 3 using 42 parameters.

C. Case III: Toric Code in a Magnetic Field
The toric code model [41] can be defined on the square lattice associated with two types of plaquette, A and B, and the spin half degrees of freedom placed on the vertices [ Fig. 7(a)]. The Hamiltonian consists of the four-site interactions acting on the two plaquettes, and can be written as where h A = i∈A σ z i and h B = i∈B σ x i . Since the model is frustration-free, i.e., [h A , h B ] = 0, ∀A, B, the ground state subspace is spanned by the states satisfying h A |ψ = h B |ψ = +|ψ . The system is in the intrinsic topologically ordered phase which is characterized by the topological degeneracies and anyonic quasiparticle statistics [52]. Similar to the SPT phases, topologically ordered phases cannot be identified using any local order parameters. Instead, they should be detected by the topological entanglement entropy, features of the entanglement spectrum, and the modular matrices extracted from the full set of the ground states respecting anyonic excitations [53][54][55][56]. Recently, it has been shown that the topological order can be classified using the virtual symmetries of the PEPS tensor [20,21].
To study the topological phase transition, we subject the model to a magnetic field pointing in the x direction, which we call the TCX model. The Hamiltonian can then be written as The ground state belongs to the charge-free sector h B = +1 for all h x ≥ 0. Using this property, the model can be mapped to the 2D TFI model [57,58], and thus the phase transtion point can be identified at 1/2.044382 ≈ 0.3285. On the optimization side, we can focus on the charge-free sector and consider the Hamiltonian without h B . Therefore, we choose the initial state as |ψ (0) = ⊗ i |+ i , which is the only product ij , one can express the α-th power of the four-site interactions h (4) ijkl as the contractions of the four four-leg tensors: where the Einstein summation convention is assumed for the contractions of virtual indices v i = 0, ..., D α − 1 [ Fig. 7(b)].
Since h A is C 4v symmetric, we choose A α vv = B α vv = C α vv = D α vv . It is then straightforward to identify n = 2 and D α = 1, ∀α with A 0 00 = I, A 1 00 = σ z . Letting A vv = ⊕A α vv , the GPEPO with bond dimension D = 2 can be constructed as [see Fig. 7(c)]. Similarly, the one-site interaction can be incorporated into the GPEPO as with C 0 = I, C 1 = σ x . Enforcing the C 4v symmetry and the normalization condition, the non-zero action of GPEPO with 4 parameters can be written as Fig. 8(a) shows the evolution of the parameters as we increase the magnetic field for the first-order GPEPS |ψ (1) = G|ψ 0 . Since σ x i |ψ 0 = |ψ 0 , c 3 and c 4 can be discarded. The phase transition point at h x ≈ 0.3115 is clearly identified by the abrupt change of the parameters. To identify the VOP, we note that at h x = 0, the exact toric code ground state corresponds to the application of the GPEPO with c 1 = c 2 = 1. In this extreme limit, the GPEPO satisfies the extra virtual symmetry , where W = σ x . It turns out that this virtual symmetry is the manifestation of the Z 2 TO phases [20] [see App. A for details]. The VOP can then be identified as the difference between the symmetry-transformed and the original GPEPS. Fig. 8(b) shows the VOP as well as the physical observable h A and σ x . One can clearly observe the happening of the phase transition at h x ≈ 0.3115. We note that given the PEPS tensor respecting the virtual Z 2 symmetries, there are several ways to calculate the topological entanglement entropies [55,59], entanglement spectrums [60], and even extract the low-lying anyonic excitations [61,62]. The optimized GPEPS can be directly served as an input tensor of those framework without any gauge transformation, which demonstrates the compatibility of GTNO with other numerical schemes.
To get a more accurate ground state, we consider the second-order GPEPS with 4 additional parameters. The VOP and physical observables show similar behavior as |ψ (1) , and the critical point is identified at h x ≈ 0.3265, which is already rather close to the result of QMC. Tab. IV shows the comparison of energies at h x = 0.2, 0.36 and the estimated critical points obtained from GPEPS and other methods. For D = 2, with merely two parameters, GPEPS clearly outperforms the perturbation results (1 parameter) on locating the critical point. Furthermore, GPEPS yields comparable energy with the VarOpt method (roughly 32 parameters), indicating that GPEPS reduces the number of parameters efficiently without losing the entanglement structure for the TCX model. For D = 4, with six parameters, GPEPS has identified a more accurate critical point than the perturbation method (15 parameters) as well. Besides, it achieves a comparable energy with the VarOpt (≈ 54 parameters), and even locates the critical point more accurately than VarOpt. All of these comparisons strongly demonstrate the power of GPEPS. We also evaluate the charge condensation order parameters introduced in Ref. [25]. The slope of the scaling close to the critical point matches the critical exponent β ≈ 0.3265 of the order parameters of the 3D Ising transition, consistent with results reported in Ref. [25](see App. A for details).

V. DISCUSSION AND OUTLOOK
We have introduced the Generic-TNO as a unifying scheme to study the general quantum many-body systems. GTNO combines the idea of ITE that the ground state can be obtained by applying some operators generated from the Hamiltonian, and VarOpt by allowing free parameters to optimize the wave   [26] and the perturbation method [43]. The critical point is located at hc = 1/3.04438(2) ≈ 0.3285 as estimated by the Monte Carlo simulations [51] using the mapping between TFI and TCX models. function. Applying GTNO to simple initial states, the resulting GTNS can be used to obtain reliable ground states with far less parameters than the usual TNS. After the optimization, symmetries of the local tensor emerge without performing any gauge transformation, providing the probes to classify quantum phases.
We demonstrate the application of GTNO by studying several models in one and two spatial dimensions. For the 1D Heisenberg model, we obtain a SU (2) symmetric GMPS with extremely small energy deviation δE ≈ 1.31×10 −5 and large correlation length ξ ≈ 205 unit cells using only six parameters. For the 1D TFI (TFC) model, we adopt the Z 2 (Z 2 × Z 2 ) symmetric GMPS to search for the ground states and detect the phase transitions using the emergent symmetries of the local tensor. The power of GTNO is further manifested when considering the 2D systems. Applying the SU (2)(U (1)) symmetric GPEPO with merely 5(16) parameters to study the 2D Heisenberg model, we obtain comparable energy with the VarOpt methods using 110 parameters for the same bond dimension D = 4. For the 2D TFI model, we allow the GPEPS favors one of the degenerate ground states in the SSB phase to compare with other numerical approaches. This shows the flexibility of GPEPO to combine with the initial states containing variational parameters. While we currently have no general recipe to deal with the Hamiltonian consisting of foursite interactions, we explicitly construct the GPEPO for the TCX model. Remarkably, the GPEPS with merely six parameters and bond dimension D = 4 already estimate more accurate critical points than the VarOpt with 54 parameters for D = 3.
There are many possible extensions of the present work. First of all, in this paper, we focus on the Hamiltonian consisting of the Pauli matrices, as they all guarantee the closed condition that {(A α |B α )} with α = 0, 1, ..., n − 1 form a basis of any polynomial function of h (2) i,j . This property also holds for several models consisting of other operators, such as the Clifford algebra Γ α satisfying {Γ α , Γ β } = 2δ αβ [63] and the clock and shift matrices defining the quantum Z n clock model [64]. As a result, the construction of GTNO can be directly applied. Furthermore, we note that the requirement of the closed condition is to guarantee that GTNO can obtain exact ground state for the frustration-free Hamiltonian. For the model not fulfilling the closed condition, one can still use the power of h (2) i,j to generate more physical operators, and we believe that the GTNO framework can still yield accurate ground state with far less parameters than the usual TNS.
Second, GTNO can be used to study several exotic and numerically challenging models whose theoretical PEPS frameworks have been developed. For example, the general stringnet and 2D SPT models have been classified using the matrixproduct operator (MPO) symmetries in the PEPS framework [21,22]. However, variationally finding the optimal ground states with desired gauge symmetries is still challenging due to the multi-site interaction of the local Hamiltonian, which limits the current numerical study to the local-filtering method [65,66]. By generalizing the GTNO to these models, one should be able to accurately study their topological phase transitions.
Third, by taking advantage of the fact that the symmetries of the local tensor automatically emerge during the optimization, GTNO can be used to explore the models whose PEPS frameworks are not yet well understood. For instance, whether PEPS can faithfully represent the chiral shin liquids (CSL) is still under heated debate. While several efforts are devoted to representing gapped CSL using PEPS, they all exhibit infinite correlation length [67][68][69][70][71][72][73][74][75], consistent with the no-go theorem proven in the free-fermion case that the parent Hamiltonian of a chiral PEPS is gapless [76]. Nonetheless, in a very recent paper, Hasik et al. argue that a faithful representation of the CSL phase is in fact possible [46]. It is interesting to apply our GTNO framework to study the similar models in Ref. [46] and see whether there will be some extra symmetries emerging in the CSL phase.
Finally, we note that GTNO can also be used to study the Kitaev's honeycomb model [77], and one can show that the construction is similar to the dimer gas operator proposed in Ref. [78]. However, the distinguishing feature of the Kitaev spin liquids is that there is a constant of motion called vortex associated with each plaquette. The current GTNO cannot be used to project the wave functions to the vortex-free sector, and one needs the loop gas operator introduced in Ref. [78]. Figuring out how to systematically construct the vortex-free projector of the several generalized Kitaev models [79,80] and applying GTNO to probe the accurate ground state prop-erties are currently under investigation. In Sec.III D and Sec.IV C, we find that the virtual Z 2 symmetries of the local GMPS and GPEPS naturally emerge in the 1D SSB and 2D TO phases, respectively. The virtual order parameters (VOP) signaling the phase transitions are then identified as the norm difference between the original and symmetry-transformed tensors. Remarkably, we can defined another quantities, called entanglement order parameters (EOPs), to detect the phase transitions and further study their critical behaviors. The idea of EOP is originally introduced in Ref. [61] as a probe to model anyons' behavior by physically deforming the toric code wave functions to the trivial phases. Since then, it has been generalized to the Abelian double models [81], string-net models [66], and even been used to identify quantum spin liquids nature of the spin-1 Kitaev model [62]. Very recently, Iqbal et al. have further shown that one can use the EOP to extract the critical exponents near the phase transitions [25]. In the following we give an intuitive explanation on why the tensors satisfying this property, called Z 2 -invariant MPS and PEPS, form a natural framework to represent the 1D Z 2 SSB and 2D TO phases [20]. After that, we motivate the construction of the EOP [25] and show the results using the tensor optimized by GMPO and GPEPO. We refer interested readers to Refs. [20,25] for details.

Z2-invariant MPS
Given a Z 2 invariant MPS satisfying u g Au g = A with g = {I, X} and u g the representation of g, one can construct the parent Hamiltonian such that its ground state subspace can be obtained by attaching a single u g on the virtual bond. To simplify the notation, we denote u I (u X ) as I(X) and the operator that anti-commutes with X as Z in the following. The degeneracy stems from the fact that the parent Hamiltonian cannot distinguish those states locally, as one can use the Z 2invariant property to translate the position of u g . Besides,  we can create a domain wall excitation on top of the ground state by attaching a single Z operator on the virtual bond of the Z 2 -invariant MPS. Since the resulting tensor transforms anti-invariantly instead of invariantly under the virtual Z 2 action, one can show that no local physical operator can create a single domain wall, which is a feature of the topologically non-trivial excitation. We summarize the framework of Z 2invariant MPS in Fig. 9(a). For the s-th order GMPS studied in Sec.III D, X (s) = I ⊗(s−1) ⊗ σ x and Z (s) = I ⊗(s−1) ⊗ σ z in the Z 2 SSB phase.
However, Z 2 -invariant MPS does not necessarily guarantee the Z 2 symmetry breaking phase, as the domain wall ansatz we create can either be confined ( DW |DW = 0) or condensed ( DW |gs = 0). Therefore, one should further examine its normalization and overlap with the ground state. In other words, the system is in the Z 2 SSB phase only when the domain wall confinement order parameter K DW = DW |DW = 0 and the domain wall condensation order parameter C DW = DW |gs = 0. Using the fact that the single Z action can be used to create a domain wall, K DW and C DW can be transformed into the virtual observables as shown in Fig. 9(b). One then use those observables to detect the phase transitions and further study the critical behavior. Now, we use the GMPS optimized by GMPO in Sec.III D to calculate the EOPs. Throughout the calculations, we find that K DW remains nonzero for all h x . This is consistent with the fact that applying the transverse field will induce the phase transition only through the condensation. Fig. 10 shows the domain wall condensation order parameter C DW as a function of transverse field for the s-th order GMPS with s = 1, ..., 4. One can clearly observe that C DW signals the occurrence of the phase transitions, and the transition points coincide with the one identified using the VOP in the main text. Besides, instead of getting flatter like the VOP, C DW becomes sharper as the bond dimension increases, showing its advantage of studying critical behaviors over the VOP. We note that the GMPS does not enjoy the virtual Z 2 symmetry in the symmetrypreserved phase, and the physical interpretation of K DW and C DW in this regime is unclear. However, it is remarkable that the framework developed under the strict assumption of the virtual Z 2 symmetry can still be applied to GMPO. Besides, throughout the calculations, we do not perform any gauge fixing condition as implemented in Ref. [25]. All of these demonstrate the power and flexibility of our GMPO framework.

Z2-invariant PEPS
Parallel to the one dimensional counterpart, given a Z 2invariant PEPS satisfying A(u g ⊗ u g ⊗ u g ⊗ u g ) = A with g = {I, X} and u g the representation of g, one can construct a parent Hamiltonian such that its ground state subspace is spanned by two non-contractible loop operators on the torus. Similar to the Z 2 -invariant MPS, the degeneracy can be understood by noting that one can use the Z 2 -invariant property to translate the loop operators. Furthermore, we can create anyonic excitations of the parent Hamiltonian on top of any four ground states. The charge anyon |e can be created by attaching a single X operator, and the flux anyon |m can be created by a half-infinite virtual Z string on the virtual bonds of the Z 2 -invariant PEPS. Finally, the fermion | can be created using the braiding rule = e × m of the Z 2 QSLs. We summarize the framework of Z 2 -invariant PEPS in Fig. 11(a). For the s-th order GPEPS studied in Sec.IV C, X (s) = I ⊗(s−1) ⊗ σ x and Z (s) = I ⊗(s−1) ⊗ σ z in the Z 2 TO phase. Similar to the Z 2 -invariant MPS, Z 2 -invariant PEPS does not promise the Z 2 TO phase, and one should investigate the normalization of the different anyonic excitations and their overlaps with the ground state. Remarkably, since the condensation of the charge(flux) anyon is always accompanied by the confinement of the flux(charge) anyon, the TO phase can be identified if the charge confinement order parameter K e = e|e is nonzero and the charge condensation order parameter C e = e|gs vanishes. Now, we use the GPEPS optimized by GPEPO in Sec.IV C to calculate the EOPs. Throughout the calculations, we find that K e remains nonzero for all h x . This is consistent with the fact that applying large enough h x will induce the phase transition through the charge condensation. Fig. 12(a) shows the charge condensation order parameter C e as a function of magnetic field for the 1st and 2nd order GPEPS. As expected, non-vanishing C e indicates the happening of the phase transition, and the estimated critical fields are consistent with the ones using VOP in Sec.IV C. To demonstrate that GPEPS can also be used to study the critical behavior even if it does not respect the virtual Z 2 symmetry in the charge condensed phase, we perform the log-log scale of C e and h x in Fig. 12(b). We find that both the slopes of 1st and 2nd order GPEPS match the the critical exponent β ≈ 0.3265 of the order parameters of the 3D Ising transition. This is consistent with the results reported in Ref. [25] even if we do not strictly impose the Z 2 symmetry and have extremely few free parameters (two and