Stacked tree construction for free-fermion projected entangled pair states

The tensor network representation of a state in higher dimensions, say a projected entangled-pair state (PEPS), is typically obtained indirectly through variational optimization or imaginary-time Hamiltonian evolution. Here, we propose a divide-and-conquer approach to directly construct a PEPS representation for free-fermion states admitting descriptions in terms of filling exponentially localized Wannier functions. Our approach relies on first obtaining a tree tensor network description of the state in local subregions. Next, a stacking procedure is used to combine the local trees into a PEPS. Lastly, the local tensors are compressed to obtain a more efficient description. We demonstrate our construction for states in one and two dimensions, including the ground state of an obstructed atomic insulator on the square lattice.


I. INTRODUCTION
Tensor network (TN) representations, like the matrix product state (MPS), provide highly efficient descriptions of quantum many-body states.An MPS is always disconnected when any of its bonds is cut.Physically, this implies the virtual Hilbert space attached to any bond in the network could be given a simple interpretation in terms of the bipartite entanglement of the state [1,2].This interpretation is valid whenever the tensor network is free of loops, and it has two major consequences.First, a loop-free TN state is constructible, in the sense that one could directly construct the TN representation of any given state, up to an error threshold, by successive bipartitions of the system [1][2][3][4].Second, a loop-free TN state is also computable, in that there exist canonical forms which enable the numerically exact contraction of TN diagrams arising from, for instance, the computation of physical observables [1,2,5].
Though powerful, the MPS ansatz is natural only for one-dimensional (1D) systems (or as a quasi-1D modeling of higher-dimensional systems).Applying the ansatz in higher dimensions through either a 1D ordering of all the sites [6][7][8] or through its extension to a tree TN state [9][10][11], which remains loop-free, would unavoidably assign certain physically neighboring sites to far-apart nodes on the network.Such TNs are generally incompatible with the physical locality of the state and therefore cannot efficiently encode part of the short-range entanglement in the system.
The projected entangled pair state (PEPS) is another natural generalization of MPS to higher dimensions in which physical locality is retained at the cost of introducing loops [2].As a result, a PEPS is generally neither constructible nor computable: variational * hcpo@ust.hkoptimization or imaginary-time Hamiltonian evolution is needed for finding the PEPS representation of a state, and approximations are invoked in evaluating physical observables through TN contractions [12].
The recent proposal on isometric TN states [13,14] has provided a fruitful avenue for attacking the computability problem of PEPS.Here, we seek to address the complementary problem concerning constructibility, namely, can one obtain a PEPS representation directly from a given state?This question had been answered in the affirmative for the ground states of special models, like those corresponding to stabilizer codes [15][16][17].However, for more general problems, indirect approaches like optimization or imaginary-time evolution remain the only tenable options so far.This is true even for free-fermion states, as is reflected in the recent bodies of work concerning fermionic Gaussian TN states [18][19][20][21][22][23].
In this work, we demonstrate the PEPS constructibility of the ground state of a free-fermion obstructed atomic insulator in two dimensions [24] Our approach follows a divide-and-conquer strategy and consists of three steps, in the order of "tree, stack, and compress."First, we derive the tree TN representations for the local descriptions of the state over small open disks.Next, we stack the tree TN states to cover the full two-dimensional space.Importantly, the patches overlap and so the resulting TN takes a PEPS form.Lastly, we compress the local tensors to obtain an efficient representation.This is achieved by applying MPS techniques to the partial contractions of the TN state along one-dimensional subregions.
We remark that, as a proof of principle, we consider here the ground states of translation-invariant free-fermion Hamiltonians.This allows us to shortcut some of the analysis through band-theory techniques like Wannierization.
Our approach can be readily generalized to any free-fermion state admitting a (possibly approximate) localized description in terms of filled Wannier functions (WFs).In our current formulation, however, the "stacking" step makes explicit use of the Gaussian nature of the local tensors; generalizing this step to an interacting state is likely nontrivial.Nevertheless, the free-fermion TN representation we constructed could still serve as a natural starting point for constructing an interacting fermionic TN state through, for instance, the Gutzwiller projection [25][26][27][28].

