## Constructive Fermionic Matrix Product States for Projected Fermi Sea

Kangle Li, 1, ∗ Yan-Bai Zhang, 1, ∗ and Hoi Chun Po 1, †

1 Department of Physics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong SAR, China

Projected wave functions offer a means for incorporating local correlation effects in gapless electronic phases of matter like metals. Although such wave functions can be readily specified formally, it is challenging to compute their associated physical observables. Tensor network approaches offer a modern numerical method for this task. In this work, we develop and demonstrate a constructive tensor-network approach for obtaining physical quantities, like fermion two-point functions and density-density correlation functions, for one-dimensional projected Fermi sea states. We benchmark our method against exact analytical results for spin1 2 / electrons subjected to the Gutzwiller projection, and then present results on spinless fermions with nearest-neighbor repulsion. For a state with two pairs of Fermi points, we reveal a correlation-tuning of the characteristic wave vector of the charge density modulation in the system.

## 1. INTRODUCTION

Projected wave functions, like the Gutzwiller ansatz [1, 2], provide a general framework for incorporating electron correlation effects into the description of gapless phases of matter [3-18]). Intuitively, such wave functions contain two key ingredients: a 'projection,' which embodies local correlation effects like the energy penalty of double electron occupancy on the same site, and a pre-projection parent state, like a filled Fermi sea, which captures the qualitative features of the lowenergy physical properties of interest. In strongly correlated systems, projected wave functions also provide a way for constructing fractionalized ground states by projecting out unphysical states in an enlarged parton Hilbert space [19-22].

Valuable insights into the interplay between local correlation effects and the characteristic features anticipated from a mean-field starting point can be obtained by investigating how the projection strength affects physical observables. Generally speaking, however, such physical observables are challenging to compute. While exact analytical results are available for certain special models [6, 7, 23-25], in general either approximate perturbative treatment or numerical Monte Carlo methods are needed for extracting the physical properties of a projected wave function [8-11].

In this work, we develop and demonstrate a constructive approach for the TN-aided evaluation of projected wave functions. We restrict our scope to projected Fermi sea states in one spatial dimension, which allows us to leverage the highly efficient matrix product state (MPS) and transfer matrix frameworks for evaluating physical observables. We report results on the fermion two-point functions, as well as the density-density and spin-spin correlation functions, of projected Fermi sea states. Our results are first benchmarked against the analytical results available in the standard Gutzwiller problem of a filled Fermi sea of spin1 / 2 electrons. Next, we apply our method to study spinless fermions subjected to nearest-neighbor repulsion, for which exact analytical results are not available. This allows us to explore the interplay between such local correlations and the fermiology of the mean-field parent state. When there are two pairs of distinct Fermi points, we reveal a competition between the parent-state fermiology and the local correlation effects in the charge density modulation. In contrast to expectations from the weak coupling limit, we find that a pronounced peak in the static structure factor develops at a wave vector which cannot be simply inferred from the Fermi wave vectors characterizing the Fermi surfaces.

Tensor network (TN) states provide an alternative starting point for describing strongly correlated phases of matter [2635]. More broadly, TN also serves as a natural framework for describing a projected wave function: once the TN form of the pre-projection parent state is obtained, the projection can be implemented readily by applying additional local operators to the physical legs of the tensors. Earlier works have demonstrated the promise of this approach in both gapped and gapless systems in which double occupancy of the electrons or partons are projected away [20-22, 36-39]. Aside from providing a new practical means to evaluate physical observables in Gutzwiller wave functions, such projected states could also be used for initializing variational TN calculations in which the competitiveness of a more entangled states [21, 36], like those featuring possible Fermi surfaces, could be explored.

∗ These authors contributed equally.

† hcpo@ust.hk

Our approach starts from constructing the MPS representation of a free-fermionic state at the level of correlation matrix [21, 34, 39-43], which encodes all the data present in the free-fermion state. Various technical approaches exist for this task, including the successive discovery and application of local rotations [34], as well as the variational optimization of ansatz local tensors [27, 29, 42]. Here, we construct fermionic Gaussian MPS based on Schmidt decomposition. This can be achieved either by sweeping across the whole finite chain directly [17], or by first decomposing the locally defined Wannier functions of the target state and then stacking them into the full MPS representation [44]. The local Gaussian tensors can then be turned into many-body ones [35], which enables the subsequent incorporation of Gutzwiller-like projections. In our approach, the fermionic nature of the tensors is maintained and the many-body bond spaces can all be interpreted as fermionic Fock spaces. This facilitates the translation of Gaussian-level tensor contractions to the many-body treatment.

Although the scope of this work is restricted to projected

Fermi sea states in one spatial dimension, we note that our general constructive approach for projected wave functions can be employed in any spatial dimensions, with the main difficulty stemming only form the usual increased complexity in evaluating physical observables in higher-dimensional TN states. Our work helps pave the way for exploring the power of TN methods in studying gapless phases of fermions in higher dimensions.

The paper is organized as follows. First we introduce the basic tools needed for construction and computation of fermionic MPS in Section 2. Then we show the explicit example of 1d Fermi sea state with two schemes in Section 3. Gutzwiller projected states are considered in Section 4, where we benchmark with some known analytical results and present our main results on the interplay between nearest-neighbor repulsion and the fermiology of the parent state.

## 2. BASIC TOOLS FOR FERMIONIC MPS

## i. Overview of the procedures

Our starting point is a free-fermion state corresponding to a filled Fermi sea state. All the data in this parent state are encapsulated in its (static) two-point functions. On a lattice, these two-point functions are typically viewed as entries of a matrix, called the covariance matrix [40, 45, 46] . If the state has a definite particle number, a more restricted form, defined below as the correlation matrix, would suffice for describing the state. Since we are interested in projected Fermi sea states, we assume the particle number of the state is conserved throughout. Correspondingly, we specialize our discussion to the case of correlation matrices, although the same operations can be readily generalized to covariance matrices. We also use 'free-fermion' and 'Gaussian' interchangeably, with the two understood to be the same when the particle number is conserved.

Next, we represent the parent state as a fermionic Gaussian MPS (fGMPS), for which all local tensors are represented by their correlation matrices[21, 34, 35, 43]. A key step in obtaining the desired MPS representation, as in the usual case, is the Schmidt decomposition of the state with respect to a bipartition of the system. For Gaussian states, the required Schmidt decompositions can be performed in terms of the correlation matrices (Section 2 ii). Conversely, the decomposed pieces can be recombined by performing suitable partial traces, which for Gaussian states can again be achieved through a simple contraction formula (Section 2 iii).

Through an MPS representation of the parent Gaussian state, local correlation effects can now be incorporated readily by applying suitable matrix product operators to the parent state. We will focus on the operators corresponding to soft projectors which either penalize double occupancy of spin1 / 2 fermions, or the simultaneous presence of spinless fermions on two adjacent sites. As such operators are interacting in nature, we need to first convert all the local tensors in the parent fGMPS into a many-body form. This can be achieved through the standard Jordan-Wigner (JW) trans- formation, which also allows one to compute the contraction of local fermionic tensors through suitably defined matrices (Section 2 iv). We can then extract physical observables from the MPS using transfer matrix techniques.

Throughout the paper, we denote ˆ c † / ˆ c as the creation/ annihilation operators of complex fermions.

## ii. Fermionic Gaussian matrix product state

Thanks to the Wick's theorem, all physical information of a free-fermion state is encoded in its correlation matrix C ij = ⟨ ˆ ˆ c c i † j ⟩ . We discuss in this subsection how the usual MPS techniques can be adapted for correlation matrices.

Schmidt Decomposition . Given a free-fermion state, the Schmidt decomposition of the state with respect to a subregion A and its complement B can be done using only the data in the correlation matrix C . Let us write C in a block-matrix form as

$$C = \left ( \begin{array} { c c } C _ { A } & C _ { A B } \\ C _ { A B } ^ { \dagger } & C _ { B } \end{array} \right ).$$

We focus on the restricted correlation matrix C A and diagonalize it with a unitary matrix, i.e. C A = U Λ U † . Then the eigenvalues λ i are related to the Schmidt weights: for λ i = 0 ( λ i = 1 ), the mode is fully filled (empty), contributing no entanglement across the bipartition[46, 47]. The rotated matrix

$$C ^ { \prime } = \left ( \begin{array} { c c } \Lambda & U ^ { \dagger } C _ { A B } \\ C _ { A B } U & C _ { B } \end{array} \right ) \quad \quad \ ( 2 )$$

is still a pure state, with the modes corresponding to the block occupied by Λ being regarded as bond modes. Since modes with λ i = 0 or 1 can only be exactly filled or empty in any Schmidt states with finite Schmidt weights, they can be stripped from the bond Hilbert space. In practice, we may set a small numerical threshold ε such that all modes with λ i &lt; ε or λ &gt; 1 -ε are treated as frozen. These frozen modes in C ′ can be stripped away by contracting with the suitable fully filled or empty states. Alternatively, one can also simply remove the rows and columns corresponding to these modes, and then incorporate an orthonormalization step to maintain the purity of the reduced state (see appendices of [44] for more details).

If the state C ′ is truncated, the unitary U becomes an isometry. For example, if a mode λ i ≈ 0 is frozen to filled, then we just remove the corresponding column vector in U ; but if a mode λ i ≈ 1 in C ′ is frozen to empty, then we need to keep it for purification, see below.

Kernel purification. The above Schmidt process returns a unitary U and a correlation matrix C ′ . We can obtain a fermionic Gaussian state corresponding to U by doing purification based on a kernel state C ker . This procedure can be viewed as a specialized form of the channel-state duality [48, 49], in which we ensure the purified state is also fermionic Gaussian. Generally, the kernel state can be chosen to be a maximally entangled state between the physical Hilbert space and an auxiliary

Hilbert space of the same dimension. Here, we use a freefermion kernel state with the correlation matrix

$$C ^ { k e r } = \frac { 1 } { 2 } \left ( \begin{array} { c c c c } \mathbb { I } & s \mathbb { I } \\ s ^ { * } \mathbb { I } & \mathbb { I } \end{array} \right ), \quad & ( 3 ) \quad & \text{dec} \\ \text{an arbitrary complex phase factor. In terms of cre-} \quad & \text{wh} \\ \quad \,. \quad \,. \quad \,. \quad \,. \quad \,. \quad \,.$$

✶ ✶ where s is an arbitrary complex phase factor. In terms of creation operators, the state is (up to a normalization factor)

$$| \psi ^ { k e r } ( s ) \rangle = \prod _ { i } ( - \hat { c } _ { i } ^ { \dagger } + s \hat { d } _ { i } ^ { \dagger } ) \, | 0 \rangle \,, \quad \ \ ( 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where we use ˆ c † i to denote the physical modes in U and ˆ d † i to denote the auxiliary modes introduced for purification. For simplicity, we choose s = 1 . Then the unitary is purified into a free-fermion state with the correlation matrix

$$C _ { | U \rangle } = \left ( \begin{array} { c } U & \\ & 1 \end{array} \right ) C ^ { k e r } \left ( \begin{array} { c } U ^ { \dagger } & \\ & 1 \end{array} \right ). \quad \ \ ( 5 )$$

✶ ✶ If U is an m × m matrix, then we get a free-fermion state with 2 m modes, the first m modes are physical modes in subregion A and the second set of m modes are assigned to the bonds.

As previously mentioned, the unitary matrix U may be truncated to an isometry when there are frozen modes in the Schmidt decomposition. Suppose the isometry, denoted by U ′ , is composed of n 1 dynamical modes, assembled into the columns of a submatrix U (1) , and another n 2 number of modes, corresponding to the submatrix U (2) , which have λ = 1 and are therefore frozen. Then, we can purify the isometry by applying to a kernel state together with some frozen modes

$$\tilde { C } ^ { k e r } = \left ( \begin{array} { c } C _ { 2 n _ { 1 } \times 2 n _ { 1 } } ^ { k e r } \\ & 1 _ { n _ { 2 } \times n _ { 2 } } \end{array} \right ), \quad \quad ( 7 )$$

$$\tilde { U } & = \left ( \begin{array} { c c c } U _ { m \times n _ { 1 } } ^ { ( 1 ) } & U _ { m \times n _ { 2 } } ^ { ( 2 ) } \\ & 1 _ { n _ { 1 } \times n _ { 1 } } \end{array} \right ) & \quad ( 6 ) & \text{ial} \, \text{is gi} \\ \text{kernel state together with some frozen modes} & \quad \tilde { C } _ { B } \end{array} \right )$$