II. SETUP
We begin by explaining how our "tree, stack, and compress" steps are carried out for a free-fermion state.We consider a free-fermion Hamiltonian Ĥ = ij h ij ĉ † i ĉj , where ĉ † i and ĉi are respectively the free-fermion creation and annihilation operators.The subscript i denotes possible degrees of freedom, like physical sites, orbitals etc.The ground state |Ψ⟩ of Ĥ, as is the case for any fermionic Gaussian state, is fully determined by its two-point correlation functions [29].
With the number conservation symmetry in our context, we only need to focus on the correlation matrix We further specialize to the case that Ĥ is translationally invariant and |Ψ⟩ can be obtained by the filling of a full set of WFs, which corresponds to an atomic insulator.The WFs can be viewed as a particularly suitable choice of Fourier transform of the filled Bloch states such that they become exponentially localized in the real space [30].For an atomic insulator, the WFs can be chosen such that they further respect all the internal and spatial symmetries of the system.Tree decomposition.Our first step is to obtain a tree TN description for the ground state over a small local subregion.This can be achieved by first focusing on a single WF, which represents a locally defined fermion mode that is occupied in the ground state.Given the exponential localization, the WF can be well-approximated by a truncation to a disk of some radius r trunc which is on the order of its localization length.This is illustrated in Fig. 1 (a), where the WF centered at the blue dot is picked, and the truncation is indicated by the green circle.
To obtain a TN representation of the truncated WF, we define a tree which specifies how the sites in the region are to be connected in the TN.For instance, as demonstrated in Fig. 1 (b), we can grow a tree with 4-fold rotation symmetry on the square lattice.
As a tree is loop-free, we can convert the wave function into a tree TN form by successively applying Schmidt decompositions.More concretely, we view the center site of the tree, which coincides with the center of the WF, as its root.Note that the center need not be occupied by a physical site, and so there may not be a physical leg attached to the center (Fig. 1).Any other sites can be given a height according to its distance from the root.We say two sites belong to the same level if they are equidistant from the root.Starting from the highest level, we perform Schmidt decomposition to obtain the local tensors defined on the sites in the level.Schmidt decompositions within the same level are independent.
For free-fermion states, Schmidt decomposition can be done at the level of correlation matrices, since the reduced density matrix of a subregion is still Gaussian and so it shares the same eigenbasis with the restricted correlation matrix.More explicitly, consider a bipartition of the system into A and B. We can diagonalize the restricted correlation matrices of C as: where  After Schmidt decomposition, one expects to obtain two free-fermion tensors T A and T B which can be contracted to reproduce the original state |Ψ⟩.The meaning of a free-fermion tensor, however, is unclear as there will generally be multiple legs with varying number of fermion modes attached to them.For bosonic systems, like qubits, a tensor can always be reinterpreted as a state through a mere reshaping of the legs; for fermions, care must be taken to ensure such reshaping is done in a consistent manner [32].In our formulation, this is achieved by purifying all the unitary operators into fermionic Gaussian states defined on a doubled space [33].In Appendix B, we show how the data contained in Eq. ( 1) can be packaged into two free-fermion states |T A ⟩ and |T B ⟩, which reproduce |Ψ⟩ upon contraction.This way, all the local tensors can be interpreted as free-fermion states.
As such, in the following we use "local tensors" interchangeably with the correlation matrix of its corresponding free-fermion state.
Upon performing all the contractions of the local tensors, as represented by the edges on the tree in Fig. 1(b), we reconstruct the single-particle state given by filling the original (truncated) WF in the current local subregion.Stacking.To reconstruct the full state |Ψ⟩, we would need to combine the locally defined tree TN states obtained from the individual WFs.Intuitively, we simply need to consider the collection of all the tree TN states, which in our context are related to each other through translation symmetry, and show that these states can be recombined into a single PEPS (Fig. 1 (c)).This step, referred to as the "stacking" procedure, can be achieved as follows.Suppose the local tensor C i,α at site α decomposed from the i th tree is represented as: where the correlation matrix is organized with respect to the physical legs and bonds.Here, i indexes the set of truncated WFs which have support on the site α, and the subscripts p vs. v indicate whether the fermion modes are associated with the physical or the virtual legs.
The stacking of the local tensors C i,α for the site α could be expressed as: where we suppose that there are m trees that contribute to the site α.
Here, the virtual spaces from different trees are independent and so the corresponding parts of their correlation matrices are simply combined as a direct sum.However, all the trees share the same physical Hilbert space at site α and so their contributions add up.As defined, Cα is not a proper correlation matrix in general, as the summing procedure defined above does not correspond to any well-defined operations on the Hilbert spaces concerned.Nevertheless, a proper correlation matrix C α can be obtained by the deformation procedure described in Appendix B 4. Intuitively, the failure of Cα to be a proper correlation matrix stems from the fact that the restrictions of the different trees to the physical Hilbert space of site α lead to modes which are not orthogonal to each other.The deformation process can then be simply interpreted as a suitable orthonormalization step.
We thus obtain a free-fermion PEPS defined by the collection of deformed local tensors.Upon contracting all the virtual legs, we obtain an approximation of the ground state |Ψ⟩.Compression.The approximate PEPS representation obtained from stacking, however, is far from optimal.In combining the individual trees, we treated their virtual Hilbert spaces as independent.
This leads to a superficially high bond dimension which grows as we increase the truncation radius r trunc used in approximating the WFs (Appendix A 2).As a last step, therefore, we perform a compression of the local tensors.
The idea is that we could first contract the local tensors in one direction to form a free-fermion state defined on an open 1D chain.We can then perform another MPS decomposition of the state while retaining only the most significant virtual modes, which correspond to a truncation to the bond dimension.More concretely, the virtual modes are retainiend according to their contribution to the von Neumann entanglement entropy: [29,[34][35][36][37], where ω i corresponds to i th diagonal entry for matrix Ω AA or Ω BB in equation 1.To reduce the bond dimension, we drop virtual modes that contribute the least to the entanglement entropy, i.e., we drop the mode i if ω i < ϵ or ω i > 1 − ϵ for a prescribed small threshold ϵ.Physically, these dropped virtual modes correspond to degrees of freedom that are well-localized within one side of the entanglement cut, i.e., they do not mediate entanglement across the cut and can therefore be dropped.
After successive Schmidt decomposition, local tensors deeply embedded in a long enough 1D chain should regain bulk properties.Therefore, we choose the local tensor in the middle as the updated C α .Repeating the above process along all possible directions for multiple times, the final fully-compressed C α is obtained, see Fig. 2.

III. EXAMPLES
We now move on to demonstrating our construction to two obstructed atomic insulator (OAI) systems as a proof of principle.OAI is a special class of atomic insulators for which the centers of the WFs cannot be chosen to coincide with any atomic sites in the system [38].
The 1D Su-Schrieffer-Heeger (SSH) model, one of the most well-known model for a topological insulator, can also be viewed as an OAI if only inversion but not chiral symmetry is retained.
The Hamiltonian for the SSH model is Ĥ = x (t − s)ĉ † O1,x ĉO2,x + (t + s)ĉ † O2,x ĉO1,x+1 + h.c., where t is the uniform hopping parameter, and s is a staggering between intra-and inter-cell hoppings.O 1 and O 2 are two different orbitals in a unit cell.In the numerics, we choose t = −1, s = 0.1.
As a second example, we construct an OAI model on the 2D square lattice protected by the four-fold rotation symmetry C 4 .
We assign three fermion modes, corresponding respectively to s, d x 2 −y 2 and p x + ip y atomic orbitals, to each of the sites.Our model is constructed by lowering the energy of a set of non-orthogonal "quasi-orbitals" f † which transform differently from all of the atomic orbitals in the system [38][39][40].This leads to a band insulator for which the WFs of the filled band are equivalent, symmetry-wise, to the quasi-orbitals we started from.More concretely, the Hamiltonian is , where the mode f † ⃗ R is localized to the center of the plaquette in unit cell ⃗ R, and transforms trivially under the C 4 rotation symmetry.An atomic insulator obtained by filling s-like WFs localized to the centers of the plaquettes is topologically distant from the innate atomic insulators in the Hilbert space, and so this Hamiltonian serves as an OAI model (Appendix A 1).In the numerical calculation, we set η = 2.
The main results are tabulated in Table I.The number of fermion modes attached to a virtual leg is denoted by b.For the 2D model (Fig. 1), we distinguish the vertical/horizontal bond b v,h and the diagonal bond b d for each square tensor.Before and after compression, we compare the difference between the correlation reproduced from PEPS C TN and the exact through the entry-wise maximum difference max i,j (|C TN − C exact | ij ).Small values of this difference, on the order of 10 −3 , were found, and so the exact and TN results are indiscernible in Figs. 3 and 4. We also compare the relative difference of ground state energy density, defined as (e TN − e exact )/δE gap where e TN and e exact are ground state energy density per unit cell obtained from TN representation and exact calculations respectively.δE gap is the minimal gap over the first Brillouin zone.

IV. DISCUSSION
In this work, we present a general scheme for constructing PEPS for free-fermion states arising from the filling of exponentially localized WFs.As a proof of principle, we demonstrate our approach for models in one and two dimensions.
Although translation invariance was used to simplify the computation, our approach can be generalized to a strictly real-space formulation for more general systems with incommensurate order or disorder.Interaction ĉ † ⃗ r 0 ,s ⟩) on a 100×100 square lattice, where ⃗ r0 = (0, 0).Data for ⃗ r along direction [11] within the range (−10, 10) and (10, 10) is demonstrated.R(⟨ĉ ⃗ r,l ĉ † ⃗ r 0 ,s ⟩) is 0 for ⃗ r outside the region.The red, green and blue colors represent the s, d x 2 −y 2 and p+ = px + ipy orbitals respectively.The dashed lines are for exact values and the markers show the reconstructed results obtained from our TN.
effects could also be incorporated by combining the free-fermion TN state with, for instance, Gutzwiller projectors [25][26][27].In closing, we remark that, in the restricted context of free-fermion states, there might be tantalizing connections between our construction and a tensor-network-based solution to the quantum marginal problem [41][42][43][44][45][46][47]: data confined to small local subregions are first handled by tree TN states, which are then patched into the full pure state through the stacking step.It is an interesting question to consider how our approach might be generalized to attack the corresponding many-body problem.
We construct a minimal 3-band "flat-band" model for an OAI on the square lattice in the main text.The word "flat-band" indicates that we only consider the dispersion of the filled band, which is obtained by filling s-like WFs, leaving the rest of bands dispersionless.Given the space group P 4 , we derive a model with the band representation of the ground state induced from the Wyckoff position 1b : ( 1 2 , 1 2 ).atomic sites of the square lattice are at Wyckoff position 1a : (0, 0).This is required by the definition of OAI, whose band representations should be induced from the unoccupied Wyckoff positions [38].Then by comparing the irreducible representations of the little groups at high-symmetry points in momentum space, we only need three different kinds of atomic orbitals from Wyckoff position 1a, namely s, d x 2 −y 2 and p x + ip y .
To obtain the Hamiltonian, we could start from constructing a quasi-orbital f † ), consisting of atomic orbitals from the nearest atom sites, as shown in Fig. 6.Column vectors w † ⃗ r indicate the local hybridization for each quasi-orbital at atomic site ⃗ r.The basis of column vectors is ĉ † s,⃗ r , ĉ † d x 2 −y 2 ,⃗ r and ĉ † px+ipy,⃗ r respectively from top to the bottom.The expression of the quasi-orbital using column vectors is f † ) for the unit cell at the origin, where ⃗ r 1 = (0, 0), ⃗ r 2 = (1, 0), ⃗ r 3 = (0, 1) and ⃗ r 4 = (1, 1).
With the quasi-orbitals all over the lattice, we construct the minimal Hamiltonian: where we summing over all unit cells ⃗ R. In the numerical calculation, we have η = 2, α = 1, β = 3 and γ = 5.The band structure across the first Brillouin zone is shown in Figure 5, where it is evident that a significant gap exists throughout the zone, providing the evidence that the model represents an insulator.

The Wannier Function
We could solve the model in eq.A1 for eigenvalue spectrum and the Bloch states ψ l ( ⃗ k) , by first applying Fourier transform and then eigenvalue decomposition.However, a direct Fourier transform to ψ l ( ⃗ k) cannot guarantee a set of well-localized WFs |Ψ l (⃗ r)⟩ in the real space, as WFs are well localized only when the Bloch functions possess a smooth gauge [30].We could smoothen the phase associated to ψ l ( ⃗ k) during Fourier One cell of the square lattice.The balck dot O is the origin, and the coordinate system is built according to the directions of the red arrows.Yellow dots at vertices are the atoms in the square lattice, which are at the Wyckoff position 1a : (0, 0).The blue dot is the center of the WF for the current unit cell and is at the Wyckoff position 1b : ( 1 2 , 1 2 ).The vectors are local wavefunctions w † ⃗ r for each atom site ⃗ r with the basis s, d x 2 −y 2 and px + ipy respectively.i.e. for the atom at ⃗ r = (0, 0), the wavefunction ŵ † transform by using trial functions ϕ ⃗ R (⃗ r) to perform a projection.ϕ ⃗ R (⃗ r) is obtained through a rough guess of the distribution of |Ψ l (⃗ r)⟩.For instance, we choose the quasi-orbital in Appendix A1 as the trial function for the 2D OAI model.We first Fourier transform ϕ ⃗ R (⃗ r) into ϕ l ( ⃗ k) , then the Bloch states with smooth gauge ψl ( ⃗ k) could be attained: Thus, WFs are: where Ω is a normalization factor.
However, it is worth noting a technical subtlety that the projection fails when the overlap ψ l ′ ( ⃗ k) ϕ l ( ⃗ k) between the trial function and Bloch state is 0. Additionally, since we use exponentially-decaying non-compact WFs, we set the threshold to 0.1 during the projection.If min( ψ l ′ ( ⃗ k) ϕ l ( ⃗ k) ) > 0.1, We assert that the symmetric WFs with exponential tails have been successfully obtained.
We plot one WF centered at 0 and ( 1 2 , 1 2 ) for 1D and 2D model respectively, as shown in Figure 7 and 8.It is observed that WFs for each orbital undergo rapid decay while retaining symmetry.
Due to the exponential localization of WFs, it is feasible to truncate them within a reasonable radius.However, truncating WFs necessitates a tradeoff between bond dimension and accuracy in the construction of PEPS.Generally speaking, larger truncation radii lead to larger bond spaces and more accurate results, as evidenced in Tables II and III