✶ to get a ( m + n 1 ) -mode purified state represented by correlation matrix C | ˜ U ⟩ = ˜ ˜ UC ker ˜ U † .

The Schmidt decomposition above is canonical on one side; if A resides on the left of B , then the fermionic state | U ⟩ is left-canonical as a tensor. Alternatively, one can also perform a left-right-canonical Schmidt decomposition, with two unitary matrices U and V diagonalizing C A and C B . Then we obtain three pieces

$$n \text{ three pieces} \\ C _ { | U \rangle } & = \left ( \begin{array} { c c c c } 1 & U \\ U ^ { \dagger } & 1 \end{array} \right ), \quad & \text{For} \quad & \text{tracl} \\ C _ { \Lambda } & = \left ( \begin{array} { c c c c } \Lambda _ { 0 } & U ^ { \dagger } C _ { 0 1 } V \\ V ^ { \dagger } C _ { 0 1 } ^ { \dagger } U & \Lambda _ { 1 } \end{array} \right ), \quad & ( 8 ) & \text{ing} \\ C _ { | V \rangle } & = \left ( \begin{array} { c c c c } 1 & V ^ { \dagger } \\ V & 1 \end{array} \right ), \quad & \text{is e} \\ \text{dual} \\ \mathbb { e } C _ { | U \rangle } & \text{is left-canonical,} C _ { | V \rangle } & \text{is right-canonical, and} C _ { \Lambda } \\ \mathfrak { d } \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

✶ where C | U ⟩ is left-canonical, C | V ⟩ is right-canonical, and C Λ is a 'diagonal' state.

The Schmidt decomposition serves as a basic step for obtaining a MPS representation of a free-fermion state, as discussed in Sec. 3.

## iii. Partial Trace and Kernel Contraction

Having discussed how a fermionic state can be Schmidtdecomposed into local states, we next address how such local states can be 'contracted' back into a global one. The strategy, which applies equally well in the many-body case as is needed later, is to perform a partial trace in the tensor product Hilbert space.

Consider two density matrices ˆ ρ PB ∈ H P ⊗ H B and ˆ ρ ′ B ∈ H B . The partial trace of the two density matrices ˜ ρ P = Tr B [ ˆ ρ PB ˆ ρ ′ B ] is another density matrix defined in Hilbert space H P . The mentioned contraction can be rephrased in terms of such partial traces. Given two density matrices ˆ ρ AB ∈ H A ⊗ H f B and ˆ ρ B C ′ ∈ H B ′ ⊗ H f C , where the two sub-Hilbert spaces H B and H B ′ has the same dimension, we define the BB ′ contracted state to be

$$\hat { \rho } _ { A C } = \text{Tr} \left [ ( \hat { \rho } _ { A B } \otimes _ { f } \hat { \rho } _ { B ^ { \prime } C } ) \rho _ { B B ^ { \prime } } ^ { k e r } \right ] \text{ \quad \ } ( 9 )$$

where ˆ ρ ker BB ′ is the density matrix of a kernel state, which implicitly defines an identification of the fermion modes in Hilbert space B with those in B ′ . We can chose this kernel state to be a pure free-fermion state with correlation matrix taking the same form as in Eq. (4), with the complex phase ˜ s . For consistency, however, the phase factor has to be chosen to be ˜ = s -s ∗ [44], where s denotes the phase chosen for the purification kernel (we choose s = 1 and hence ˜ = s -1 ). We call this ˆ ρ ker BB ′ the contraction kernel .

Specializing to free-fermion states, the discussion above also returns a free-fermion state after contraction. The partial trace formula for correlation matrices [38, 40, 42, 44, 45] is given as follows: denote the two states as C for ˆ ρ AB and ˜ C B for ˆ ˜ ρ B with

$$C = \left ( \begin{array} { c c } C _ { A } & C _ { A B } \\ C _ { A B } ^ { \dagger } & C _ { B } \end{array} \right ),$$

then the contracted state C ′ is given by

$$\underset { \begin{array} { c } \text{ide} ; \\ \end{array} } { C ^ { \prime } } & = C _ { A } - 2 C _ { A B } \left [ 2 C _ { B } - 1 + ( 2 \tilde { C } _ { B } - 1 ) ^ { - 1 } \right ] ^ { - 1 } C _ { A B } ^ { \dagger }, \\ \text{mm} \\ \text{ini} - \quad } \text{and if } \tilde { C } _ { B } \text{ is a pure state, this can be further simplified by}$$

✶ For two Gaussian states C AB , C B C ′ with B,B ′ to be contracted, our kernel contraction formula is obtained by replacing the C by C AB ⊕ C B C ′ , and replacing ˜ C B by C ker . We also remark that the kernel contraction of ˆ ρ AB with a state ˆ ρ B ′ is equivalent to the partial trace of ˆ ρ AB and the particle-hole dual of ˆ ρ B ′ , which has correlation matrix ✶ -C B ′ .

and if ˜ C B is a pure state, this can be further simplified by (2 ˜ C B -✶ ) 2 = ✶ , C ′ = C A -C AB [ C B + ˜ C B -] -1 C † AB . (11)

## iv. Many-Body fermionic matrix product states

To move on to correlated states, we promote the fGMPS to a many-body representation. Given a set of pure Gaussian states

as local tensors of fMPS, we can get many-body local tensors by writing down the wave functions of each Gaussian states. Starting from the correlation matrix, this can be achieved by taking the Slater determinant of all the filled modes. Equivalently, in second-quantized form we can write

$$| \Psi \rangle = \hat { c } _ { f _ { 0 } } ^ { \dagger } \hat { c } _ { f _ { 1 } } ^ { \dagger } \cdots \hat { c } _ { f _ { \nu - 1 } } ^ { \dagger } | \text{vac} \rangle, \text{ \quad \ \ } ( 1 2 )$$

where | vac ⟩ denotes the fermionic vacuum state, ν denotes the charge of the state, and ˆ c † f i for i = 0 1 , , . . . , ν -1 are the filled modes. For numerical calculations, we may represent all the locally defined fermionic operators using a JW transformation [50], and with that | Ψ ⟩ can be represented as a 2 N dimensional vector, where N ≥ ν denotes the number of fermion modes in the local Hilbert space.

As previously discussed, the kernel contraction of local tensors is defined through the partial trace. Since all local tensors are viewed as pure fermionic states, this partial trace process can be done in terms of the wave functions directly. Suppose the wave functions of two local states are represented numerically by two vectors Φ AB and Ψ B C ′ in occupation number basis. The contracted state is given by

$$( [ \hat { P } _ { A } ] ) ^ { N _ { B } } \Phi _ { A B } [ K ] _ { B B ^ { \prime } } \Psi _ { B ^ { \prime } C }, \quad \ \ ( 1 3 ) \quad \text{or} \, \ u$$

where [ ˆ ] P C is the matrix form of a parity operator of the modes in sub-Hilbert space H A and its power N B is the number of modes in H B . This extra parity operator is due to the JW representation and does not effect if N B is even (see Appendix E for a derivation).

The kernel purification can be extended to the many-body setting. Given a many-body operator ˆ O , the kernel purification is | ˆ O ⟩ = ( ˆ O ⊗ ✶ ) | ψ ker ⟩ . This can be reconciled with the thermofield double state if ˆ O is the square root of a density matrix. This purification preserves fermionic structure naturally, and the operation of ˆ O to a quantum state | φ ⟩ is equal to the contraction of purified | ˆ O ⟩ and φ , up to an overall sign which depends on the parity of the state (Appendix E). This allows us to treat fermionic matrix product operators in terms of local fermionic states as well.

Computing Observables . Transfer matrix is used to compute observables. After tracing out the physical modes of a local tensor and regrouping elements, we are left with a transfer matrix M i whose rows/ columns are corresponding to the left/ right bonds. Combining with contraction kernels they can be multiplied together, and the expectation value is the trace of the product