Appendix B: Operations on fG States with Diagrams
Fermionic Gaussian (fG) states could be expressed using density operators in the thermal form:   where Ĥ is a fermionic quadratic Hamiltonian and Z = Tr e − Ĥ .
The properties of fG states are totally defined by the correlation functions.Since Wick's theorem holds for fG states, it suffices to only take the two-point correlation functions into account.Furthermore, we impose the number-conserving condition, thus Ĥ becomes a free-fermion Hamiltonian: is sufficient in describing fG states, as the two correlation functions are related through For simplicity, when we refer to the correlation functions later in the text, we explicitly refer to C ij .For a pure state |Ψ⟩, the correlation function could be expressed as C ij = ⟨Ψ|ĉ i ĉ † j |Ψ⟩, while more generically, the correlation function with respect to the density matrix ρ of the state is C ij = Tr ĉi ĉ † j ρ .The fG Hamiltonian could be rewritten in terms of the correlation function: which indicates that the eigenvalues of C are bounded between 0 and 1.In the extreme case, when the eigenvalue is 0 (1), the corresponding eigenstate is completely filled (empty) since it should always (never) be occupied according to the energy.

Reshaping Operators
Operators acting on states could be equivalently represented as contraction among states, whereby operators are reshaped into states with suitable auxiliary space.This could be understood in analogue to the thermal field double in general bosonic cases, wherein the bras are reversed to kets making density operators reshaped into a state.Meanwhile, this process of reshaping could also be understood as purification in the context of quantum information [33].However, in the fG case, more care should be considered with appending auxiliary space.

a. Bosonic Case
Reshaping is straightforward in bosonic systems, if there's no further restriction.For instance, we can simply choose the auxiliary space to be identical as the original Hilbert space.To illustrate, we start from purifying an identity operator 1 to an identity state |1⟩: If the strategy works, the contraction K(|1⟩ , |1 ′ ⟩) between two identity states |1⟩ and |1 ′ ⟩ should be an identity state, since 1 |1⟩ = |1⟩.Tracing off the auxiliary space, we obtain the contraction result: which is as expected.However, the expression of |1⟩ is not unique, since the unitary transformation of basis |i⟩ will change the expression of |1⟩.To hold the consistency, we fix the basis in the reshaping of the identity, and derive the reshaping of operators based on the reshaped identity, using the relation for any arbitrary operator Â: The reshaping is demonstrated schematically in Figure 9.The rules of TN diagrams are as following: dashed lines indicate the contraction; bare solid lines represent identity operator; bare lines with arrows are identity states; squares with arrows out are states, while circles with arrows in and out are operators.

b. fG case
Following the same steps in Appendix B 1 a, we start from the reshaping of a fG identity operator.According to eq.B1, the identity for a m-fermion system is ρ1 = 1/2 m , which is a maximally entangled state.The corresponding correlation function is C 1 = 1/2.In contrast to the bosonic case, simply doubling C 1 does not give a pure fG state.To satisfy the requirement of a number-conserving fG state, we have to put extra constrains, i.e. the eigenvalues of C |1⟩ should be 0 or 1.The ways of purification are not unique.The convention we stick to is: where the first part of basis corresponds to the physical legs, while the second part belongs to the auxiliary space.Then based on equation B5, we could reshape any unitary or isometric operator Û in complex-fermion basis: The schematics in Figure 9 also work for the fG case, while all squares representing fG states.The proof is in Appendix B 3.