$$\langle \hat { O } _ { i } \rangle = \frac { \text{Tr} \left [ M _ { 0 } ^ { ( K ) } M _ { 1 } ^ { ( K ) } \dots M _ { \hat { O } _ { i } } ^ { ( K ) } M _ { i + 1 } ^ { ( K ) } \dots M _ { L - 1 } ^ { ( K ) } \Xi ^ { N ^ { \pm } + 1 } \right ] } { \text{Tr} \left [ \prod _ { i } M _ { i } ^ { ( K ) } \Xi ^ { N ^ { \pm } + 1 } \right ] }, \quad \text{Giv} \\ \text{where } M ^ { ( K ) } \text{ is a ```comnosite`` transfer matrix with the suner-} \quad \text{As an}$$

where M ( K ) is a 'composite' transfer matrix with the superscript ( K ) denoting that it has been multiplied with the kernel matrix. M ( K ) ˆ O i is the transfer matrix with an insertion of operator ˆ O i in the physical legs. One may worry about potential

̂

̂

̂

̂

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-5b9acbb1-8e2f-4cce-9de5-449f10dc0f79.jpg)

̂

̂

Figure 1. Simplifying the insertion of Jordan-Wigner string

sign problems since, for an odd number of bond modes, Eq. (13) has an extra parity matrix. In fact this does not influence the expectation values for the bonds inside the chain, and only possible signs will happen at the two ends of the chain if the fMPS chain is closed. In Eq. (14), we denote the parity matrix of boundary bond modes as Ξ , and the exponent N t is the half of the total number of bond modes in the MPS chain, or equivalently, the total occupation number of all kernel states used for contractions. The derivation of Eq.(14) is in Appendix E.

For the computation of non-local correlation functions, we need to insert suitable JW strings. For instance, to compute the two-point function ⟨ ˆ ˆ c c † i i + n ⟩ , the usual way with onedimensional MPS is to insert a JW string from the site i to the site i + n -1 . In our fMPS representation, since all local tensors are particle-number conserved fermion states, this JW string can be simplified. The expectation value is

$$\langle \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + n } \rangle & = \text{Tr} _ { \text{phys} } \left [ \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + n } \text{Tr} _ { \text{bond} } ( \hat { \rho } _ { | \psi \rangle } \hat { \rho } _ { K } ) \right ] \\ & = \text{Tr} _ { \text{bond} } \text{Tr} _ { \text{phys} } \left [ \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + n } \hat { \rho } _ { | \psi \rangle } \hat { \rho } _ { K } \right ],$$

where we use ˆ ρ ψ to represent the density operator of the product state of all local tensors and ˆ ρ K the product state of all contraction kernels. Then the JW string is going through all physical modes and bond modes from the mode i to the mode i + n , as shown in Fig. 1. But all local tensors are particle-number conserved, so the action of JW-string, which is a string of local parity operators, on a local tensor simply returns an overall sign; see Fig. 1.

## 3. CONSTRUCTION OF GAUSSIAN FMPS

Given a fermionic Gaussian state, we want to construct an MPS representation from the correlation matrix[34, 35, 37]. As an illustration, consider a single band Hamiltonian of spinless fermions

$$\text{er-} _ { \text{ial} } \quad \hat { H } = - t \sum _ { i } \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + 1 } + h. c. = - \sum _ { k } ( 2 t \cos k ) \hat { c } _ { k } ^ { \dagger } \hat { c } _ { k }, \quad ( 1 6 )$$

̂

̂

̂

which has the ground state

$$\left | \psi \right \rangle = \prod _ { \left | k \right | \leq \pi / 2 } \hat { c } _ { k } ^ { \dagger } \left | 0 \right \rangle. \quad \quad \left ( 1 7 \right )$$

This ground state is fully described by its correlation matrix

$$C _ { i, j } = \langle \psi | \, \hat { c } _ { i } \hat { c } _ { j } ^ { \dagger } \, | \psi \rangle = \Psi \Psi ^ { \dagger } \, \text{ \quad \ \ } ( 1 8 )$$

where Ψ contains the column eigenvectors of the singleparticle Hamiltonian corresponding to all the filled modes.

We discuss two schemes for obtaining the MPS representation of the state. As usual, one can obtain an MPS representation by successively performing Schmdit decomposition along the chain [21, 42, 43] from, say, left to right. We will refer to this as 'Scheme I.' Alternatively, we can also construct a translation-invariant MPS representation directly, which we refer to as 'Scheme II.' As we will discuss in details below, Scheme I can faithfully capture the Fermi surface singularity of the original filled Fermi sea, but due to the lack of translation invariance our calculations are limited to relatively small system sizes. In contrast, Scheme II enables the evaluation of physical observables on much longer chains, but, due to the finite bond dimension, the resulting translation-invariant MPS has a smeared Fermi surface although the Fermi wave vector can still be readily recognized [51, 52].

## i. Scheme I: Successive Schmidt Decomposition

The first scheme is to successively perform the Schmidt decomposition discussed in Sec. 2 ii. This can be summarized in

Table I. Summary of performance in obtaining matrix product representation of a filled Fermi sea in Scheme I. The second column is the entanglement threshold for truncation, and the third column is the number of maximal bond modes.The fidelity is |⟨ ψ fMPS | ψ exact ⟩| . The last two columns are energy densities (in units of t ) computed from fMPS results and exact ground states.

|   L c | ent th   |   m.b.m |   fidelity |   E M /L |   E gs /L |
|-------|----------|---------|------------|----------|-----------|
|     1 | 10 - 2   |       4 |      0.58  |  -0.63   |   -0.6365 |
|     2 | 10 - 2   |       4 |      0.662 |  -0.6317 |   -0.6365 |
|     1 | 10 - 3   |       7 |      0.953 |  -0.6359 |   -0.6365 |
|     2 | 10 - 3   |       7 |      0.974 |  -0.6362 |   -0.6365 |

The advantage of this first scheme is that it is more straightforward, and is also more accurate in describing the lowenergy behavior. For instance, the momentum distribution shows a sharp Fermi surface with a discontinuity. The shortcoming is that the obtained tensors are not translationinvariant, and it requires significantly more computational resources for large the system size in the subsequent many-

- 1. Choose a subregion of length L c . Define a local Hilbert space H x by combining all the existing bond modes with the physical modes in the sub region. Diagonalize the corresponding reduced correlation matrix ( C ) | i,j ∈H x with unitary matrix U . Denoting the eigenvalues as Λ , we have ( C ) | i,j ∈H x = U Λ U † .
- 2. Freeze the eigenmodes in H L with eigenvalues λ &lt; ε or λ &gt; 1 -ε , where ε is a chosen entanglement threshold.
- 3. Use the kernel state to purify the unitary (or the isometry if some modes are frozen) U to get a Gaussian state | M x ⟩ .
- 4. Take the transformed correlation matrix as the new state with new bond modes defined from the Schdmit decomposition above. Repeat from step 1 until we have swept through the entire chain.

By doing these steps successively, we can obtain a set of local fermionic Gaussian tensors { M ,M ,..., M 0 1 ℓ } , ℓ = L/L c . These local fermionic Gaussian tensors include left bonds modes, physical modes and right bond modes. Contracting all bond modes we can obtain the whole fGMPS state. Table I summarizes some numerical results on the MPS approximation of the exact state ψ . The momentum distribution n k = ⟨ ˆ ˆ c † k c ⟩ computed from the correlation matrix elements is shown in Fig. 2(a).

body calculation after performing projection. Nevertheless, this construction provides some information about the correlation functions in a finite-size system, and can be regraded as a crosscheck before computing a larger system through Scheme II below.

## ii. Scheme II: Translation-invariant tensors

Alternatively, one could also construct translation-invariant MPS for the filled Fermi sea. This can be achieved by optimizing an MPS defined by local tensors with fixed, predefined bond dimensions to minimize the energy of a model Hamiltonian [37-39, 52]. Here, we adopt a more constructive approach in which the set of translation-invariant local tensors are obtained using the method introduced in Ref. [44], which is dubbed 'tree-stack-compress.'

This approach returns a local tensor defined on a supercell of two physical sites, i.e., a local state defined in a Hilbert space of two physical modes together with some bond modes defined on both the left and right sides. The resulting MPS state is then obtained by repeating this local tensor, and so it has discrete translation invariance in units of two sites. This

Figure 2. (a) The momentum distribution of the fermionic Gaussian matrix product state (fGMPS) constructed from Schmidt decomposition (Scheme I). The entanglement threshold is 10 -3 and the chain is of length 128. (b) The momentum distribution of fGMPS constructed from tree-stack method. The number of bond modes is 14, and the length of the fGMPS is 500. (c) The crosses are the correlation functions of the translation-invariant fGMPS and the dashed lines are exact correlation functions. (d) The reconstructed dispersion for modes at each momentum from translation-invariant tensors (Scheme II), in the folded Brillouin zone. The zoomed in region shows the quadratic dispersion near the Fermi point.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-eb6075c8-9d50-4cf0-9e28-05a80f66db24.jpg)

subsequently saves computational cost when computing physical observables.

We first treat the 1d spinless fermion chain using two-site unit cells. The original Brillouin zone is folded, and we get a two-band dispersion. The lower band is fully filled, which enables a Wannierization of the lower band. The treestack-compress method (Fig. 3) proceeds by first finding a fGMPS representation for each of the Wannier mode, and then stacking all their translation-related copies to obtain the full fGMPS. The bond dimension of the resulting fGMPS can then be brought down through a compression step. Applied to the problem of a filled Fermi sea state, however, the Wannier modes are only algebraically instead of exponentially decaying. Nevertheless, one can still truncate the Wannier mode down to a finite region of some radius about the Wannier center, which amounts to approximating the filled Fermi sea by an atomic insulator. Such approximation works well for the highenergy (or short-distance) properties of the state, but is inherently incapable of capturing the Fermi surface singularity. Nevertheless, similar to existing works on the study of critical states using translation-invariant MPS [39, 52], a smeared Fermi surface is still retained and its sharpness is generically controlled by the bond dimension of the tensor-network repre- sentation. In our tree-stack-compress method, this is in turned controlled by the truncation radius used in obtaining the approximate Wannier functions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-cb83c002-03e0-424f-a1c2-0b61f7f5faf4.jpg)

. . . . . .

Figure 3. Schematic for the tree-stack construction of a free-fermion matrix product state. Each row represents a tensor network for a Wannier state with truncation radius 2. The colored region shows how the tensors are stacked to give a translation-invariant MPS approximation to the full state.

We summarize some numerical results here. For example, given the truncation radius of the Wannier function to be 16 ,

we can obtain a local tensor with 14 bond modes on each side. The MPS representation is then further compressed such that each bond contains only four fermion modes. For the compressed tensor, the difference of energy density from the exact ground state is 5 × 10 -4 . Its associated momentum distribution is shown in Fig. 2(b), which shows a smeared Fermi surface discontinuity about k F = π/ 2 . A comparison between the exact and MPS-approximated two-point functions is shown in Fig. 2(c), which demonstrates good agreement in the short-distance. The momentum-resolved energy spectrum for the fGMPS state is represented in Fig. 2(d), which again shows good agreement with the exact result in the high-energy part of the spectrum, i.e., away from the Fermi surface. However, as shown in the inset, there is an inherent rounding of the linear dispersion about the Fermi surface, corresponding to its smearing.

The Fermi surface smearing is a result of the critical nature of the filled Fermi sea state, which leads to a logarithmic contribution to the entanglement area law [51, 53, 54]. As a result, the approximation of this critical state using an MPS with a finite bond dimension leads to an unavoidable failure in the long-wavelength limit. Nevertheless, the degree of smearing, i.e., the sharpness of the Fermi surface, can generally be controlled through the bond dimension [52, 55]. In Appendix F, we demonstrate how the smearing width can be systematically reduced at the cost of an increasing bond dimension.

## 4. GUTZWILLER PROJECTED WAVE FUNCTIONS

Having obtained the MPS approximation of the filled Fermi sea state, we can move on to construction of many-bodylocal tensors for interacting problems. One typical example is the Gutzwiller projected wave function [1]. The standard Gutzwiller projector is

$$\hat { G } _ { \text{double} } = \prod _ { r } ( 1 - ( 1 - g ) \hat { n } _ { r, \uparrow }, \hat { n } _ { r, \downarrow } ), \quad ( 1 9 ) \quad \text{and } t$$

which for g &lt; 1 penalizes configurations with double occupancy on any site. Here, we consider repulsive interactions between the fermions and so take g ∈ [0 , 1] . The (unnormalized) Gutzwiller projected wave functions are defined as ˆ G ψ | ⟩ where | ψ ⟩ is the filled Fermi sea state. In one spatial dimension, physical observables of the standard Gutzwiller projected states can be evaluated analytically for a finite g when | ψ ⟩ is half-filled [23-25]. This case then serves as a benchmark for our numerical methods.

We would also consider another Gutzwiller-type projector, suitable for spinless fermion, which penalizes configurations based on nearest-neighbor (nnbr) repulsion:

$$\hat { G } _ { \text{nnbr} } = \prod _ { x } ( 1 - ( 1 - g ) \hat { n } _ { x } \hat { n } _ { x + 1 } ). \quad \quad ( 2 0 )$$

Notice that, as ˆ G double acts on different sites independently, its action on a MPS can be computed readily. In contrast, ˆ G nnbr imprints inter-site correlations on the state and does

Figure 4. (a) The momentum distribution of the ˆ G double -projected state of 40 sites. The entanglement threshold of the fGMPS is 10 -2 and the maximal number of bond modes in the chain is 3. (b) The momentum distribution for ˆ G nnbr -projected state in a chain of 64 sites

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-f9139a9f-5d41-4de1-8d1e-5c3364426d2f.jpg)

not act in a strictly local manner. Nevertheless, one can readily derive its matrix product representation, as is shown in Appendix D. We note that earlier works have employed the Gutzwiller approximation to study this problem of spinless fermions projected by ˆ G nnbr [56-59]. Here, we show that the tensor networks can provide an alternative numerical means for extracting physical observables like momentum distributions and static structure factors from the projected wave function.

## i. Many-body MPS from Scheme I

We begin by studying the projected Fermi sea states obtained by applying the Gutzwiller-type projectors to the fGMPS obtained through successive Schmidt decomposition (c.f. Sec. 3 i). We compute various correlation functions for these projected states. First, we benchmark our results in

the spinful Gutzwiller projector ˆ G double acting on the halffilled Fermi sea. The two-point functions ⟨ ˆ c † i, ↑ ˆ c j, ↑ ⟩ of a spinful fermion chain of L = 40 sites for varying g were computed. The momentum distribution n k , obtained through Fourier transform of the two-point functions, is shown in Fig. 4(a). In agreement with the exact results, we find that the Fermi surface discontinuity survives for finite g , although the size of the jump, corresponding to the quasi-particle weight, is reduced. In addition, the Gutzwiller projector modifies the momentum distribution far away from the Fermi surface; as we will see, such short-distance properties of the projected state can be well-reproduced using tensor-network methods. Notice that, due to the exact degeneracy at the Fermi surface, we choose the occupation number per spin species to be L 2 +1 . In contrast, the analytical result is for exact half-filling. The difference in fermion count, which is a finite-size effect, leads to a more significant deviation between the computed and the exact analytical result as the projection strength increases (i.e., for smaller values of g ).

Likewise, in Fig. 4(b) we present the momentum distribution of nearest-neighbor Gutzwiller projected state ˆ G nnbr | ψ ⟩ of a 64 -site chain. Similar to the spinful case, the Fermi surface persists for finite values of g , with the jump size is again renormalized to a smaller value than one. In the hardprojection limit of g = 0 , the projected state is exactly a cat state 1 √ 2 ( 0101 | ... ⟩ ± | 1010 ... ⟩ ) if the system is exactly half-

Another natural quantity to consider for the Gutzwiller projected wave function is the average double occupancy ⟨ ˆ n n ↑ ˆ ↓ ⟩ average = 1 L ∑ L l =1 ⟨ ˆ n l, ↑ ˆ n l, ↓ ⟩ , which is suppressed from the original filled Fermi sea. Fig. 5(c) shows the exact average double occupancy in [24] compared with our numerical results for g from 0 1 . to 1 with step size of 0 1 . , which accurately approached the exact results. Other physical correlation functions, such as spin-spin correlation and densitydensity correlation,

$$C ^ { S S } ( r ) & = \frac { 1 } { L } \sum _ { i } ( \langle \hat { S } _ { i } ^ { z } \hat { S } _ { i + r } ^ { z } \rangle - \langle \hat { S } _ { i } ^ { z } \rangle \langle \hat { S } _ { i + r } ^ { z } \rangle ) \quad \text{(21)} \quad \text{can I} \\ C ^ { N N } ( r ) & = \frac { 1 } { L } \sum _ { i } ( \langle \hat { n } _ { i } \hat { n } _ { i + r } \rangle - \langle \hat { n } _ { i } \rangle \langle \hat { n } _ { i + r } \rangle ) \quad \text{(21)} \quad \text{see F} \\ & \text{have}$$

̸

can also be computed. After Fourier transformation, we obtain their representation in momentum space C SS ( q ) and C NN ( q ) , known as the static structure factor for spin and charge (see, for example, [60]). In Fig. 5(b), the fMPS and exact spin-spin correlation in momentum space for g = 0 1 . , 0 3 . , 0 5 . , 1 are computed. When the momentum | k | &lt; π = 2 k F , the fMPS results accurately fit the exact results with all computed g . When | k | approaches π , the exact spinspin correlation for g = 1 becomes a kink, which is the 2 k F singularity. The fMPS results quantitatively are still close to the exact data but the cusps at q = π are rounded. Fig. 5(e) is the comparison between fMPS and exact density-density correlations. Unlike the case for spin-spin correlations, when filled. Correspondingly, the momentum distribution will be exactly 1 2 over the entire Brillouin zone. Our results are qualitatively similar to [57] where the author studied the n k by Gutzwiller approximation. We also included the case of exact half-filling in Fig. 5(f) for comparison with the subsequent calculations using Scheme II (see Sec. 4 ii).

## ii. Many-body MPS from Scheme II

While faithfully reproducing the Fermi surface discontinuity in the momentum distribution, the lack of translation invariance in the resulting tensor from the successive Schmidt decomposition (Scheme I) limits the system sizes we can reach. As an alternative, we now consider the projected Fermi sea state obtained from the translation-invariant local tensors described under Scheme II. First, we again benchmark our results against the standard Gutzwiller-projected state ˆ G double , where we use a local tensor with 4 bond modes on each side. The momentum distribution is shown in Fig. 5(a). We can compare it with the exact result as shown Fig. 5(d). We see they match well in the high-energy region far from the Fermi points. Near the k F = π/ 2 , the fMPS results have smeared Fermi surfaces, while the exact results has a sharp Fermi surface.

̸

g = 1 the fMPS results accurately resolve the density-density correlations for the whole momentum space.

After benchmarking against the exact analytic results in the standard one-dimensional Gutzwiller projected Fermi sea, we move on to consider other projected Fermi sea states. The momentum distribution of the nearest-neighbor Gutzwiller projector ˆ G nnbr | ψ ⟩ for spinless fermions is presented in Fig. 5. Up to the mentioned smearing of the Fermi surface singularity, the results are quantitatively consistent with the computations done with Scheme I, except that a much longer chain can now be assessed. The density-density correlation function is computed for strong projections, which shows a peak at 2 k F = π , which indicates the symmetry breaking at g = 0 ; see Fig. 5(f).

A similar but more interesting problem is the projected wave function of a ground state with two pairs of distinct Fermi points. As an example, we take the | ψ ⟩ to be the ground state of the following next-nearest-neighbour hopping Hamiltonian,

$$\hat { H } = - \sum _ { i = 1 } ^ { L } ( t _ { 1 } \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + 1 } + t _ { 2 } \hat { c } _ { i } ^ { \dagger } \hat { c } _ { i + 2 } + h. c. ) \quad ( 2 2 )$$

with t 2 = 5 t 1 . The solutions and correlation functions of this model are summarized in Appendix C. The dispersion relation of this Hamiltonian is ε k = -( t 1 cos k + t 2 cos 2 k ) , as shown in Fig. 6(a). By adjusting the chemical potential we get a ground state exactly at half-filling. For the chosen parameters here, the two distinct Fermi points in [0 , π ]

Figure 5. (a) The momentum distribution for spinful fermions with finite value of Gutzwiller projection g , based on spinless local tensor with 4 bond modes on each side, and its comparison with the analytical results. The system size is L = 1000 . (b) The average double occupancy compared with the analytical results. (c,d) The density-density correlation and spin-spin correlation in momentum space compared with the analytical results. (e) The momentum distribution of nearest-neighbor projected state. The crosses are computed with Scheme I (64 spinless fermions) and dots are computed from Scheme II (1000 spinless fermions). Both are exactly at half filling. (f) The density-density correlation function of nearest-neighbor projected state

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-a0ec0592-6465-49f6-a5ad-eb7c80d8bbdc.jpg)

are k F;1 ≈ 0 273 . π, k F;2 ≈ 0 773 . π . We use the method of Scheme II to construct the translation-invariant local tensor, which has four physical modes and four bond modes for each side. With these Gaussian tensors, we compute the correlations for ˆ G nnbr | ψ ⟩ . In Fig. 6(b) the momentum distribution is shown for the projected wave functions with different g . Similar to the cases of a single Fermi sea, these projected states with two Fermi seas also maintains smeared Fermi surface discontinuities.

We also study the density-density correlation function. For different values of t 2 /t 1 , the free fermion state before projection shows different features, see Fig. 8(b) in Appendix C. For different g , the computation of density-density correlation functions is shown in Fig. 6(c). At the Gaussian limit ( g = 1 ), the two peaks are at π/ 2 . As g decreases, corresponding to increasing projection and hence correlation strength, the two peaks sharpen and develop in magnitude. Their positions also shift towards q = π . Such correlation-driven enhancement in the static structure factor can be interpreted as a precursor for a charge-density wave instability. Interestingly, we find that the characteristic wavelength of this charge density modulation is continuously renormalized by the correlation strength, although the location of the Fermi points, as shown in Fig. 6(b), remains unchanged from the parent free-fermion state.

The behavior of the system in the limit of vanishing g warrants further discussion. In fact, when g = 0 the exact pro- jected state vanishes, as the input free-fermion state does not have any weights in the two states | 1010 ... ⟩ and | 0101 ... ⟩ (Appendix B). However, in our fGMPS representation of the input state, small weights in these two components are generated due to the approximations involved. As g vanishes, the relative weights of these unphysical components are generally amplified, which renders these results unreliable. For completeness, we show the behavior of C NN ( q ) for g &lt; 0 01 . in Appendix B, but for the mentioned reasons these results in the g → 0 limit should be interpreted with caution.

## 5. DISCUSSION

In this work, we formulate a constructive approach for evaluating physical observables in one-dimensional projected Fermi sea states using matrix product states. As a benchmark, we show that our results are in good agreement with the exact analytical results available for half-filled spin1 / 2 electrons subjected to the Gutzwiller projector which penalizes double occupancy. Notably, when translation-invariant states are used, the Fermi surface singularity is generally smeared out, as expected due to the short-range entangled nature of the matrix product ansatz. Nevertheless, the remnant of the Fermi surface is still discernible in physical observables like momentum distribution function and spin-spin or density-density cor-

Figure 6. (a) The dispersion relation of Eq. (22) for t 2 = 5 t 1 . (b) The momentum distribution of ˆ G nnbr projected wave function of 240 spinless fermions. (c) The density-density correlation function for ˆ G nnbr projected spinless fermions with two Fermi surfaces, for larger g and smaller values.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-015a4fac-19a7-4060-9d2b-0567bebed724.jpg)

relation functions. We further apply our approach to study the projected Fermi sea states of spinless fermions at half filling and subjected to nearest-neighbor repulsion. For the case of a single pair of Fermi points, the density-density correlation shows a peak at q = π as the projection becomes stronger, whereas for the case of two pairs of Fermi points, the two peaks originating from the pre-projection parent state develop in magnitude and shift towards q = π .

While we focus on the evaluation of physical observables in a pre-defined projected Fermi sea state in this work, our results can serve as the starting point for further variational calculations. Conceptually, the Gutzwiller-type projectors we considered can be readily generalized to any other short-range entangled operators admitting a matrix product representation.

The projected states we constructed, which retain smeared Fermi surfaces, can be taken as the initial state of wellestablished variational energy minimizers, like the densitymatrix renormalization group approach [21, 36, 43]. We note that in TN methods like the density-matrix renormalization group, the energy of a Hamiltonian is minimized over the variational TN ansatz. For one-dimensional fermions subjected to local interactions, one generically expects the non-interacting ground state, with a sharp Fermi surface, to be superseded by a Luttinger liquid with only a kink singularity instead of a discontinuity in the the momentum distribution function. In contrast, our approach is constructive in nature and corresponds to a TN representation of the projected Fermi sea states. As a result, the Fermi surface discontinuity is approximately retained, in the sense that a sharp transition of momentum occupation at the Fermi surface, with the sharpness controlled by the bond dimension of the TN, is reproduced in our method.

We also note that our approach can be applied to study gapless phases of fermions in higher dimensions, possibly in the context of a projected parton state [19-22, 36, 43]. For critical states with a vanishing Luttinger volume, our approach [44] for constructing a tensor-network representation of a gapless Gaussian state proceeds similarly in high dimensions; for more general states with a Fermi sea, however, our current approach is limited to commensurate cases for which the Luttinger volume vanishes under suitable Brillouin zone folding. Nevertheless, it is an interesting question to study the efficacy of the outlined approach in describing highly entangled quantum phases of matter, say topologically ordered system with spinon Fermi surfaces. Along another line, one may also combine the approach with entanglement renormalization [31, 61-63], which can provide a more faithful tensor-network representation of gapless states with sharp Fermi surfaces.

Note added. In finalizing the manuscript, we became aware of the work [64], which also discussed how to construct fermionic Gaussian MPS from correlation matrices with mode truncation.

## ACKNOWLEDGMENTS

This work is supported by the Hong Kong Research Grant Council (ECS 26308021) and the Croucher Foundation (CIA23SC01). We thank Hong-Hao Tu for helpful discussions and comments.

## Appendix A: Exact solutions of 1d Gutzwiller projected state

In the main part, we compared the numerical results with exact solutions for different correlation functions. Here we list the related formulas of the exact solutions based on [2325]. The exact momentum distribution formula used in Fig.

5(d) is:

$$n _ { k } ( g ) & = \frac { g ^ { 2 } + 4 g + 1 } { 2 ( 1 + g ) ^ { 2 } } + \frac { g ^ { 2 } } { ( 1 + g ^ { 2 } ) ^ { 2 } } \times \\ & \times \frac { 4 } { \pi \sqrt { ( 2 - G ) ^ { 2 } - ( \tilde { k } G ) ^ { 2 } } } \mathcal { K } ( \frac { G \sqrt { 1 - \tilde { k } ^ { 2 } } } { \sqrt { ( 2 - G ) ^ { 2 } - ( \tilde { k } G ) ^ { 2 } } } ) } \quad ( A 1 ) \quad \text{The\suprct} \\ & \quad \text{here}$$

where ˜ = 2 k | k /π | ≤ 2 , G = 1 -g 2 and

$$\mathcal { K } ( k ) = \int _ { 0 } ^ { 2 \pi } d \psi [ 1 - k ^ { 2 } \sin ^ { 2 } ( \psi ) ] ^ { - 1 / 2 } \quad ( A 2 ) \quad \tilde { \ a c o l }$$

is the first kind of complete elliptic integral.

The following is the formula for the spin-spin correlation function in real space C SS ( r ) = 1 L ∑ L l =1 ( ⟨ ˆ S z r + l ˆ S z l ⟩ G -⟨ ˆ S z r + l ⟩ G ⟨ ˆ S z l ⟩ G ) for electron density n = n ↑ + n ↓ = 1 :

$$C _ { j \geq 1 } ^ { S S } ( n = 1 ) = - \frac { 1 } { \pi j } \int _ { 0 } ^ { 1 } d y \frac { \sin ( \pi j y ) } { F ( y ) } \quad \text{ (A3)} \quad \text{ is m}$$

where the function F x ( ) ≡ 1 -(1 -g 2 ) x . The formula comes from the Fourier transform of the spin-spin correlation in momentum space:

$$C ^ { S S } ( k ) = \begin{cases} - \frac { 1 } { 1 - g ^ { 2 } } \ln F ( \frac { | k | } { \pi } ), 0 \leq | k | \leq 2 k _ { \text{F} } \quad & \text{name} \colon \\ - \frac { 1 } { 1 - g ^ { 2 } } \ln F ( n ), 2 k _ { \text{F} } \leq | k | \leq \pi \quad & \text{name} \colon \\ \text{ping t} \end{cases} \quad ( A 4 ) \quad \text{name} \colon$$

The density-density correlation in real space is defined as C NN ( r ) = 1 L ∑ L l =1 ( ⟨ ˆ n z r + l ˆ n z l ⟩ G - ⟨ ˆ n z r + l ⟩ G ⟨ ˆ n z l ⟩ G ) . Upon Fourier transform to the momentum space, one finds

$$C ^ { N N } ( k ; g, n = 1 ) = \frac { g ^ { 2 } } { 1 - g ^ { 2 } } \ln F ( - \frac { Q } { g ^ { 2 } } ) \quad ( \mathbf A 5 ) \quad \text{Frm} _ { / 1 0 1 }$$

the exact correlation in real space is:

$$C _ { j } ^ { N N } = & - \frac { \bar { p } } { ( \pi j ) ^ { 2 } } \{ \sin \left ( \pi j \right ) \ln g ^ { 2 } + \sin \bar { p } [ C i ( \bar { p } ) - C i ( p ) ] + \\ & + \cos \bar { p } [ S i ( p ) - S i ( \bar { p } ) ] \}$$

(A6)

where p = πj/ (1 -g 2 ) , ¯ p = g 2 p = πjg / 2 (1 -g 2 ) , Si x ( ) = ∫ x 0 sin t t dt is the sine integral function and C i ( x ) = -∫ i x nfty cos t t dt is the cosine integral function.

The formula for double occupancy formula in real space is:

$$C ^ { D } & = \frac { 1 } { L } \sum _ { l = 1 } ^ { L } \langle \hat { n } _ { l, \uparrow } \hat { n } _ { l, \downarrow } \rangle _ { G } = \sum _ { \substack { ( A 7 ) \quad \longrightarrow \\ \text{overla} \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } } \sum _ { \substack { ( A 9 ) \quad \longrightarrow \\ \text{On} \text{e} } }$$

$$\frac { g ^ { 2 } } { 2 ( 1 - g ^ { 2 } ) ^ { 2 } } [ - ( 1 - g ^ { 2 } ) ( n - m ) + \ln \frac { 1 - ( 1 - g ^ { 2 } ) m } { 1 - ( 1 - g ^ { 2 } ) n } ] \ ( A 8 ) \quad \text{contai} \\ \text{are ar}$$

where m = ( n ↑ -n ↓ ) is the magnetization and n is the electron density.

## Appendix B: ˆ G nnbr -Projected state at g = 0

The nearest-neighbor Gutzwiller projector ˆ G nnbr tends to suppress the weight of configurations with doubly occupied nearest-neighbor sites in the state. At the limit g = 0 , only configurations like | 0101 ... ⟩ and | 1010 ... ⟩ will survive. Such single configuration has exponentially small overlap with the pre-projected state, which we can show by a simple estimation as follows. Given the fermionic Gaussian state | ψ ⟩ = ∏ k ˆ c † k | 0 ⟩ is at half-filling, the component | 0101 ... ⟩ will have a coefficient determined by the Slater determinant

$$\langle 0 1 0 1 \dots | \psi \rangle = \frac { 1 } { L ^ { L / 4 } } \sum _ { \sigma \in \mathcal { S } _ { L / 2 } } \text{sgn} ( \sigma ) e ^ { i \sum _ { j } x _ { j } \sigma ( k _ { j } ) }, \quad ( \mathbf B 1 )$$

where σ ∈ S L/ 2 denotes the permutation of L/ 2 objects, and σ k ( j ) denotes the permutation of all momentum indices of modes in the Fermi sea. Because the norm of each summand is not greater than 1, we see

$$\Sigma _ { 0 } | \langle 0 1 0 1 \dots | \psi \rangle | \leq \frac { 2 ^ { L } } { L ^ { L / 4 } } = e ^ { L ( \ln 2 - \frac { 1 } { 4 } \ln L ) } < e ^ { - \chi _ { 0 } L }, \quad ( \text{B2} )$$

as long as L is large enough there always exist a positive constant χ 0 satisfying the inequality. For the case with only nearest-neighbor hopping, there is only one pair of Fermi points, which does contain the components | 0101 ... ⟩ and | 1010 ... ⟩ . But for the case with next-nearest-neighbor hopping the thing is different by the following claim.

Claim : The overlap of | 1010 ... ⟩ and the two-pair-Fermipoint state | ψ ⟩ as the ground state of the next-nearest-hopping model is zero.

Proof : First we notice the eigenmodes of the two-pairFermi-point state are ˆ c † k = ∑ k e ikx ˆ c † x . The coefficient ⟨ 1010 ... | ψ ⟩ is computed from the Slater determinant of the submatrix containing x = 0 2 4 , , , ... . Restricting the coefficients of ˆ c † k to these sites, we see they form a vector

$$( 1, e ^ { 2 i k }, e ^ { 4 i k }, \dots ) ^ { T }.$$

This vector is linear dependent to another vector with k ′ = k + π . So if there are at least one pair of such k, k ′ in the Fermi sea, then the determinant will be exactly zero. Now we notice an important fact: the next-nearest hopping model has always a pair of such modes in the two disconnected pieces of Fermi sea, ( k, k ′ ) = (0 , π ) , whatever t 1 /t 2 is chosen as long as there are two pairs of Fermi points. Similarly for the overlap ⟨ 0101 ... | ψ ⟩ .

One the other hand, our Wannierized wave functions might contain some component of such | 1010 ... ⟩ states. These are artificial components introduced in the construction of fGMPS. As discussed in the main text, the projected state for g &lt; 0 01 . may have amplified unphysical components. The plot of the C NN ( q ) is shown in Fig. 7.

Figure 7. The static structure factor for g &lt; 0 01 . as a zoomed-in style.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-d4d6d4a0-e5ca-46d6-bf71-f3f8cbeaf4c8.jpg)

## Appendix C: The next-nearest-neighbor hopping model

In momentum space, the model in Eq. (22) reads

$$\hat { c } _ { k } ^ { \dagger } & = \frac { 1 } { \sqrt { L } } \sum _ { i } e ^ { i k x _ { i } } \hat { c } _ { i } ^ { \dagger }, & \quad \begin{array} { c c c c } \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i \\ \hat { H } & = - \sum _ { k } ( t _ { 1 } \cos k + t _ { 2 } \cos 2 k ) \hat { c } _ { k } ^ { \dagger } \hat { c } _ { k }. & \end{array} \begin{array} { c c c c } \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i } \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h h i & \ t h i & \ t h h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i &\ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i &\ t h i & \ t h i & \t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h j & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h h } \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h j } \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h k i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ th i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t g i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t h i & \ t )$$

Our motivation for introducing the next-nearest-neighbor hopping is to obtain a filled Fermi sea state with multiple pairs of Fermi points. To this end, we can first solve for the critical point for which the local band minimum at k = π becomes degenerate in energy with the original Fermi level at k = π/ 2 at half-filling

$$t _ { 1 } \cos \frac { \pi } { 2 } + t _ { 2 } \cos \pi = t _ { 1 } \cos \pi + t _ { 2 } \cos 2 \pi \Rightarrow t _ { 1 } = 2 t _ { 2 }. \ ( C 2 ) \quad \text{the $BR$} \quad..$$

This means we need t 1 &lt; 2 t 2 to obtain a new pair of Fermi points. Denote t 1 /t 2 = t &lt; 2 , and assume the chemical potential is µ , we can solve the following equation to determine the Fermi points:

$$t _ { 1 } \cos k + t _ { 2 } ( 2 \cos ^ { 2 } k - 1 ) + \mu = 0 \quad \ \ ( C 3 )$$

which leads to

$$\cos k _ { 1, 2 } = \frac { - t \pm \sqrt { t ^ { 2 } + 8 ( 1 - \mu _ { t } ) } } { 4 }, \ \mu _ { t } = \mu / t _ { 2 }. \ \ ( C 4 ) \quad \text{Fermi}$$

Suppose the roots of this equation are 0 &lt; k F;1 &lt; k F;2 , the condition for having two pairs of Fermi points is

$$k _ { F ; 2 } - k _ { F ; 1 } = \frac { \pi } { 2 } \Longleftrightarrow \cos k _ { F ; 1 } \cos k _ { F ; 2 } = - \sin k _ { F ; 1 } \sin k _ { F ; 2 }, \quad \text{where} \\ \dots \colon \dots \dots \colon \dots \colon \dots \end{cases}$$

which can be simplified into

$$\mu _ { t } = \frac { t ^ { 2 } } { 4 }. \quad \quad \quad \text{(C6)}$$

Figure 8. (a) The band structure Eq. (C1) in folded Brillouin zone for Wannierization in tree-stack construction of fGMPS; (b) The densitydensity correlation C NN ( q ) for different t 2 with fixed t 1 = 1 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-62a68122-aacd-4f9b-abac-63ac315a8db9.jpg)

Then we know the solution should be

$$\text{now one solu on solu } \mathfrak { c } \\ \cos k _ { 1, 2 } = \frac { - t \pm \sqrt { 8 - t ^ { 2 } } } { 4 }. \quad \quad ( C 7 )$$

The condition of having two different roots is

$$\cos k _ { \mathrm F ; 2 } \leq 0 \leq \cos k _ { \mathrm F ; 1 }, \mu _ { t } \geq t - 1, \quad ( \mathrm C 8 )$$

the latter is coming from the -1 ≤ cos k ≤ 1 , so t ≤ 2 is obtained, consistent with the computation of critical point. The Fermi seas are consisting of modes of | k | &lt; k F;1 and k F;2 &lt; k | | &lt; π , so the ground state is | ψ ⟩ = ∏ k ∈ FS ˆ c † k | 0 ⟩ . The correlation matrix in real space can be computed through the Wick's theorem. The result is

̸

$$\langle \hat { c } _ { i } \hat { c } _ { i + z } ^ { \dagger } \rangle = \frac { 1 } { 2 } \delta _ { z = 0 } - \delta _ { z \neq 0 } \frac { 1 } { L } \sum _ { k \in D F S } e ^ { i k z }. \quad \text{(C9)}$$

By Wick's theorem, the density-density correlation function can be computed.

To study this in tree-stack method[44], we need to first fold the Brillouin zone twice so that all filled modes consist of two fully filled bands (Fig. 8(a)). Then do Wannierization and apply tree-stack construction, we can obtain a local tensor with four spinless fermions.

## Appendix D: Nearest-neighbor Gutzwiller Projector

Besides the Gutzwiller projected ground state of a spinful nearest-neighbor interaction chain, we also considered a free Fermi state projected by a nearest-neighbor Gutzwiller-like projector:

$$\hat { G } _ { \text{nbr} } = \prod _ { x } e ^ { - \alpha \hat { n } _ { x } \hat { n } _ { x + 1 } }, \text{ \quad \ \ } ( \text{D} 1 )$$

where we could impose periodic boundary conditions. To unpack, we may focus on a single term in the product, and notice

↦

$$e ^ { - \alpha \hat { n } _ { x } \hat { n } _ { x + 1 } } \mapsto \begin{cases} 1 & \text{for} ( n _ { x }, n _ { x + 1 } ) = ( 0, 0 ), ( 0, 1 ), ( 1, 0 ) ; \\ e ^ { - \alpha } & \text{for} ( n _ { x }, n _ { x + 1 } ) = ( 1, 1 ) \\ \end{cases}$$

where we denote g := e -α for convenience, and we may write

$$e ^ { - \alpha \hat { n } _ { x } \hat { n } _ { x + 1 } } = 1 + ( g - 1 ) \hat { n } _ { x } \hat { n } _ { x + 1 }. \quad \text{(D3)} \quad T ^ { | }$$

This Gutzwiller-like projector can be understood as follows: consider an 'island' of charge in the state |· · · 011 · · · 10 · · ·⟩ , where we have a string of 1 's with length ℓ . This 'island' contains ℓ -1 nearest-neighbor occupied pairs, and so the projector gives a factor of g ℓ -1 . The overall factor one obtains on any state (in the particle number basis) is the product of the factor from all of such islands of charges. Because our MPS language is based on the fermionic system, we need to turn the projector into a fermionic matrix product form. Its matrix product form is shown in Fig. 9. The matrix can be written in

0

1

1

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-04ecec38-d4ec-4cba-b12d-09ad748d9147.jpg)

1

1

Figure 9. Matrix product form for a Gutzwiller-like projector corresponding to nearest-neighbor density-density interaction. The vertical legs correspond to the physical fermionic Fock space (up = 'out' and down = 'in'). and the horizontal legs are virtual bosonic Hilbert space (of a qubit). All the configurations not shown have 0 weight.

a compact form:

$$\hat { T } ^ { [ x ] } = \sum _ { i _ { x }, j _ { x } } \, \frac { 1 } { \longleftarrow \left ( \sum _ { j _ { x } } ^ { 1 } \right ) _ { j _ { x } } } = \sum _ { i _ { x }, j _ { x } } | i _ { x } \rangle \, \hat { O } _ { i _ { x } j _ { x } } ^ { [ x ] } \, \langle j _ { x } | \,, \quad \text{(D4)}$$

where

$$\hat { O } _ { i _ { x } j _ { x } } ^ { [ x ] } = ( a _ { L } \hat { n } _ { x } ) ^ { i _ { x } } ( a _ { R } \hat { n } _ { x } ) ^ { j _ { x } } ; \ \ a _ { L } a _ { R } = g - 1. \quad ( \text{D5} )$$

where i x j x ∈ { 0 1 , } . As a convenient choice, set a L = ( g -1) , a R = 1 :

$$\hat { O } _ { 0 0 } ^ { [ x ] } = \hat { 1 } ; \ \hat { O } _ { 0 1 } ^ { [ x ] } = \hat { n } _ { x } ; \ \ & \cdot \ \hat { O } _ { 1 0 } ^ { [ x ] } = \hat { O } _ { 1 1 } ^ { [ x ] } = ( g - 1 ) \hat { n } _ { x }$$

Then a continuous chain of nearest-neighbor Gutzwiller projectors is recovered by the product of ˆ T [ x ] :

$$& \text{y write} \quad \text{Then a continuous chain of nearest-neighbor Gutzwiller pro-} \\ & \text{jectors is recovered by the product of $\tilde{T}^{[x]}$} \colon \\ & \quad \hat { T } ^ { [ 0 ] } \hat { T } ^ { [ 1 ] } \dots \hat { T } ^ { [ \ell ] } \\ & \quad = \quad \sum _ { i _ { 0 }, i _ { 0 }, \cdots i _ { \ell }, j _ { \ell } } | i _ { 0 } \rangle \, \hat { O } ^ { [ 0 ] } _ { i _ { 0 } j _ { 0 } } \, \langle j _ { 0 } | \, i _ { 1 } \rangle \hat { O } ^ { [ 1 ] } _ { i _ { 1 } j _ { 1 } } \, \langle j _ { 1 } | \, i _ { 2 } \rangle \cdots \hat { O } ^ { [ \ell ] } _ { i _ { \ell } j _ { \ell } } \, \langle j _ { \ell } | \\ & \text{allows} \colon \\ & \text{to} \quad \sum _ { i _ { 0 }, j _ { \ell } } | i _ { 0 } \rangle \, ( ( g - 1 ) \hat { n } _ { 0 } ) ^ { i _ { 0 } } \\ & \text{to} \quad \begin{array} { c } \sum _ { i _ { 0 }, j _ { \ell } } | i _ { 0 } \rangle \, \hat { O } ^ { [ 0 ] } _ { i _ { 0 } j _ { 1 } } \, \hat { O } ^ { [ 1 ] } _ { i _ { 1 } j _ { 2 } } \cdots \hat { O } _ { i _ { \ell } j _ { \ell } } \, \langle j _ { \ell } | \\ \end{array} \\ & \text{to} \quad \begin{array} { c } \sum _ { i _ { 0 }, i _ { 0 }, i _ { 2 } \cdots i _ { \ell }, j _ { \ell } } | i _ { 0 } \rangle \, \hat { O } ^ { [ 0 ] } _ { i _ { 0 } j _ { 1 } } \, \hat { O } ^ { [ 1 ] } _ { i _ { 1 } j _ { 2 } } \cdots \hat { O } _ { i _ { \ell } j _ { \ell } } \, \langle j _ { \ell } | \\ \end{array} \\ & \text{there the contraction of bond states satisfies} \\ & \quad \langle j _ { x } | \, i _ { x + 1 } \rangle = \delta _ { j _ { x }, i _ { x + 1 } } \quad \begin{array} { c } ( ( g - 1 ) \hat { n } _ { n } ) ^ { i _ { n } } \end{array} \right ] ( \hat { n } _ { \ell } ) ^ { j _ { \ell } } \, \langle j _ { \ell } | \\ \end{array}$$

where the contraction of bond states satisfies

$$\langle j _ { x } | \, i _ { x + 1 } \rangle = \delta _ { j _ { x }, i _ { x + 1 } } \text{ \quad \ \ } \text{(D8)}$$

In fermionic form the purification of the operator ˆ T [ x ] is a composition of the purified operator and bond states. The purified operator is defined as

$$| \hat { O } _ { i _ { x } j _ { x } } ^ { [$$

where the reference state | D [ x ] ⟩ is defined on the Fock space of the physical fermion at x together with its auxiliary. Note that D [ x ] has charge 1 and is odd-parity. The same is true for ˆ O [ x ] i x x j , as all the ˆ O [ x ] i x x j are number conserving. The fermionic version of the bond states should satisfy

$$\varkappa ^ { \dots } \cdots \\ \langle K _ { x,$$

To ensure the bond states follow the fermionic statistics, we use two fermions to represent the bond:

$$\text{$\text{rmmos} \, \text{or} \, \text$$

which matches our convention of contraction kernel

$$\text{$num$ means you common on conquenent} \\ | K \rangle =$$

Combining the discussions above, we finally obtain the fermionic state corresponding to the local tensor

$$\begin{array} {$$

## Appendix E: Contraction and Transfer matrix

In this section we first derive the formula Eq. (13), and then we show the derivation of Eq. (14). For simplicity, we use the the vector symbol a † , p † to denote the product of a sequence of creation operators, for instance, a † 0 = (0 1 1) , , † ≡ ˆ a † 0;2 ˆ a † 0;3 representing the creation of mode 2 and 3 in the bond space labelled by a 0 . The conjugate of these symbols a 0 represent the conjugate of the creation operators ˆ a 0;3 ˆ a 0;2 .

## 1. Contraction via matrix product

To derive the Eq. (13), we consider a kernel contraction of local tensors,

$$\dots \dots \dots, & \quad \langle K _ { b, b ^ { \prime } } | \left ( | \Psi _ { a, b } \rangle \otimes | \Phi _ { b ^ { \prime }, c } \rangle \right ) & \quad \text{$r$-\quad} \\ = & K _ { b, b ^ { \prime } } ^ { * } \Psi _ { a b } \Phi _ { b ^ { \prime } c } \left \langle 0 | b ^ { \prime } b a ^ { \dagger } b ^ { \dagger } b ^ { \prime \dagger } c ^ { \dagger } \left | 0 \right \rangle \\ = & K _ { b, b ^ { \prime } } ^ { * } \Psi _ { a b } \Phi _ { b ^ { \prime } c } ( - 1 ) ^ { n _ { a } N _ { b } } a ^ { \dagger } c ^ { \dagger } \left | 0 \right \rangle & \quad \text{$thei$-\quad} \\ = & ( P _ { a } ) ^ { N _ { b } } \Psi _ { a b } K _ { b, b ^ { \prime } } \Phi _ { b ^ { \prime } c } a ^ { \dagger } c ^ { \dagger } \left | 0 \right \rangle. & \quad \text{$state$-\quad} \\ \tau \cdots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

To show the kernel purification is compatible with kernel contraction, we need some properties of the contraction kernel | K ⟩ and purification kernel | D ⟩ . The matrix product satisfies the following identity

$$\sum _ { m ^ { \prime } } D _ { m, m ^ { \prime } } K _ { m ^ { \prime }, \mathfrak { a } } = ( - 1 ) ^ { n _ { \mathfrak { a } } ( N _ { m } ) } = [ \hat { P } ] _ { \mathfrak { m a } } ^ { N _ { \mathfrak { m } } }. \quad ( E 2 ) \quad \ e q. ( \mathfrak { l } _ { \mathfrak { m a r i } }$$

Now consider an operator ˆ O m = O γ γ ˆ , ˆ γ being a product of some Majorana operators, acting on the modes m . Then

$$\text{ some way to run} \, \text{s, using on one moves} \, m \, \text{m} \\ \sum _ { m ^ { \prime } a } \langle K _ { m ^ { \prime } a } | \hat { O } _ { m } | D _ { m m ^ { \prime } } \rangle \otimes _ { f } | \psi _ { a b } \rangle \\ = & K _ { m ^ { \prime } a } O _ { \gamma } D _ { m m ^ { \prime } } \psi _ { a b } \left \langle 0 | \text{am} ^ { \prime } \gamma _ { m } ^ { \dagger } m ^ { \prime } \text{a} ^ { \dagger } b ^ { \dagger } | 0 \right \rangle \\ = & O _ { \gamma } D _ { m m ^ { \prime } } K _ { m ^ { \prime } a } \psi _ { a b } ( - 1 ) ^ { N _ { \gamma } ( n _ { a } + n _ { m ^ { \prime } } ) } \\ & \times ( - 1 ) ^ { n _ { m } ( n _ { m ^ { \prime } } + n _ { a } ) } \hat { \gamma } _ { m ^ { \dagger } b ^ { \dagger } | 0 \right \rangle, \quad \text{with} \, \cdot \\ \text{in} \, \infty \, \infty \, + \infty \, - \, N \text{ is fixed and the only non-synchasing} \quad \text{ossil}$$

since n m + n m ′ = N m is fixed, and the only non-vaninshing elements in DK are diagonal elements, thus n a = n m , the expression can be simplified to

$$& ( - 1 ) ^ { N _ { \gamma } N _ { m } } O _ { \gamma } \psi _ { \text{ab} } ( - 1 ) ^ { n _ { \text{a} } ( N _ { m } ) } ( - 1 ) ^ { n _ { \text{m} } N _ { m } } \hat { \gamma } \text{m} ^ { \dagger } b ^ { \dagger } | 0 \rangle \\ = & ( - 1 ) ^ { N _ { \gamma } N _ { m } } O _ { \gamma } \psi _ { \text{mb} } \hat { \gamma } \text{m} ^ { \dagger } b ^ { \dagger } | 0 \rangle \\ & & \text{where}$$

So given the operator has an even parity, or the number of modes N m it acts on is even, the result above is equal to the ˆ O m | ψ mb ⟩ . If ˆ O has an odd parity and N m is odd, then it is up to a negative sign. For a generic ˆ O with terms of mixed parities, we can distinguish the parity sectors and compensate the signs accordingly.

## 2. Transfer matrix

The transfer matrix in this paper is defined as the regrouped matrix of a reduced density matrix tracing out the physical modes. Let | ϕ L,P,R ⟩ be a local fermion state with N L + N P + N R modes, the reduced density matrix is

$$\hat { \rho } _ { L R } = \text{Tr} _ { P } \left [ \left | \phi _ { L, P, R } \right \rangle \left \langle \phi _ { L, P, R } \right | \right ] \text{ \quad \ \ } \text{(E5)}$$

In the computational basis, it is a 2 N N L R × 2 N N L R matrix. Then regrouping the elements leads to a transfer matrix M ll ′ ,rr ′ = (ˆ ρ LR lr,l r ) ′ ′ . This definition is different from the usual transfer matrix defined in the literature, as the partial trace of a fermionic density matrix involves some sign differences. We now show the computation of observables can be represented as computing the trace of a product of transfer matrices, without the need of inserting the parity operators ˆ P N B in Eq. (13). Suppose there are L local tensors | Φ a p a 0 1 1 ⟩ , | Φ ˜ a p a 1 2 2 ⟩ , ... | Φ ˜ a L -1 p L a L ⟩ and denote their tensor product as | Φ = Φ ⟩ | a p a 0 1 1 ⟩ ⊗ | Φ ˜ a p a 1 2 2 ⟩ ⊗ ... ⊗| Φ ˜ a L -1 p L a L ⟩ , we first consider the computation of normalization

$$\mathcal { N } _ { \Phi } = \langle \Phi | K \rangle \langle K | \Phi \rangle = \text{Tr} \left [ \langle K | \Phi \rangle \langle \Phi | K \rangle \right ],$$

where | K ⟩ = | K a a 1 ˜ 1 ⟩⊗| K a a 2 ˜ 2 ⟩⊗ ... ⊗| K a L a 0 ⟩ is the product state of all contraction kernels. The states | K ⟩ and | Φ ⟩ are denoted as

$$| K \rangle & = ( K _ { \mathbf a _ { 0 } \dots } ) \mathbf a _ { 1 } ^ { \dagger } \tilde { \mathbf a } _ { 1 } ^ { \dagger } \mathbf a _ { 2 } ^ { \dagger } \tilde { \mathbf a } _ { 2 } ^ { \dagger } \cdots \mathbf a _ { L } ^ { \dagger } \mathbf a _ { 0 } ^ { \dagger } \left | 0 \right \rangle ; \\ | \Phi \rangle & = ( \Phi _ { \mathbf a _ { 0 } \dots } ) \mathbf a _ { 0 } ^ { \dagger } \mathbf p _ { 1 } ^ { \dagger } \mathbf a _ { 1 } ^ { \dagger } \cdots \tilde { \mathbf a } _ { L - 1 } ^ { \dagger } \mathbf p _ { L } ^ { \dagger } \mathbf a _ { L } ^ { \dagger } \left | 0 \right \rangle,$$

where the coefficients are also labelled by these vectors. The Eq. (E6) is done by first obtaining the local reduced density matrices and then trace with the kernel state,

$$\text{c} \text{ of } \quad \mathcal { N } _ { \phi } & = \langle K | \text{ Tr} _ { P } \left [ \langle \Phi | \Phi \rangle \right ] | K \rangle \\ & = K _ { \mathfrak { a } _ { 0 } \cdots } K _ { \mathfrak { a } _ { 0 } \cdots } ^ { * } M _ { \mathfrak { a } _ { 0 } \cdots } ^ { \mathfrak { a } _ { 0 } \cdots } \langle 0 | \underbrace { \mathfrak { a } _ { 0 } \mathfrak { a } _ { L } \cdots \tilde { \mathfrak { a } } _ { 1 } \mathfrak { a } _ { 1 } } _ { K } \underbrace { \mathfrak { a } _ { 0 } ^ { \dagger } \tilde { \mathfrak { a } } _ { 1 } ^ { \dagger } \tilde { \mathfrak { a } } _ { 1 } ^ { \dagger } \tilde { \mathfrak { a } } _ { 1 } ^ { \dagger } \tilde { \mathfrak { a } } _ { 1 } ^ { \dagger } \cdots \mathfrak { a } _ { L } ^ { \dagger } } _ { \Phi \, \text{reduced } } | 0 \rangle \\ & \quad \times \langle 0 | \underbrace { \mathfrak { a } _ { L } \cdots \tilde { \mathfrak { a } } _ { 1 } \mathfrak { a } _ { 1 } ^ { \prime } } _ { \Phi \, \text{reduced } } _ { K } ^ { \tilde { \mathfrak { a } } ^ { \prime \dagger } \mathfrak { a } ^ { \prime \dagger } _ { 1 } \cdots \mathfrak { a } ^ { \prime \dagger } _ { L } \mathfrak { a } ^ { \prime \dagger } _ { 0 } } _ { K } | 0 \rangle \,, \\ & \quad \text{with all ransact} \text{ in} \text{ in} \text{ in} \text{ with } \text{ Now } \text{the only}$$

with all repeated indices summed up implicitly. Now the only possibly non-trivial signs coming from these creation and annihilation operators are those due to the a 0 / a † 0 and a L / a † L , which is

$$& ( - 1 ) ^ { n _ { \mathfrak { a } _ { 0 } } ( N _ { \mathfrak { a } _ { 1 } } + \dots + N _ { \mathfrak { a } _ { L - 1 } } ) + n _ { \mathfrak { a } _ { 0 } } n _ { \mathfrak { a } _ { L } } } \\ & \quad \times \left ( - 1 \right ) ^ { n _ { \mathfrak { a } _ { 0 } ^ { \prime } } ( N _ { \mathfrak { a } _ { 1 } ^ { \prime } } + \dots + N _ { \mathfrak { a } _ { L - 1 } ^ { \prime } } ) + n _ { \mathfrak { a } _ { 0 } ^ { \prime } } n _ { \mathfrak { a } _ { L } ^ { \prime } } }$$

where we used the condition that each kernel state is particlenumber-conserving. A further simplification is n a 0 ≡ n a L + N a 0 mod 2, so the above sign is

$$( - 1 ) ^ { ( n _ { \mathfrak { a } _ { 0 } } + n _ { \mathfrak { a } _ { 0 } ^ { \prime } } ) ( N _ { \mathfrak { a } } ^ { t } + 1 ) }, \quad \quad \quad ( \mathbb { E } 1 0 )$$

where N t a is one half of the total number of bond modes of | Φ ⟩ . So the computation of the normalization is reduced to

$$\mathcal { N } _ { \Phi } & = M _ { \mathfrak { a } _ { 0 } ^ { \prime }, \mathfrak { a } _ { 1 } ^ { \prime } } ^ { \mathfrak { a } _ { 1 }, \mathfrak { a } _ { 1 } } \mathcal { K } _ { \mathfrak { a } _ { 1 } ^ { \prime }, \mathfrak { a } _ { 1 } ^ { \prime } } ^ { \mathfrak { a } _ { L - 1 }, \mathfrak { a } _ { L } } \cdots M _ { \mathfrak { a } _ { L - 1 } ^ { \prime }, \mathfrak { a } _ { L } ^ { \prime } } ^ { \mathfrak { a } _ { L - 1 }, \mathfrak { a } _ { L } } \mathcal { K } _ { \mathfrak { a } _ { L } ^ { \prime }, \mathfrak { a } _ { 0 } } ^ { \mathfrak { a } _ { L }, \mathfrak { a } _ { 0 } } \left [ \Xi ^ { N _ { \mathfrak { a } } ^ { t } + 1 } _ { \mathfrak { a } _ { 0 } ^ { \prime } } \right ] _ { \mathfrak { a } _ { 0 } ^ { \prime } } ^ { \mathfrak { a } _ { L } } \\ \text{ed} \colon & = \text{tr} \left [ M _ { 0 } \mathcal { K } _ { 0, 1 } M _ { 1 } \mathcal { K } _ { 1, 2 } \cdots M _ { L } \mathcal { K } _ { L, 0 } \Xi ^ { N _ { \mathfrak { a } } ^ { t } + 1 } \right ],$$

where we use K = K ⊗ K ∗ to denote the Kronecker product of the matrices K a i , ˜ a i and K ∗ a ′ i , ˜ a ′ i , and Ξ to denote the Kronecker product of the parity matrices ˆ P a 0 , a 0 and ˆ P a ′ 0 , a ′ 0 . The same derivation can be done for the computation of observables. With the insertion of some operators, the coefficients of Φ will be modified, but as long as the transfer matrix is defined as the partial trace of physical modes in local tensors, the sign issue is always the same as derived above.

## Appendix F: Finite-bond scaling

In the main text we discussed the tree-stack construction of translation-invariant local tensors. This process involves finding a fGMPS representation, which includes a radius r truncation around the Wannier center. The truncation controls the number of bonds of local tensor, affecting the momentum distribution n k near the Fermi surface. We define ∆ k to describe the smearing:

$$\Delta k ( \epsilon ) \coloneqq \frac { N _ { \epsilon } ( n _ { k } ) } { L }$$

where N n ϵ ( k ) is the range of momentum k satisfying ϵ &lt; n k &lt; 1 -ϵ and L is the system size. A smaller ∆ ( ) k ϵ

Figure 10. The computation here is based on a half-filled Fermi sea of spin1 2 / electrons in one dimension system. (a) The momentum density n k distribution performance near π/ 2 for different stack tree radius. (b) The ∆ ( k ϵ = 0 01) . for truncation radius r = 8 16 32 64 , , , for a system with system size L = 512 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00858v2-figure-0c48b19d-0739-4d96-a067-39bf1dfe4a7d.jpg)

- [1] M. C. Gutzwiller, Effect of correlation on the ferromagnetism of transition metals, Phys. Rev. Lett. 10 , 159 (1963).
- [2] M. C. Gutzwiller, Correlation of electrons in a narrow s band, Phys. Rev. 137 , A1726 (1965).
- [3] W. F. Brinkman and T. M. Rice, Application of gutzwiller's variational method to the metal-insulator transition, Phys. Rev. B 2 , 4302 (1970).

represents a sharper Fermi surface.

In Fig. 10, we show the momentum density distribution and ∆ k for different truncation length in a system with size L = 512 model. As the truncation radius r increase, the momentum density n k near k = π/ 2 moves closer to n k = 0 1 , , resulting in a decrease in ∆ ( k ϵ = 0 01) . . Ultimately, in the thermodynamic limit ∆ k is expected to vanish.

- [4] D. Vollhardt, Normal 3 He : an almost localized fermi liquid, Rev. Mod. Phys. 56 , 99 (1984).
- [5] F.-C. Zhang, C. Gros, T. M. Rice, and H. Shiba, A renormalised hamiltonian approach to a resonant valence bond wavefunction, Superconductor Science and Technology 1 , 36 (1988).
- [6] D. Vollhardt, P. G. van Dongen, F. Gebhard, and W. Metzner, Gutzwiller-type wave functions for correlated fermions, Mod-

- ern Physics Letters B 4 , 499 (1990).
- [7] W. Metzner and D. Vollhardt, Ground-state properties of correlated fermions: Exact analytic results for the gutzwiller wave function, Phys. Rev. Lett. 59 , 121 (1987).
- [8] B. Edegger, V. N. Muthukumar, and C. Gros, Gutzwiller-rvb theory of high-temperature superconductivity: Results from renormalized mean-field theory and variational monte carlo calculations, Advances in Physics 56 , 927 (2007).
- [9] X. Deng, X. Dai, and Z. Fang, Lda+ gutzwiller method for correlated electron systems, Europhysics Letters 83 , 37008 (2008).
- [10] J. N. Zhuang, L. Wang, Z. Fang, and X. Dai, Fast impurity solver based on gutzwiller variational approach, Phys. Rev. B 79 , 165114 (2009).
- [11] X. Deng, L. Wang, X. Dai, and Z. Fang, Local density approximation combined with gutzwiller method for correlated electron systems: Formalism and applications, Phys. Rev. B 79 , 075114 (2009).
- [12] J. B¨ unemann, The Gutzwiller Variational Theory and Related Methods for Correlated Electron Systems , Ph.D. thesis, Habilitationsschrift, Marburg (2009).
- [13] Z. Cheng and C. A. Marianetti, Foundations of variational discrete action theory, Physical Review B 103 , 195138 (2021).
- [14] Z. Cheng and C. A. Marianetti, Variational discrete action theory, Physical Review Letters 126 , 206402 (2021).
- [15] Z. Cheng and C. A. Marianetti, Precise ground state of multiorbital mott systems via the variational discrete action theory, Physical Review B 106 , 205129 (2022).
- [16] Z. Cheng and C. A. Marianetti, Gauge constrained algorithm of variational discrete action theory at n= 3 for the multiorbital hubbard model, Physical Review B 108 , 035127 (2023).
- [17] R. Or´ us, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of physics 349 , 117 (2014).
- [18] R. Or´ us, Tensor networks for complex quantum systems, Nature Reviews Physics 1 , 538 (2019).
- [19] Z.-X. Liu, Y. Zhou, H.-H. Tu, X.-G. Wen, and T.-K. Ng, Gutzwiller projected wave functions in the fermionic theory of s= 1 spin chains, Physical Review B 85 , 195144 (2012).
- [20] Y.-H. Wu, L. Wang, and H.-H. Tu, Tensor network representations of parton wave functions, Phys. Rev. Lett. 124 , 246401 (2020).
- [21] G. Petrica, B.-X. Zheng, G. K.-L. Chan, and B. K. Clark, Finite and infinite matrix product states for gutzwiller projected meanfield wave functions, Phys. Rev. B 103 , 125161 (2021).
- [22] S. Budaraju, D. Poilblanc, and S. Niu, Simulating chiral spin liquids with fermionic projected entangled paired states, arXiv preprint arXiv:2403.01243 (2024).
- [23] F. Gebhard and D. Vollhardt, Correlation functions for interacting fermions in the gutzwiller ansatz, Physical Review B 38 , 6911 (1988).
- [24] W. Metzner and D. Vollhardt, Analytic calculation of groundstate properties of correlated fermions with the gutzwiller wave function, Physical Review B 37 , 7382 (1988).
- [25] M. Kollar and D. Vollhardt, Exact analytic results for the gutzwiller wave function with finite magnetization, Physical Review B 65 , 155121 (2002).
- [26] U. Schollw¨ ock, The density-matrix renormalization group, Rev. Mod. Phys. 77 , 259 (2005).
- [27] U. Schollw¨ ock, The density-matrix renormalization group in the age of matrix product states, Annals of physics 326 , 96 (2011).
- [28] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions, arXiv preprint cond-mat/0407066 (2004).
- [29] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, arXiv preprint quantph/0608197 (2006).
- [30] F. Verstraete, V. Murg, and J. I. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Advances in physics 57 , 143 (2008).
- [31] Z.-C. Gu, M. Levin, and X.-G. Wen, Tensor-entanglement renormalization group approach as a unified method for symmetry breaking and topological phase transitions, Phys. Rev. B 78 , 205116 (2008).
- [32] R. Or´ us, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349 , 117 (2014).
- [33] J. I. Cirac and F. Verstraete, Renormalization and tensor product states in spin chains and lattices, Journal of physics a: mathematical and theoretical 42 , 504004 (2009).
- [34] M. T. Fishman and S. R. White, Compression of correlation matrices and an efficient method for forming matrix product states of fermionic gaussian states, Physical Review B 92 , 075132 (2015).
- [35] A. N¨ ußeler, I. Dhand, S. F. Huelga, and M. B. Plenio, Efficient construction of matrix-product representations of many-body gaussian states, Phys. Rev. A 104 , 012415 (2021).
- [36] H.-K. Jin, H.-H. Tu, and Y. Zhou, Density matrix renormalization group boosted by gutzwiller projected wave functions, Phys. Rev. B 104 , L020409 (2021).
- [37] H.-K. Jin, H.-H. Tu, and Y. Zhou, Efficient tensor network representation for gutzwiller projected states of paired fermions, Phys. Rev. B 101 , 165135 (2020).
- [38] Q. Yang, X.-Y. Zhang, H.-J. Liao, H.-H. Tu, and L. Wang, Projected d-wave superconducting state: A fermionic projected entangled pair state study, Physical Review B 107 , 125128 (2023).
- [39] J.-W. Li, J. von Delft, and H.-H. Tu, U (1)-symmetric gaussian fermionic projected entangled paired states and their gutzwiller projection, Physical Review B 107 , 085148 (2023).
- [40] C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Fermionic projected entangled pair states, Physical Review A 81 , 052338 (2010).
- [41] J. Dubail and N. Read, Tensor network trial states for chiral topological phases in two dimensions and a no-go theorem in any dimension, Phys. Rev. B 92 , 205307 (2015).
- [42] N. Schuch and B. Bauer, Matrix product state algorithms for gaussian fermionic states, Phys. Rev. B 100 , 245121 (2019).
- [43] H.-K. Jin, R.-Y. Sun, Y. Zhou, and H.-H. Tu, Matrix product states for hartree-fock-bogoliubov wave functions, Phys. Rev. B 105 , L081101 (2022).
- [44] Y. He, K. Li, Y. Zhang, and H. C. Po, Stacked tree construction for free-fermion projected entangled pair states, Physical Review Research 6 , L022016 (2024).
- [45] S. Bravyi, Lagrangian representation for fermionic linear optics, arXiv preprint quant-ph/0404180 (2004).
- [46] I. Peschel, Entanglement in solvable many-particle models, Brazilian Journal of Physics 42 , 267 (2012).
- [47] J. Surace and L. Tagliacozzo, Fermionic gaussian states: an introduction to numerical approaches, SciPost Physics Lecture Notes , 054 (2022).
- [48] A. Jamiołkowski, Linear transformations which preserve trace and positive semidefiniteness of operators, Reports on mathematical physics 3 , 275 (1972).
- [49] M.-D. Choi, Completely positive linear maps on complex matrices, Linear algebra and its applications 10 , 285 (1975).
- [50] We use a JW convention such that, taking four fermion modes as an example, the many-body states are ordered according to

| 0000 ⟩ , | 0001 ⟩ , | 0010 ⟩ , ... .

- [51] M. M. Wolf, Violation of the entropic area law for fermions, Phys. Rev. Lett. 96 , 010404 (2006).
- [52] Q. Mortier, N. Schuch, F. Verstraete, and J. Haegeman, Tensor networks can resolve fermi surfaces, Physical Review Letters 129 , 206401 (2022).
- [53] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in quantum critical phenomena, Phys. Rev. Lett. 90 , 227902 (2003).
- [54] D. Gioev and I. Klich, Entanglement entropy of fermions in any dimension and the widom conjecture, Phys. Rev. Lett. 96 , 100503 (2006).
- [55] S. D. Geraedts, M. P. Zaletel, R. S. Mong, M. A. Metlitski, A. Vishwanath, and O. I. Motrunich, The half-filled landau level: The case for dirac composite fermions, Science 352 , 197 (2016).
- [56] P. Fazekas, Wigner crystallization as a possibility for the metalnonmetal transition in magnetite, Solid State Communications 10 , 175 (1972).
- [57] P. Fazekas, Variational studies of strongly correlated electron systems, Physica Scripta 29 , 125 (1989).
- [58] G. Seibold and E. Sigmund, Gutzwiller and slave-boson methods for intersite coulomb interactions, Zeitschrift f¨r Physik B u Condensed Matter 101 , 405 (1997).
- [59] M. Sherafati, S. Satpathy, and D. Pettey, Gutzwiller variational method for intersite coulomb interactions: The spinless fermion model in one dimension, Phys. Rev. B 88 , 035114 (2013).
- [60] G. D. Mahan, Many-particle physics (Springer Science &amp; Business Media, 2013).
- [61] G. Vidal, Entanglement renormalization, Phys. Rev. Lett. 99 , 220405 (2007).
- [62] X.-L. Qi, Exact holographic mapping and emergent space-time geometry, arXiv preprint arXiv:1309.6282 (2013).
- [63] S. L. Wong, K. C. Pang, and H. C. Po, Zipper entanglement renormalization for free fermions, arXiv preprint arXiv:2206.11761 (2022).
- [64] T. Liu, Y.-H. Wu, H.-H. Tu, and T. Xiang, Efficient conversion from fermionic gaussian states to matrix product states, arXiv preprint arXiv:2408.01155 (2024).
Kangle Li Yan-Bai Zhang Hoi Chun Po
Constructive Fermionic Matrix Product States for Projected Fermi Sea
Projected wave functions offer a means for incorporating local correlation
effects in gapless electronic phases of matter like metals. Although such wave
functions can be readily specified formally, it is challenging to compute their
associated physical observables. Tensor network approaches offer a modern
numerical method for this task. In this work, we develop and demonstrate a
constructive tensor-network approach for obtaining physical quantities, like
fermion two-point functions and density-density correlation functions, for
one-dimensional projected Fermi sea states. We benchmark our method against
exact analytical results for spin-1/2 electrons subjected to the Gutzwiller
projection, and then present results on spinless fermions with nearest-neighbor
repulsion. For a state with two pairs of Fermi points, we reveal a
correlation-tuning of the characteristic wave vector of the charge density
modulation in the system.