Schmidt Decomposition for Pure fG States
Bipartition of a pure state |Ψ⟩ into two subsystems A and B could be represented as a Schmidt decomposition: FIG. 10.The Schmidt decomposition of a state Ψ.The circled UA and UB are isometries correspond to the subsystems A and B after one bipartition.The squared UA and UB are purified states representing the subsystems A and B respectively.The Ψo contains the Schmidt coefficients of the Schmidt decomposition.
where λ i is the real and non-negative Schmidt coefficient, which is valid measure of the entanglement strength.i A(B) belongs to orthonormal sets of subsystem A (B). Schmidt decomposition is equivalent to the singular value decomposition (SVD): Ψ = U A ΣV † B mathematically, where U A (V † B ) contains an orthonormal set for subsystem A (B). Σ is a diagonal matrix with diagonal entries the Schmidt coefficients.
To see the relation between the diagonalization of a submatrix C A of C in subregion A and Schmidt decomposition, we consider the reduced Then it is clear the eigenbasis of C A are just the Schmidt basis of ρ and Schmidt weights λ i can be computed as Given the correlation matrix of a pure fG state C, the eigenvalues of C is either 0 or 1, i.e.U † CU = 1 0 0 0 , where U is the unitary that diagonalizes C with the columns referring to each eigenfunctions respectively.Thus, the correlation matrix could be rewritten using the eigenfunctions Φ with eigenvalues 1 as C = ΦΦ † , where those modes in Φ are referred to as the empty modes in our convention.Dividing elements in Φ into two parts based on their basis belonging to subsystem A or B, the state |Φ⟩ and correspondingly the correlation C become: We could further apply SVD to wavefunction Φ A and Φ B respectively, where: Substituting SVD result into eq.B11, we could obtain the expression for the Schmidt decomposition: (B14) However, there exists a phase ambiguity.Although the diagonal entries of the middle matrix are semi-positive and thus fixed, the off-diagonal entries of the middle matrix are not.
We claim canonical Schmidt decomposition by requiring the off-diagonal entries non-negative, which could be substantiated by imposing V A = V B .Thus, the canonical Schmidt decomposition is: which corresponds to the first equality in Figure 10.

Contraction
In this subsection we first present the contraction formula of fG states at the level of correlation matrices, and then show a brief derivation of the formula with Grassmann variables.For a fG state ρ, which lives in the composite Hilbert space H c ⊗ H c ′ , and another fG state in H c , we could obtain the fG state ρc ′ belonging to the Hilbert space H c ′ by taking contraction between the state ρ and another fG state ρc in the space H c .The contraction is equal to Tr c [ρ c ρ] up to a normalization factor.This kind of contraction is referred as "complete contraction", meaning the all degrees of freedoms in ρc are traced out.We further define more general "incomplete contraction" by first regarding ρc ′′ as in a different Hilbert space H c ′′ and taking tensor product ρ′ = ρ ⊗ ρc ′′ ∈ H c ⊗ H c ′ ⊗ H c ′′ , and then utilizing complete contraction with a "contraction kernel" state C kernel in H c ′ ⊗ H c ′′ .The contraction kernel is chosen such that the complete contraction of ρ and ρc is equal to the incomplete contraction of ρ and particle-hole dual of ρc ′′ .In most cases of this paper we use incomplete contractions, which is consistent as a reverse operation of Schmidt decomposition, i.e., the two states obtained from Schmidt decomposition can be contracted back to the full state.
In the language of correlation matrix, the complete contraction between ρ and ρc is given as : With this complete contraction formula, we can do incomplete contractions of two fG states C ρ and C ρc ′′ by first taking the direct sum C ρ′ = C ρ ⊕ C ρc ′′ and then applying complete contraction of C ρ′ and C kernel , which is the particle-hole dual of C |1⟩ : Next, we could check the equality in Figure 9 (a), where we employ the complete contraction between C ρ′ = C |1⟩ ⊕ C |1⟩ and C kernel : Then we can calculate: Similarly, we could also verify the contraction in Figure 10, where we apply the complete contraction of where the final result is exactly the expression in equation B15.Now we present a derivation of the formula with Grassmann integrations, based on the Grassmann representation of partial trace in [48].Consider a density matrix ρ of m fermionic degrees of freedoms, the Grassmann representation is where η is a column vector of real Grassmann numbers, and M ij = − i 2 ⟨[γ i , γj ]⟩ is usually called covariance matrix with Majorana operators defined through ĉk = 1 2 (γ 2k + iγ 2k+1 ).By doing a bit algebra we can turn to a complex Grassmann representaion, such that where ξ is the complex conjugate of ξ.For two states ρ, ρ c as defined before, the complex Grassmann representation are: (B27) Putting the form of C ρ and C ρc together, the right hand side turns out to be the following Gaussian integration: Using Gaussian integration formula, we obtain a Gaussian state which is represented by For the case that ρ c is pure, we can further simplify this result by noticing that (2C ρc − 1) 2 = 1, so it becomes the contraction formula eq.(B16).

Deformation
In general the correlation matrix C = ĉi ĉ † j of a fermionic state has eigenvalues between 0 and 1.In particular, if the state is pure, then the eigenvalues of the correlation matrix have to be 0 or 1 if and only if the state is Gaussian.In the main text, we met a question of stacking the local tensor of Wannier states to get a local tensor of the whole ground state, which physically means taking a product of these Wannier states.At the level of correlation matrices, it corresponds to the sum of each individual correlation matrices, each regarded as being embedded in a local Hilbert space of all physical and bond degrees of freedom.But the literal sum of two correlation matrices may not be a valid correlation matrix for a physical state, i.e., the eigenvalues of the sum could be negative.To obtain a physically valid correlation matrix of a pure free-fermion state, we need a deformation to make sure all eigenvalues that are not equal to 1 to be exactly equal to 0.
The reason why the sum of two correlation matrices is not a valid correlation matrix is the possible non-orthogonality of the two states to be summed.To see this, let C 1 and C 2 be two correlation matrices of pure fG states.The fermionic degrees of freedom are ĉ † 1 , ..., ĉ † k , k ∈ N + .Then the eigen-modes of C 1 could be represented by some normalized row vectors Ψ 1 , such that where r 1 is the rank of C 1 .Each row vector in Ψ 1 denotes an eigen-mode d † i = j m i,j ĉ † j , so the Dirac symbol of the corresponding product state has coefficients being the Slater determinants of any r-column submatrix of Ψ 1 , or say, the minors of order r.Similar one can obtain row vectors Ψ 2 of rank r 2 for the correlation matrix C 2 .We see that each time we add a row vector to Ψ 1 , we are adding an occupation of a mode to the state |Ψ 1 ⟩.However, for the relation C = 1 − Ψ † Ψ to work, there is a necessary condition: all normalized row vectors should be orthogonal to each other (we assume that all row vectors are already linearly independent).For instance, given two modes d † Since Q is unitary, d † j ≡ Q l,j ĉ † l is a canonical basis, thus this state is equivalent to i d † i up to a normalization.We also see from the QR decomposition that the if the row vectors Ψ are not orthonormal, the Ψ † Ψ = Q * R * R T Q T will have the same non-vanishing eigenvalues as R * R T which are only guaranteed to be positive but not necessarily equal to identity.
Therefore, we see the naive literal "sum" of C 1 and C 2 will not work since in general the eigen-modes of C 1 and those of C 2 are not orthogonal.Suppose the joint row vectors Ψ = Ψ 1 Ψ 2 has rank r 1 + r 2 < k, then we can find the eigenbasis of C = 1 − Ψ † Ψ, that is, C = U † Σ C U with eigenvalues in Σ C = diag[λ 1 , ..., λ r1+r2 , 1, ..., 1].Notice that the eigenvectors corresponding to eigenvalue 1 will be the empty modes, which will not be changed under the Gram-Schmidt process.So taking the eigenvectors except those corresponding to the eigenvalue 1 as the Ψ, we can obtain the a valid correlation matrix representing the product state |Ψ 1 ⟩ |Ψ 2 ⟩.This is equivalent to deforming all the non-unity eigenvalues of C to zero.
Numerically one may worry that the exact unity may be not easy to distinguish, if some occupied modes of non-orthogonal Ψ make some eigenvalues of C = 1− Ψ † Ψ close to identity numerically.Luckily this situation will not happen in the main text stacking the local tensors from Wannier functions.Physically the eigenvalues of Ψ † Ψ will be close to zero only when some row vectors in Ψ are close to be collinear, i.e., ⃗ m 1 ≈ λ ⃗ m 2 , λ ̸ = 0.A Wannier function of the state is not translational invariant, so local tensors at different locations are not possibly collinear.In fact, the spectrum of eigenvalues in the stacked tensors turn out to exhibit clear gaps separating the unity and non-unity values.Thus there is no numerical ambiguity in the deformation process.

FIG. 1 .
FIG. 1.The overall procedure of the "tree, stack and compress".(a) Exponentially decaying WFs centered at the centers of plaquettes of the square lattice, while atomic sites are the vertices.(b) The tree decomposition of a single WF, where the center of the WF is indicated by a blue dot.The green circle encloses the truncated region for the WF.The dashed lines indicate the square lattice, while the solid lines are legs of the tensors.The legs pointing out-of-page are the physical legs, whereas those in-plane are virtual legs/bonds of local tensors.Blue legs connect the diamond in the center to the neighboring sites, and black legs connect among the square tensors.(c) PEPS obtained by stacking the trees over the whole lattice using translation symmetry.Note that both blue (between square and diamond tensors) and black (among square tensors) are present along the diagonal direction.(d) The local tensors all over the lattice after the compression.

FIG. 2 .
FIG. 2. Compression procedure.The thickness of the legs indicates the size of the bond space.Red dashed lines indicate the successive Schmidt decomposition performed along the red arrow.The hollow arrow points to the final result of each compression step, which is the middle piece among all the decomposed local tensors.The ordering of compression direction is first along (a) horizontal direction, then (b) vertical direction and lastly along the two diagonal directions (c) and (d).
FIG. 3. Correlation function ⟨ĉ x,l ĉ † 0,O 1 ⟩ for a 100-site SSH model.Only data inside the range (−20, 20) is shown as the values are practically zero outside of this range.The blue and red colors represent the two orbitals O1 and O2 respectively.The dashed lines represent exact values while the markers represent the compressed ones obtained from the constructed TN.

FIG. 5 .
FIG. 5. Band structure of the 2D OAI model on a 100 × 100 lattice.The smallest band gap over the first Brillouin zone is 0.2286.Bands are shifted according to the Fermi energy -0.1143.
and from the discussion of Gaussian state, ρA = 1 Z A m e −ξm ĉm ĉ † m where ξ m are the eigenvalues of restricted correlation matrix in subregion A, since the restricted correlation matrix and density matrix are diagonalized simultaneously.

the product state is d † 1 d † 2 =
N −1 (−a 2 b 1 |110⟩ + a 1 b 3 |101⟩ + a 2 b 3 |011⟩),with normalization factor N = a .It is clear to get the corresponding correlation matrix C, which is which does not work.Instead, the proper row vectors should be the ones after doing the Gram-Schmidt process and normalization, consider a set of non-orthonormalized modes { d † 1 , d † 2 , ..., d † r }, from which we get the non-orthonormalized row vectors Ψ, the corresponding orthonormalized row vectors Ψ will give the correct correlation matrix of the product state r i=1 d † i |0⟩.To see that, we make use of the QR decomposition of the transpose matrix ΨT k×r = Q k×k R k×r where Q is unitary and R is upper triangular in the first r rows and zero in the remaining k − r rows.So the state i d † i |0⟩ = i l mi,l c † l |0⟩ represented by Ψ is equivalent to i j≤i,l C AA , C AB , and C BB are the corresponding submatrices of C. Ω AA and U A are respectively the eigenvalues and diagonalizing unitary of the restricted correlation matrix C AA ; similarly for Ω BB and U B .Ω AB would generally depend on the arbitrary phases in U A

TABLE II .
Truncation error list of PEPS for 100-site SSH model.

TABLE III .
Truncation error list of PEPS for 2D OAI on a 20 × 20 square lattice.