## Efficient conversion from fermionic Gaussian states to matrix product states

Tong Liu, 1, 2 Ying-Hai Wu, 3 Hong-Hao Tu, 4, ∗ and Tao Xiang 1, 5, 2, †

1 Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China 3 School of Physics and Wuhan National High Magnetic Field Center, Huazhong University of Science and Technology, Wuhan 430074, China 4 Physics Department and Arnold Sommerfeld Center for Theoretical Physics,

Ludwig-Maximilians-Universit¨t a M¨ unchen, 80333 Munich, Germany

5 Beijing Academy of Quantum Information Sciences, Beijing 100190, China

Fermionic Gaussian states are eigenstates of quadratic Hamiltonians and widely used in quantum many-body problems. We propose a highly efficient algorithm that converts fermionic Gaussian states to matrix product states. It can be formulated for finite-size systems without translation invariance, but becomes particularly appealing when applied to infinite systems with translation invariance. If the ground states of a topologically ordered system on infinite cylinders are expressed as matrix product states, then the fixed points of the transfer matrix can be harnessed to filter out the anyon eigenbasis, also known as minimally entangled states. This allows for efficient computation of universal properties such as entanglement spectrum and modular matrices. The potential of our method is demonstrated by numerical calculations in two chiral spin liquids that have the same topological orders as the bosonic Laughlin and Moore-Read states, respectively. The anyon eigenbasis for the first one has been worked out before and serves as a useful benchmark. The anyon eigenbasis of the second one is, however, not transparent and its successful construction provides a nontrivial corroboration of our method.

## I. INTRODUCTION

At the fundamental level, particles in a quantum manybody system interact with each other. However, the essential physics in many cases are captured by a free particle description, notably in the framework of mean-field theory. A plethora of fascinating states can be generated when free fermions roam in diverse lattice potentials and repel each other solely according to Pauli's exclusion principle. In general, the Hamiltonians of free fermions are quadratic in terms of the creation and annihilation operators so they can be diagonalized easily using singleparticle basis transformations. The solutions can be expressed as fermionic Gaussian states (FGSs) whose correlations obey the Wick's theorem [1-4]. If strong correlations invalidate mean-field treatment, many-body problems become very difficult due to the exponential growth of the Hilbert space.

projected entangled pair states (PEPSs) [6-11]. The two methods have their respective strengths and weaknesses but both can provide successful variational Ans¨tze. Due a to the apparent differences, these approaches have been almost independent from each other when they were developed and applied in concrete problems.

In a variety of strongly correlated phases, mean-field theory can still provide useful insight when it is applied on fictitious partons that emerge from breaking each physical object into multiple pieces. FGSs are still useful when the partons are fermionic and described by quadratic Hamiltonians [5]. It is necessary to impose physical constraint using Gutzwiller projection such that partons are glued together to become physical objects. While FGSs of partons are simple, the final wave functions are generally not easy to analyze. Another fruitful approach for strongly correlated systems is tensor networks that include matrix product states (MPSs) and

∗ h.tu@lmu.de

† txiang@iphy.ac.cn

It is natural to explore if FGSs and tensor networks can be combined to harness their advantages in collaboration. One attempt is to design hybrid ansatz that contains both of them as ingredients [12, 13]. Another route is converting (physically motivated or post-optimized) FGSs into tensor networks that may serve as initial states for numerical simulations to facilitate better convergence. Along this line of thought, several methods have been developed to convert FGSs into MPSs [14-19] and they are indeed useful in density matrix renormalization group (DMRG) calculations [19-24]. In particular, feeding a DMRG program with suitable initial states could simplify the search for the degenerate ground states in a two-dimensional (2D) system with intrinsic topological order [21, 22, 24]. The PEPS representations of FGSs have also been investigated [25-28] and their application in optimization algorithms is anticipated [29-36].

In this work, we develop an algorithm that efficiently converts FGS to MPS by making full usage of the FGS correlation matrix formalism. While this method is capable of dealing with cases without translation symmetry, its potential can be most clearly demonstrated using translationally invariant infinite matrix product states (iMPSs). It can be applied when each unit cell contains multiple sites cells such as 2D systems on infinite cylinders. We have achieved significant reduction of computational cost and constructed reliable iMPSs for very large systems. This method also leads to another im-

2

portant discovery about spontaneous symmetry breaking or topological order. It is well-known that the ground states of such systems are degenerate and their MPSs are 'non-injective' in the sense that local matrices generally have more than one blocks in the canonical form [37]. Nevertheless, injective MPSs can be obtained from noninjective ones by projecting onto each irreducible block in the virtual space. This progress enriches our toolbox for characterizing topological order in Guzwiller pojected FGSs. For a given iMPS representation, ground-state degeneracy manifests itself in the multiplicity of the leading eigenvectors ('fixed points') of the iMPS transfer matrix. If we perform the projection to constructed an injective MPS, its bipartite entanglement entropy would be minimized. For 2D topological order, injective MPSs with this property are called anyon eigenstates (or 'minimally entangled states' [38, 39]), which are particularly useful for extracting topological data [40-42]. As benchmark examples, this procedure is employed to find the anyon eigenbasis in two chiral spin liquid models and characterize their topological orders. The first system has Abelian SU(2) 1 topological order and its anyon eigenbasis has not been constructed before. The second system has nonAbelian SU(2) 2 topological order but its anyon eigenbasis is unknown and clearly illustrates the utility of our approach.

The rest of this paper is organized as follows. In Sec. II, we describe the algorithm in detail. It has three main steps: (i) Get the correlation matrix of the FGS to be converted; (ii) Decompose the correlation matrix of the whole system into local projectors and virtual bonds; (iii) Perform mode decimation and reconfigure local projectors and virtual bonds as local tensors. In Sec. III, the algorithm is applied on two chiral spin liquid models to construct their anyon eigenbases from which entanglement spectra are computed. This paper is summarized with an outlook in Sec. IV.

## II. METHOD

In this section, we present our algorithm in detail. This section is organized as follows. We first introduce FGSs and the correlation matrix formalism in Sec. II A. Next we derive and solve the recursion equation for the correlation matrix of the MPS local tensor in Sec. II B. Then, we discuss how to generate local tensors from the local correlation matrix and perform truncation in Sec. II C. Finally, we propose the mode decimation scheme to accelerate the computation of the MPS local tensor in Sec. II D and discuss how to implement the anyon eigenbasis filtration using iMPS in Sec. II E.

## A. Fermionic Gaussian states and correlation matrix formalism

Since we will extensively use FGSs in this work, we first give a brief review of the formulation and introduce necessary concepts. For notational simplicity, we restrict ourselves in Sec. II to particle-number-conserving FGSs [termed as U(1)-FGSs]. The extension to the more general case with pairing is straightforward but requires some adaptions of the formalism (see, e.g., Ref. [19]).

We start with a fermionic lattice system. The fermionic creation and annihilation operators are written as c † j and c j , respectively. When presenting the algorithm, we treat j as the site index, with j = 1 , . . . , N . For cases where fermions have internal degrees of freedom (e.g., spin or orbital), one may view j as the combination of site and internal state indices, e.g., j = ( i, σ ).

Let us consider a pure U(1)-FGS | ψ ⟩ , which is fully characterized by its correlation matrix [3, 43]

$$\mathcal { C } _ { j, j ^ { \prime } } = 2 \, \langle \psi | c _ { j } ^ { \dagger } c _ { j ^ { \prime } } | \psi \rangle - \delta _ { j, j ^ { \prime } } \,. \quad \quad ( 1 )$$

C is an N × N Hermitian matrix satisfying C 2 = ✶ N . Diagonalizing the correlation matrix C gives

✶ where M is a unitary matrix. Q is the number of occupied modes in | ψ ⟩ , and the corresponding eigenvectors (first Q rows of the unitary matrix M † ) can be used to write | ψ ⟩ explicitly as a Fermi sea state

$$M ^ { \dagger } \mathcal { C } M & = \begin{pmatrix} 1 _ { Q } & 0 \\ 0 & - 1 _ { N - Q } \end{pmatrix}, \quad & ( 2 ) \\ M \text{ is a unitary matrix. } Q \text{ is the number of oc-}$$

$$| \psi \rangle = \prod _ { q = 1 } ^ { Q } f _ { q } ^ { \dagger } | 0 \rangle \,,$$

where f † q = ∑ N j =1 M † q,j c † j with q = 1 , . . . , Q ( q = Q + 1 , . . . , N ) are occupied (unoccupied) modes in | ψ ⟩ , and | 0 ⟩ is the vacuum, c j | 0 ⟩ = 0 ∀ j .

For our algorithm, we need the Schmidt decomposition of U(1)-FGSs. Let us partition the first N L sites (index j = 1 , . . . , N L ) of the system into the left part (denoted as 'L') and remaining N R = N -N L sites into the right part (denoted as 'R'). Accordingly, the correlation matrix in Eq. (1) is decomposed into four blocks

$$\mathcal { C } = \begin{pmatrix} \mathcal { C } ^ { \text{LL} } & \mathcal { C } ^ { \text{LR} } \\ ( \mathcal { C } ^ { \text{LR} } ) ^ { \dagger } & \mathcal { C } ^ { \text{RR} } \end{pmatrix},$$

where C LL and C RR are N L × N L and N R × N R Hermitian matrices, respectively. Actually, C LL is the correlation matrix of the reduced density operator ρ L for the left subsystem

$$\mathcal { C } _ { j, j ^ { \prime } } ^ { L L } = 2 \, \text{tr} _ { L } ( \rho ^ { L } c _ { j } ^ { \dagger } c _ { j ^ { \prime } } ) - \delta _ { j, j ^ { \prime } } \, \text{ \quad \ \ } ( 5 )$$

with 1 ≤ j, j ′ ≤ N L and 'tr L ' being the trace over degrees of freedom in the left subsystem. Similar definition holds for C RR in the right subsystem.

Substituting Eq. (4) into C 2 = ✶ N gives ( C LL 2 ) + C LR ( C LR ) † = ✶ N L and ( C LR ) † C LR + ( C RR 2 ) = ✶ N -N L . It becomes clear that the two unitary matrices obtained from the singular value decomposition (SVD) of C LR can diagonalize C LL and C RR . For instance, we have

$$U ^ { \dagger } \mathcal { C } ^ { \text{LL} } U = \begin{pmatrix} \Lambda & 0 & 0 \\ 0 & \mathbb { 1 } & 0 \\ 0 & 0 & - \mathbb { 1 } \end{pmatrix} \,, \quad & ( 6 ) & \text{for} \\ \underset { - 1 } { \text{is a diagonal matrix with diagonal entries } } \Lambda _ { \mu } \in \begin{pmatrix} \Lambda & 0 & 0 \\ \Lambda & 0 \\ \Lambda & 0 \end{pmatrix} \,,$$

✶ where Λ is a diagonal matrix with diagonal entries Λ µ ∈ ( -1 1). , The size of the three diagonal blocks in Eq. (6) depends on the entanglement structure of | ψ ⟩ . We shall denote these blocks as 'occupied', 'entangled', and 'unoccupied' (abbreviated as O, E, and U, respectively), for the reason that will become clear soon. The unitary matrix U defines a single-particle basis rotation in the left subsystem

$$d _ { \text{L}, \mu } ^ { \dagger } = \sum _ { j = 1 } ^ { N _ { \text{L} } } U _ { \mu, j } ^ { \dagger } c _ { j } ^ { \dagger } & & ( 7 ) & & \text{(} \mu \\ & & \text{Sc} \\ & & \text{u} \colon & &$$

with µ = 1 , . . . , N L . In this new basis, the three blocks in the diagonalized subsystem correlation matrix [Eq. (6)] translates into

$$\langle \psi | d _ { \mathrm L, \mu } ^ { \dagger } d _ { \mathrm L, \mu } | \psi \rangle = \begin{cases} 1 & \text{if } \mu \in \mathrm O \\ ( 1 + \Lambda _ { \mu } ) / 2 & \text{if } \mu \in \mathrm E \\ 0 & \text{if } \mu \in \mathrm U \end{cases},$$

where µ ∈ O E U distinguish , , three type of modes: d † L ,µ ∈ O ( d † L ,µ ∈ U ) modes are completely occupied (unoccupied) in | ψ ⟩ and have no entanglement with the right subsystem, and d † L ,µ ∈ E are only partially occupied in | ψ ⟩ and hence entangle with certain modes in the right subsystem. In this context, O , E U should not be viewed as , particular sets; they are rather labels for modes (we will use them for the right subsystem, too).

Similarly, the diagonalization of C RR gives

$$V ^ { \dagger } \mathcal { C } ^ { \text{RR} } V = \begin{pmatrix} - \Lambda & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{pmatrix} \,, \quad & ( 9 ) \\ \text{the unitary matrix } V \text{ defines an orthonormal } \text{ the } \text{ua}$$

✶ where the unitary matrix V defines an orthonormal single-particle basis for the right subsystem

$$d _ { \text{R}, \nu } ^ { \dagger } = \sum _ { j = N _ { \text{L} } + 1 } ^ { N } V _ { \nu, j } ^ { \dagger } c _ { j } ^ { \dagger } & & ( 1 0 ) \quad \text{wh} \\ & \text{for}$$

with ν = 1 , . . . , N R . The occupation of these modes in | ψ ⟩ is given by

$$\langle \psi | d _ { R, \nu } ^ { \dagger } d _ { R, \nu } | \psi \rangle = \begin{cases} 1 & \text{if $\nu\in O$} \\ ( 1 - \Lambda _ { \nu } ) / 2 & \text{if $\nu\in E$} \\ 0 & \text{if $\nu\in U$} \end{cases} \quad ( 1 1 ) \quad \text{not $\cosh$} \\ \text{or $p$} \end{cases}$$

Combining that the unitary matrices U † and V SVDdiagonalize C LR [i.e., ( U † C LR V ) µ,ν = ⟨ ψ d | † L ,µ d R ,ν | ψ ⟩ ∝ δ µ,ν ] and that C 2 = ✶ N (valid also in the new basis), we obtain

$$\langle \psi | d _ { \text{L}, \mu } ^ { \dagger } d _ { \text{R}, \nu } | \psi \rangle = \frac { \sqrt { 1 - \Lambda _ { \mu } ^ { 2 } } } { 2 } \delta _ { \mu, \nu } \text{ \quad \ \ } ( 1 2 )$$

for µ, ν ∈ E. Together with Eqs. (8) and (11), we can reconstruct | ψ ⟩ in terms of the orthonormal modes defined in two subsystems as follows:

$$| \psi \rangle & = \prod _ { \mu \in \mathbb { E } } \left [ \sqrt { \frac { 1 + \Lambda _ { \mu } } { 2 } } d _ { \mathbb { L }, \mu } ^ { \dagger } + \sqrt { \frac { 1 - \Lambda _ { \mu } } { 2 } } d _ { \mathbb { R }, \mu } ^ { \dagger } \right ] \\ & \quad \times \prod _ { \mu \in \mathcal { O } } d _ { \mathbb { L }, \mu } ^ { \dagger } \prod _ { \nu \in \mathcal { O } } d _ { \mathbb { R }, \nu } ^ { \dagger } | 0 \rangle \,. \\ \dots & \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots$$

When multiplying out the product for entangled modes ( µ ∈ E), this gives the Schmidt decomposition of | ψ ⟩ .

For later discussions, it will be useful to encode the Schmidt vectors in two U(1)-FGSs living in enlarged Hilbert spaces

$$\dots & \dots \dots \dots \dots \dots \\ & | \phi ^ { L } \rangle = \prod _ { \mu = 1 } ^ { N _ { L } } \left [ \frac { 1 } { \sqrt { 2 } } \left ( d _ { L, \mu } ^ { \dagger } + b _ { N _ { L }, \mu } ^ { \dagger } \right ) \right ] | 0 \rangle _ { p, \nu }, \\ & | \phi ^ { R } \rangle = \prod _ { \nu = 1 } ^ { N _ { R } } \left [ \frac { 1 } { \sqrt { 2 } } \left ( d _ { R, \nu } ^ { \dagger } + b _ { N _ { L + 1 }, \nu } ^ { \dagger } \right ) \right ] | 0 \rangle _ { p, \nu }, \\. & \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 \ 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.$$

where every subsystem orthonormal mode ( d † L ,µ and d † R ,ν ) entangles with one virtual mode ( b † N ,µ L and b † N L +1 ,ν ). | 0 ⟩ p v , is the shared vacuum of physical ('p') and virtual ('v') modes. The virtual modes, indicated by their site indices N L and N L+1 , locate at the 'entangling edge' of the Schmidt decomposition. Introducing a virtual bond state between virtual modes

$$| S \rangle = \prod _ { \mu \in \mathbb { E } } \left [ \sqrt { \frac { 1 - \Lambda _ { \mu } } { 2 } } b _ { N _ { L }, \mu } ^ { \dagger } - \sqrt { \frac { 1 + \Lambda _ { \mu } } { 2 } } b _ { N _ { L } + 1, \mu } ^ { \dagger } \right ] \\ & \times \prod _ { \mu \in \mathbb { U } } b _ { N _ { L }, \mu } ^ { \dagger } \prod _ { \nu \in \mathbb { U } } b _ { N _ { L } + 1, \nu } ^ { \dagger } | 0 \rangle _ { \nu } \,,$$

the original state | ψ ⟩ is restored by contracting the virtual modes (depicted in Fig. 1):

$$\left | \psi \right \rangle = \left \langle S \right | \left ( \left | \phi ^ { L } \right \rangle \otimes \left | \phi ^ { R } \right \rangle \right ) \,, \quad \quad \left ( 1 6 \right )$$

where an unimportant overall factor is ignored. This form of the Schmidt decomposition will be used in Sec. II B. As a minor remark, the tensor product notation in Eq. (16) is not entirely rigorous in mathematical sense, since fermionic systems need a Z 2 -graded tensor product for suitably bookkeeping the fermionic signs. Nevertheless, for our purpose, this notation is intuitive and does not cause sign inconsistency. Below we shall continue with such notations (and sometimes even drop the tensor product symbol).

FIG. 1. Schematic of the Schmidt decomposition of a fermionic Gaussian state | ψ ⟩ via introducing virtual fermionic modes and the associated virtual bond state ⟨ S | .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-b51cdf41-6356-4d80-90e1-cd8e48364099.jpg)

## B. Calculation of the MPS local correlation matrix

In this subsection, we illustrate how to obtain an MPS representation of | ψ ⟩ . The elementary building blocks of MPSs are single-site tensors. The key part of our algorithm is the calculation of a central-site tensor of the MPS. For finite-size systems without translation invariance, MPS tensors at other lattice sites can be obtained by iterating this algorithm. For one-site translation invariant systems, one can gradually increase system sizes such that the central-site tensor is converged; this gives an iMPS (by copying the single-site tensor to every site) in the thermodynamic limit. For systems with multi-site unit cells (e.g., cylinders), iMPSs can be obtained by calculating single-site tensors separately for all sites within one unit cell.

Let us choose a central site with index m m ( ∼ N/ 2). We first follow the procedure in Sec. II A and decompose | ψ ⟩ in the form of Eq. (16):

$$| \psi \rangle = \langle S _ { m, m + 1 } | \left ( | \phi ^ { [ \leq m ] } \rangle \otimes | \phi ^ { [ > m ] } \rangle \right ), \quad ( 1 7 ) \quad \text{and}$$

where ⟨ S m,m +1 | is the virtual bond state, and | ϕ [ ≤ m ] ⟩ and | ϕ [ &gt;m ] ⟩ ) are U(1)-FGSs encoding the Schmidt vectors of the left (sites 1 to m ) and right (sites m + 1 to N ) subsystems, respectively. For the derivations below, we only need | ϕ [ ≤ m ] ⟩ :

$$| \phi ^ { [ \leq m ] } \rangle = \prod _ { \mu = 1 } ^ { m } \left [ \frac { 1 } { \sqrt { 2 } } \left ( d _ { [ \leq m ], \mu } ^ { \dagger } + b _ { m, r, \mu } ^ { \dagger } \right ) \right ] | 0 \rangle _ { p, \nu } \,, \quad ( 1 8 ) \quad \begin{array} { c } \dot { c } _ { m } ^ { \dagger } \, \text{a} \\ b _ { m, r, \mu } ^ { \dagger } \end{array}$$

where d † [ ≤ m,µ ] are orthonormal modes for the left subsystem [see Eq. (7)], and b † m,r,µ are virtual modes on site m (index ' r ' means right side), maximally entangling with d † [ ≤ m,µ ] [see Eq. (14)].

To get a local state at site m (playing the role of an

FIG. 2. Further decomposition of the fermionic Gaussian state | ψ ⟩ to obtain the local state | A [ m ] ⟩ at site m .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-be2521f6-09c8-45ae-96ef-724c964bca57.jpg)

MPS tensor), we decompose | ψ ⟩ further as

$$| \psi \rangle = & \left ( \langle I _ { m - 1, m } | \otimes \langle S _ { m, m + 1 } | \right ) \\ & \times \left ( | \phi ^ { [ < m ] } \rangle \otimes | A ^ { [ m ] } \rangle \otimes | \phi ^ { [ > m ] } \rangle \right ) \quad ( 1 9 )$$

by requiring

$$| \phi ^ { [ \leq m ] } \rangle = \langle I _ { m - 1, m } | \left ( | \phi ^ { [ < m ] } \rangle \otimes | A ^ { [ m ] } \rangle \right ) \quad ( 2 0 )$$

in Eq. (17). Here | ϕ [ &lt;m ] ⟩ is the state encoding the Schmidt vectors of the new left subsystem (including sites 1 to m -1)

$$\underset { \gamma = 1 } { s e } | \phi ^ { [ < m ] } \rangle = \prod _ { \gamma = 1 } ^ { m - 1 } \left [ \frac { 1 } { \sqrt { 2 } } \left ( d _ { [ < m ], \gamma } ^ { \dagger } + b _ { m - 1, r, \gamma } ^ { \dagger } \right ) \right ] | 0 \rangle _ { v } \,, \quad ( 2 1 )$$

and | I m -1 ,m ⟩ is a maximally entangled virtual bond state between sites m -1 and m :

$$\underset { \text{to} } { m } \Big | \Big > = \prod _ { \alpha = 1 } ^ { m - 1 } \left [ \frac { 1 } { \sqrt { 2 } } \left ( b _ { m - 1, r, \alpha } ^ { \dagger } + b _ { m, l, \alpha } ^ { \dagger } \right ) \right ] | 0 \Big > _ { \nu }, \quad ( 2 2 )$$

where b † m,l,α are virtual modes on site m (index ' ' means l left side). The state | A [ m ] ⟩ , appearing in Eq. (20), is a U(1)-FGS at site m , which consists of the physical mode c † m and virtual modes b † m,l,α [ α = 1 , . . . , ( m -1)] and b † m,r,µ ( µ = 1 , . . . , m ). The explicit form of | A [ m ] ⟩ needs to be determined. The graphical representation of the decompositions in Eqs. (17), (19) and (20) is depicted in Fig. 2. To determine | A [ m ] ⟩ , we observe that a result in Ref. [26] [Eq. (A14) thereof] can be used to translate Eq. (20) into a relation of correlation matrices of the U(1)-FGSs involved, from which the correlation matrix of | A [ m ] ⟩ can be calculated. Before using that result,

one should first decompose the correlation matrices into blocks corresponding to certain modes. For instance, the correlation matrix of | A [ m ] ⟩ , denoted by C [ m ] , takes a 3 × 3 block form

$$\mathcal { C } ^ { ( m, l, \alpha ^ { \prime } ) } = ( m ) \begin{pmatrix} \mathcal { C } ^ { [ m ] } _ { 1, 1 } & \mathcal { C } ^ { [ m ] } _ { 1, 2 } & \mathcal { C } ^ { [ m ] } _ { 1, 3 } \\ [ \mathcal { C } ^ { [ m ] } _ { 1, 2 } ] ^ { \dagger } & \mathcal { C } ^ { [ m ] } _ { 2, 2 } & \mathcal { C } ^ { [ m ] } _ { 2, 3 } \\ [ \mathcal { C } ^ { [ m ] } _ { 1, 3 } ] ^ { \dagger } & [ \mathcal { C } ^ { [ m ] } _ { 2, 3 } ] ^ { \dagger } & \mathcal { C } ^ { [ m ] } _ { 3, 3 } \end{pmatrix} \quad ( 2 3 ) \quad \searrow \searrow \mathcal { C }$$

with

$$\text{with} & & \text{$\quad \colon$} \quad \colon \\ \left | \phi ^ { [ \leq m ] } \right \rangle, \\ ( \mathcal { C } ^ { [ m ] } _ { 2, 2 } ) _ { ( m ), ( m ) } & = 2 \langle A ^ { [ m ] } | b ^ { \dagger } _ { m, l, \alpha } b _ { m, l, \alpha ^ { \prime } } | A ^ { [ m ] } \rangle - \delta _ { \alpha, \alpha ^ { \prime } }, \\ ( \mathcal { C } ^ { [ m ] } _ { 2, 2 } ) _ { ( m ), ( m ) } & = 2 \langle A ^ { [ m ] } | c ^ { \dagger } _ { m } c _ { m } | A ^ { [ m ] } \rangle - 1, \\ ( \mathcal { C } ^ { [ m ] } _ { 3, 3 } ) _ { ( m, r, m u ), ( m, r, \mu ^ { \prime } } & = 2 \langle A ^ { [ m ] } | b ^ { \dagger } _ { m, r, \mu } b _ { m, r, \mu ^ { \prime } } | A ^ { [ m ] } \rangle - \delta _ { \mu, \mu ^ { \prime } }, \\ ( \mathcal { C } ^ { [ m ] } _ { 1, 2 } ) _ { ( m, l, \alpha ), ( m ) } & = 2 \langle A ^ { [ m ] } | b ^ { \dagger } _ { m, l, \alpha } c _ { m } | A ^ { [ m ] } \rangle, \\ ( \mathcal { C } ^ { [ m ] } _ { 1, 3 } ) _ { ( m, l, \alpha ), ( m, r, \mu ^ { \prime } } & = 2 \langle A ^ { [ m ] } | b ^ { \dagger } _ { m, l, \alpha } b _ { m, r, \mu ^ { \prime } } | A ^ { [ m ] } \rangle, \\ ( \mathcal { C } ^ { [ m ] } _ { 2, 3 } ) _ { ( m ), ( m, r, \mu ^ { \prime } } & = 2 \langle A ^ { [ m ] } | c ^ { \dagger } _ { m } b _ { m, r, \mu ^ { \prime } } | A ^ { [ m ] } \rangle. \quad \text{(24)} & & \mathcal { C } ^ { [ \leq } \\ \text{For} \text{ ease of reading. } we have indicated the blocks in}$$

For ease of reading, we have indicated the blocks in Eq. (23) with the corresponding indices of the fermionic modes (this will also be used for other correlation matrices below). For rows and columns in Eq. (23) involving the physical mode, ( m ) is an index with dimension one (physical modes are 'spinless'); nevertheless, we keep it to have a uniform format (it will be convenient for extensions to cases where physical modes do have internal degrees of freedom).

The correlation matrix of | ϕ [ &lt;m ] ⟩ [Eq. (21)], denoted by C [ &lt;m ] , has a 2 × 2 block form

$$\mathcal { C } ^ { [ < m ] } = \underset { ( m - 1, r, \gamma ) } { ( j ^ { \prime } < m ) } & \quad ( m - 1, r, \gamma ^ { \prime } ) \\ ( m - 1, r, \gamma ) & \quad \mathcal { C } ^ { [ < m ] } _ { 1, 2 } \\ [ \mathcal { C } ^ { [ < m ] } _ { 1, 2 } ] ^ { \dagger } & 0 \ \end{array} \longrightarrow \, \longrightarrow _ { - }$$

where C [ &lt;m ] 1 2 , has entries ( C [ &lt;m ] 1 2 , ) ( j ) ( , m -1 ,r,γ ′ ) = 2 ⟨ ϕ [ &lt;m ] | c b † j m -1 ,r,γ ′ | ϕ [ &lt;m ] ⟩ . Here j and γ ′ range from 1 to ( m -1). Following Eq. (7), we write the subsystem orthonormal modes as d † [ &lt;m ,γ ] = ∑ m -1 j =1 ( U [ &lt;m ] ) † γ,j c † j , where the unitary matrix U [ &lt;m ] diagonalizes the subsystem correlation matrix from site 1 to ( m -1) [c.f. Eq. (6)]. The definition of | ϕ [ &lt;m ] ⟩ in Eq. (21) gives 2 ⟨ ϕ [ &lt;m ] | d † [ &lt;m ,γ ] b m -1 ,r,γ ′ | ϕ [ &lt;m ] ⟩ = δ γ,γ ′ , which leads to

$$\mathcal { C } _ { 1, 2 } ^ { [ < m ] } = U ^ { [ < m ] }. \quad \text{(26)}$$

Moreover, we obtain 2 ⟨ ϕ [ &lt;m ] | c c † j j ′ | ϕ [ &lt;m ] ⟩ -δ j,j ′ = 2 ⟨ ϕ [ &lt;m ] | b † m -1 ,r,γ b m -1 ,r,γ ′ | ϕ [ &lt;m ] ⟩ -δ γ,γ ′ = 0, which correspond to the vanishing diagonal blocks in Eq. (25).

FIG. 3. Physical and virtual modes forming the vectors | ϕ [ ≤ m ] ⟩ , | ϕ [ &lt;m ] ⟩ , | A [ m ] ⟩ , and | I m -1 ,m ⟩ in Eq. (20).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-40f1986f-0123-4fc9-9c90-60fbeebd5792.jpg)

The correlation matrix of | ϕ [ ≤ m ] ⟩ [Eq. (18)] is given by

$$( j ^ { \prime } < m ) \quad ( m ) \quad ( m, r, \mu ^ { \prime } ) \\ ) _ { \mathcal { C } ^ { [ \leq m ] } } = \quad ( m ) \begin{pmatrix} 0 & 0 & \mathcal { C } _ { 1, 3 } ^ { [ \leq m ] } \\ 0 & 0 & \mathcal { C } _ { 2, 3 } ^ { [ \leq m ] } \\ [ \mathcal { C } _ { 1, 3 } ^ { [ \leq m ] } ] ^ { \dagger } & [ \mathcal { C } _ { 2, 3 } ^ { [ \leq m ] } ] ^ { \dagger } & 0 \end{pmatrix} \quad ( 2 7 ) \\ _ { \mathcal { C } }$$

with ( C [ ≤ m ] 1 3 , ) ( j ) ( , m,r,µ ′ ) = 2 ⟨ ϕ [ ≤ m ] | c b † j m,r,µ ′ | ϕ [ ≤ m ] ⟩ and ( C [ ≤ m ] 2 3 , ) ( m , m,r,µ ) ( ′ ) = 2 ⟨ ϕ [ ≤ m ] | c † m m,r,µ b ′ | ϕ [ ≤ m ] ⟩ , where j ranges from 1 to ( m -1) and µ ′ from 1 to m . The remaining derivation is very similar to that for | ϕ [ &lt;m ] ⟩ above and will not be repeated. If we write the subsystem orthonormal modes as d † [ ≤ m,µ ] = ∑ m j =1 ( U [ ≤ m ] ) † µ,j c † j ≡ ∑ m -1 j =1 ( U [ ≤ m ] 1 ) † µ,j c † j +( U [ ≤ m ] 2 ) † c † m , we obtain

$$\begin{pmatrix} \mathcal { C } _ { 1, 3 } ^ { [ \leq m ] } \\ \mathcal { C } _ { 2, 3 } ^ { [ \leq m ] } \end{pmatrix} = U ^ { [ \leq m ] } = \begin{pmatrix} U _ { 1 } ^ { [ \leq m ] } \\ U _ { 2 } ^ { [ \leq m ] } \end{pmatrix}, \quad \quad ( 2 8 )$$

where U [ ≤ m ] 1 ( U [ ≤ m ] 2 ) is an ( m -1) × m matrix (a 1 × m row vector).

The correlation matrix of the virtual bond state | I m -1 ,m ⟩ [Eq. (22)] is most easily obtained

$$\mathcal { C } ^ { I } = \frac { ( m - 1, r, \gamma ) } { ( m, l, \alpha ) } \binom { 0 } { 1 } _ { m - 1 }, \quad ( 2 9 ) \\ \intertext { \i e r e the non-vanishing off-diagonal block comes from } r = \frac { 1 } { 1 } \vdash \frac { 1 } { 2 } + r \sum \tilde { c } = \tilde { r } _ { m - 2 } \dots \dots \dots }$$

where the non-vanishing off-diagonal block comes from 2 ⟨ I m -1 ,m | b † m -1 ,r,γ b m,l,α ′ | I m -1 ,m ⟩ = δ γ,α ′ . By using Eq. (A14) in Ref. [26], Eq. (20) translates into the fol-

lowing correlation matrix equation:

$$& \text{long correlation matrix equation:} \quad \text{basis} \, \cdot \\ & \quad \left ( \begin{matrix} 0 & 0 & \mathcal { C } _ { 1, 3 } ^ { [ \leq m ] } \\ 0 & 0 & \mathcal { C } _ { 2, 3 } ^ { [ \leq m ] } \\ [ \mathcal { C } _ { 1, 3 } ^ { [ \leq m ] } ] ^ { \dagger } & [ \mathcal { C } _ { 2, 3 } ^ { [ \leq m ] } ] ^ { \dagger } & 0 \end{matrix} \right ) = \left ( \begin{matrix} 0 & 0 & 0 \\ 0 & \mathcal { C } _ { 2, 2 } ^ { [ m ] } & \mathcal { C } _ { 2, 3 } ^ { [ m ] } \\ 0 & [ \mathcal { C } _ { 2, 3 } ^ { [ m ] } ] ^ { \dagger } & \mathcal { C } _ { 3, 3 } ^ { [ m ] } \end{matrix} \right ) \\ & - \left ( \begin{matrix} \mathcal { C } _ { 1, 2 } ^ { [ < m ] } & 0 \\ 0 & [ \mathcal { C } _ { 1, 2 } ^ { [ m ] } ] ^ { \dagger } \\ 0 & [ \mathcal { C } _ { 1, 3 } ^ { [ m ] } ] ^ { \dagger } \end{matrix} \right ) \left ( \begin{matrix} 0 & 1 \\ 1 & \mathcal { C } _ { 1, 1 } ^ { [ m ] } \end{matrix} \right ) ^ { - 1 } \left ( \begin{matrix} [ \mathcal { C } _ { 1, 2 } ^ { [ < m ] } ] ^ { \dagger } & 0 & 0 \\ 0 & \mathcal { C } _ { 1, 2 } ^ { [ m ] } & \mathcal { C } _ { 1, 3 } ^ { [ m ] } \end{matrix} \right ), \quad \text{where} \quad \text{$(30)$} \quad \text{$(n_{0}$} \}, \\ & \text{$(30)$} \quad \text{$instan} \\ & \text{from which the solution for the correlation matrix $\mathcal{C}_{\dots.d.e.}$} \quad \text{$end} \quad a _ { A ^ { [ n _ { 0 } } } }$$

from which the solution for the correlation matrix C [ m ] of the local state | A [ m ] ⟩ reads

$$\mathcal { C } _ { 1, 1 } ^ { [ m ] } & = 0, \quad \mathcal { C } _ { 1, 2 } ^ { [ m ] } = 0, \quad \mathcal { C } _ { 2, 2 } ^ { [ m ] } = 0, \quad \mathcal { C } _ { 3, 3 } ^ { [ m ] } = 0, \quad \text{taker} \\ \mathcal { C } _ { 1, 3 } ^ { [ m ] } & = - [ \mathcal { C } _ { 1, 2 } ^ { [ < m ] } ] ^ { \dagger } C _ { 1, 3 } ^ { [ \leq m ] }, \quad \mathcal { C } _ { 2, 3 } ^ { [ m ] } = C _ { 2, 3 } ^ { [ \leq m ] }. \quad \text{(3 1)} \quad \text{phys} _ { e } ^ { \dagger }.$$

Together with Eqs. (26) and (28), the solution for Eq. (23) is given by

$$\mathcal { C } ^ { [ m ] } = \begin{pmatrix} 0 & 0 & - ( U ^ { [ < m ] } ) ^ { \dagger } U _ { 1 } ^ { [ \leq m ] } \\ 0 & 0 & U _ { 2 } ^ { [ \leq m ] } \\ - ( U _ { 1 } ^ { [ \leq m ] } ) ^ { \dagger } U ^ { [ < m ] } & ( U _ { 2 } ^ { [ \leq m ] } ) ^ { \dagger } & 0 \\ 0 & & & ( 3 2 ) & D \ ( D \\ & & & & \text{total $fc$} \end{pmatrix}$$

Thus, once the unitary matrices diagonalizing subsystem correlation matrices are determined, C [ m ] as the local correlation matrix of the MPS is also obtained.

## C. Calculation of the MPS local tensor

With the local correlation matrix C [ m ] in Eq. (32), we now demonstrate how to obtain a concrete form of the its corresponding Gaussian state | A [ m ] ⟩ and write it explicitly as an MPS local tensor. As the correlation matrix of a U(1)-FGS, C [ m ] can be diagonalized as

$$S ^ { \dagger } \mathcal { C } ^ { [ m ] } S = \begin{pmatrix} 1 _ { Q _ { m } } & 0 \\ 0 & - 1 _ { 2 m - Q _ { m } } \end{pmatrix}, \quad & ( 3 3 ) \quad & \text{MP} \\ \text{e the total number of local modes (including both} \quad & \text{app} \end{pmatrix}$$

✶ where the total number of local modes (including both virtual and physical modes) in | A [ m ] ⟩ is 2 m and S is a unitary matrix. Then, we can write | A [ m ] ⟩ as

$$| A ^ { [ m ] } \rangle = \prod _ { p = 1 } ^ { Q _ { m } } d _ { p } ^ { \dagger } | 0 \rangle _ { p, \nu }, \quad \quad ( 3 4 ) \quad \begin{matrix} \text{on} \\ \text{tio} \\ \text{cor} \\ \text{cui} \end{matrix}$$

where the Q m occupied modes are given by

$$d _ { p } ^ { \dagger } = S _ { p, 0 } ^ { \dagger } c _ { m } ^ { \dagger } + \sum _ { \alpha = 1 } ^ { m - 1 } S _ { p, ( l, \alpha ) } ^ { \dagger } b _ { m, l, \alpha } ^ { \dagger } + \sum _ { \beta = 1 } ^ { m } S _ { p, ( r, \beta ) } ^ { \dagger } b _ { m, r, \beta } ^ { \dagger }. \quad \text{$\substack{$\text{move}$} \\ \text{ter $G$} \\ \text{comp} \\ \text{In $1$} \\ \text{e$} \dots$} \colon$$

The 'column index' 0 in S † p, 0 is used for denoting the local physical mode. If we rewrite Eq. (34) in the Fock basis of physical and virtual modes, it becomes

$$| A ^ { [ m ] } \rangle = & \sum _ { \{ n _ { 0 }, n _ { l, \alpha }, n _ { r, \beta } \} } A ^ { \{ n _ { 0 } \} } _ { \{ n _ { l, \alpha } \}, \{ n _ { r, \beta } \} } \\ & \times ( c ^ { \dagger } _ { m } ) ^ { n _ { 0 } } ( b ^ { \dagger } _ { m, l, \alpha } ) ^ { n _ { l, \alpha } } ( b ^ { \dagger } _ { m, r, \beta } ) ^ { n _ { r, \beta } } | 0 \rangle _ { p, v }, \quad ( 3 6 )$$

where n 0 , n l,α and n r,β are occupation numbers of the local physical, left virtual, and right virtual modes, respectively. Here, we have used a binary list { n 0 } { , n l,α } { , n r,β } to label the fermionic states. For instance, { 0101 } represents a state with only the second and the fourth modes occupied. The coefficient A { n 0 } { n l,α } { , n r,β } in Eq. (36) is the Slater determinant defined with S [ { n 0 } { , n l,α } { , n r,β } ], which is the sub-matrix taken form the unitary matrix S by choosing its first Q m rows and the columns chosen with respect to occupied physical and virtual fermionic modes. The coefficients for virtual states, such as | I m,m +1 ⟩ and | S m,m +1 ⟩ , can be calculated in the same way.

The main bottleneck for this step is that for large systems, converting the local correlation matrix into a local tensor (or calculating all the coefficients) is not manageable, as the bond dimension would grow exponentially. For obtaining an MPS with a feasible bond dimension D ( D = 10 3 ∼ 10 ), 4 we keep only D states of virtual fermions, corresponding to D elements in binary lists { n l,α } and { n r,β } . The selection criterion for the kept elements is that they have dominant contributions to the bipartite entanglement according to the Schmidt decompsition in Eq. (15) (this is equivalent to the SVD truncation [18]). For FGSs with relatively small entanglement (or with small system sizes), this method is already sufficient. However, for FGSs that require large system sizes to obtain a converged MPS local tensor, this method is not efficient enough. In Sec. II D below, we describe an improved method for calculating the MPS local tensor.

It is worth noting that the 'single SVD truncation' scheme described above (i.e., without further DMRG sweeps) is not guaranteed to provide the most accurate MPS approximation for a fixed bond dimension D , as we simultaneously truncate all the virtual bonds. Nevertheless, we observe that this scheme already provides a good approximation as long as the bond dimension D is not too small.

As a further note, we mainly focus in this manuscript on Gutzwiller projected U(1)-FGSs as trial wave functions for spin systems, where a local fermion-number constraint restores the spin Hilbert space (e.g., singly occupied two-component fermions at each site representing spin-1/2). With a fixed number of local physical modes implemented by the Gutzwiller projection, virtual states at each sites also have a fixed particle number after Gutzwiller projection. This can obviate the need to compute unnecessary MPS local tensor entries.

In the construction of a fermionic tensor network, the fermionic sign should be carefully managed. For practical calculations, one needs to contract the virtual modes in Eq. (36) with those in the virtual bond states to ob-

tain an ordinary MPS, where fermionic signs arise due to moving and contracting the virtual modes. However, these signs can be absorbed into the definition of the MPS local tensor (see Ref. [44]).

## D. Mode decimation

In Sec. II A, II B, and II C, we have described a complete conversion algorithm, as we start from a U(1)-FGS | ψ ⟩ and arrive at an explicit form of the MPS local tensor. There, the elements of the MPS local tensor are given by the Slater determinants of some sub-matrix of S [Eqs. (33) and (35)], so the computational complexity scales as O N ( 3 ), where N is the system size. This scaling is polynomial but may still prevent us from simulating large systems. As reaching large system sizes is crucial for obtaining iMPS, the present subsection provides an improved method to compute the MPS local tensor.

If we follow the SVD truncation scheme in Sec. II C and restrict the number of virtual states in Eq. (36) to D (by choosing elements in { n l,α } and { n r,β } with dominant entanglement support), some virtual modes ( b † m,l,α and b † m,r,β ) are occupied or unoccupied [see Eq. (15); for these modes Λ µ are close to ± 1]. We estimate that only around log ( 2 D ) modes are 'active' (with Λ µ not close to ± 1) in these D virtual states. With this observation in mind, we find that it is possible to define another local U(1)-FGS | ˜ A [ m ] ⟩ which has less virtual modes (and thus a smaller correlation matrix ˜ C [ m ] ) but gives exactly the same result when calculating the MPS local tensor with bond dimension D . As the size of ˜ C [ m ] is of order log 2 ( D ), the cost of generating a single MPS local tensor element can be reduced from O N ( 3 ) to O ((log 2 ( D )) 3 ). This is also more efficient than O N ( 2 log 2 ( D )), where the latter requires a fast update strategy for Slater determinants (see Ref. [18]).

In the following, we describe how to construct | ˜ A [ m ] ⟩ by mode decimation and generate the MPS local tensor from it. Hereafter, we use a tilde hat for quantities where the mode decimation has been performed.

Due to the SVD truncation, some entangled modes in the Schmidt decomposition become occupied (unoccupied) modes. To simplify the notation, we transfer this kind of entangled modes from set E to set O (U) and denote the resulting set as E, O and U. ˜ ˜ ˜ Then, the local Gaussian state under mode decimation can be expressed as

$$| \tilde { A } ^ { [ m ] } \rangle = ( \langle I ^ { ( m, l ) } | \otimes \langle I ^ { ( m, r ) } | ) ( | P ^ { ( m, l ) } \rangle \otimes | A ^ { [ m ] } \rangle \otimes | P ^ { ( m, r ) } \rangle ) \quad \text{the co}$$

FIG. 4. Schematic of the Gaussian states (a) | ˜ A [ m ] ⟩ and (b) ⟨ S m,m +1 | . The solid dots represent fermionic modes. The modes connected by a spiral line belong to the same Gaussian state. Different dots or modes in a single dashed rectangular are contracted.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-d7361ee6-0c3a-4c24-892c-4d1699764636.jpg)

with

$$C _ { \alpha } & \text{ with } \\ | P ^ { ( m, l ) } \rangle = \prod _ { \beta \in \tilde { U } } a _ { m, l, \beta } ^ { \dagger } \prod _ { \alpha \in \tilde { E } } ( f _ { m, l, \alpha } ^ { \dagger } - a _ { m, l, \alpha } ^ { \dagger } ) | 0 \rangle, \\ \mathcal { r } _ { \alpha } & | I ^ { ( m, l ) } \rangle = \prod _ { \alpha = 1 } ^ { m - 1 } ( a _ { m, l, \alpha } ^ { \dagger } + b _ { m, l, \alpha } ^ { \dagger } ) | 0 \rangle, \\ \mathfrak { i } _ { \alpha } & | I ^ { ( m, r ) } \rangle = \prod _ { \alpha = 1 } ^ { m } ( b _ { m, r, \alpha } ^ { \dagger } + a _ { m, r, \alpha } ^ { \dagger } ) | 0 \rangle, \\ | P ^ { ( m, r ) } \rangle = \prod _ { \alpha \in \tilde { E } } ( a _ { m, r, \alpha } ^ { \dagger } - f _ { m, r, \alpha } ^ { \dagger } ) \prod _ { \beta \in \tilde { U } } a _ { m, r, \beta } ^ { \dagger } | 0 \rangle. \quad ( 3 8 ) \\ \mathfrak { h } ), & \text{ Here we have introduced two new types of virtual modes,}$$

Here we have introduced two new types of virtual modes, denoted by f † and a † , their positions are depicted in Fig. 4(a). The projectors P ( m,r ) and P ( m,l ) can help fix or freeze modes belonging to U or O. ˜ ˜ As f † µ ∈ ˜ , E the number of f modes can be smaller than that of b modes, so | ˜ A [ m ] ⟩ has less virtual modes.

To calculate tensor elements using | ˜ A [ m ] ⟩ , we should also perform mode decimation on the virtual states. For instance, if we start with the virtual state { 101100 } and find that the first and the last modes should be fixed, then the virtual state for | ˜ A [ m ] ⟩ should be { 0110 } . Then | P ( m,r ) ⟩ together with ⟨ I ( m,r ) | can be treated as a 'signal processor'. If the input on the right side is { 0110 } (configuration of f † m,r,α ), the output on the left side is { 101100 } (configuration of b † m,r,α ). Following this procedure, one can easily verify that calculating the local tensor using | ˜ A [ m ] ⟩ or | A [ m ] ⟩ gives the same result.

The block structure of ˜ C [ m ] is much more complicated than that of C [ m ] as some blocks no longer vanish due to mode decimation. To derive ˜ C [ m ] , we should first obtain the correlation matrices for those truncation projectors. The correlation matrix of | P ( m,r ) ⟩ can be easily calcu-

lated:

where

$$\mathcal { C } ^ { ( m, r ) } = \frac { a } { f } \begin{pmatrix} a & f \\ C _ { 1, 1 } ^ { ( m, r ) } & C _ { 1, 2 } ^ { ( m, r ) } \\ [ C _ { 1, 2 } ^ { ( m, r ) } ] ^ { \dagger } & C _ { 2, 2 } ^ { ( m, r ) } \end{pmatrix}, \quad \text{ (3 9)}$$

$$\mathcal { C } _ { 2, 2 } ^ { ( m, r ) } = 0,$$

$$[ C _ { 1, 2 } ^ { ( m, r ) } ] _ { \alpha, \beta } = - \delta _ { \alpha \beta },$$

$$[ C ^ { ( m, r ) } _ { 1, 1 } ] _ { \alpha, \beta } = \delta _ { \alpha \beta } \begin{cases} - 1, & \alpha \in \tilde { \mathbf U } \\ 1, & \alpha \in \tilde { \mathbf O } \, & \quad ( 4 2 ) & \text{Con} \\ 0, & \alpha \in \tilde { \mathbf E } & \quad \text{ each} \\ \text{a 3} \end{cases}$$

and the correlation matrix of | P ( m,l ) ⟩ reads

$$C ^ { ( m, l ) } = \frac { f } { a } \begin{pmatrix} f & a \\ C _ { 1, 1 } ^ { ( m, l ) } & C _ { 1, 2 } ^ { ( m, l ) } \\ [ C _ { 1, 2 } ^ { ( m, l ) } ] ^ { \dagger } & C _ { 2, 2 } ^ { ( m, l ) } \end{pmatrix},$$

where

$$\Re ^ { ( m, l ) } _ { 1, 1 } & = 0, \\ [ C ^ { ( m, l ) } _ { 1, 2 } ] _ { \alpha, \beta } & = - \delta _ { \alpha \beta }, \\ [ C ^ { ( m, l ) } _ { 2, 2 } ] _ { \alpha, \beta } & = \delta _ { \alpha \beta } \begin{cases} 1, & \alpha \in \tilde { \mathbb { U } } \\ - 1, & \alpha \in \tilde { O } \\ 0, & \alpha \in \tilde { E } \end{cases}. \quad \ \ ( 4 4 ) \\ \dots \quad - \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \alpha \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Lambda \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \Delta \quad \epsilon \quad \quad \quad$$

Combining Eq. (A14) in Ref. [26] and Eq. (37), we obtain each block of the correlation matrix ˜ C [ m ] [which also takes a 3 × 3 block structure as C [ m ] in Eq. (23)]:

$$\tilde { c } ^ { [ m ] } _ { 1, 1 } & = \tilde { c } ^ { ( m, t ) } _ { 1, 2 } \{ 1 - \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1, 2 } ^ { ( m, t ) } ] \tilde { t } } _ { 1, 2 } ^ { ( m, t ) } ] \tilde { t } } _ { 1, 2 } ^ { ( m, t ) } ] \tilde { t } } _ { 1, 2 } ^ { \dagger }, \\ \tilde { c } ^ { [ m ] } _ { 1, 2 } & = \tilde { c } ^ { ( m, t ) } _ { 1, 2 } \{ 1 - \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1, 3 } ^ { \dagger } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1, 3 } ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1, 3 } ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1 } \} ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } [ \tilde { c } ^ { ( m ] } _ { 1, 3 } ] \tilde { t } } _ { 1, 3 } ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { [ m ] } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 1 } \} ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 2, 2 } \} ^ { - 1 } \tilde { c } ^ { ( m, t ) } _ { 1, 3 } \tilde { c } ^ { ( m, r ) } _ { 1, 2 } \} ^ { - 1 }, \\ \tilde { c } ^ { [ m ] } _ { 1, 2 } & = \tilde { c } ^ { ( m ( m,$$

The mode decimation in the virtual bond states can be performed in a similar way. Here we only discuss ⟨ ˜ S m,m +1 | , and the calculation for ⟨ ˜ I m,m +1 | is similar. To make sure that ⟨ ˜ S m,m +1 | gives the same result as ⟨ S m,m +1 | when generating the local matrix on the virtual bond, we require that ⟨ S m,m +1 | = ( ⟨ I ( m,r ) | ⊗ ⟨ ˜ S m,m +1 | ⊗ ⟨ I ( m +1 ) ,l | ) ( | P ( m,r ) ⟩ ⊗ | P ( m +1 ) ,l ⟩ ); see Fig. 4(b). It turns out that the virtual bond state after mode decimation is given by

In the language of MPS, the MPS representation of any superposition of anyon eigenstates is non-injective and the leading eigenvalues of its transfer matrix are degenerate. The same is true for systems with spontaneous symmetry breaking.

$$| \tilde { S } _ { m, m + 1 } \rangle = \prod _ { \alpha \in \tilde { E } } \left [ \sqrt { \frac { 1 - \Lambda _ { \alpha } } { 2 } } f _ { m, r, \alpha } ^ { \dagger } - \sqrt { \frac { 1 + \Lambda _ { \alpha } } { 2 } } f _ { m + 1, l, \alpha } ^ { \dagger } \right ] | 0 \rangle. \begin{array} { c } \text{the iMP} \mathfrak { s } \\ \text{trix for t} \\ \text{state} \end{array} } \text{$the an} \text{$me$} \text{$a$}$$

## E. Filtration of anyon eigenbasis

Many Gutzwiller projected FGSs in two dimensions have intrinsic topological order, and this is also so for the two examples which we will consider subsequently. Generally, the ground states for topologically ordered systems on a cylinder are the superposition of anyon eigenstates.

If the Gutzwiller projected FGSs have already been converted to iMPSs, we can filter out the degenerate ground states by using the fixed points of the iMPS transfer matrix. Let us denote the transfer matrix by T . If the iMPS has a multi-site unit cell, T is the transfer matrix for the whole unit cell. Starting with some 'product state' Ans¨ atze E trial , we apply the transfer matrix on them and iterate many times to search for distinct fixed points (as left eigenvectors of T ), lim n →∞ E trial T n → E ζ with ζ = 1 , . . . , d deg , where d deg is the number of distinct fixed points (equivalently, the number of degenerate ground states). These fixed points can be normalized such that E E ζ † ζ ′ = δ ζ,ζ ′ .

Once all distinct fixed points have been found, we use the following steps to obtain the iMPS local tensor for each ground state: (i) Diagonalize each fixed point E ζ , E ζ = U S U ζ ζ † ζ with S ζ ≥ 0. (ii) Truncate S ζ (remove its

FIG. 5. (a) Diagonalizing the fixed points E ζ and truncating its vanishing eigenvalues to obtain the isometries V ζ . (b) Use the isometries to obtain the iMPS local tensor for each degenerate ground state.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-360697d6-ce41-41c5-bf34-152f14bd3352.jpg)

zero eigenvalues) to obtain E ζ = V S V ζ ˜ ζ † ζ , where ˜ S ζ &gt; 0 encodes positive eigenvalues. V ζ are isometries satisfying V V ζ † ζ = ✶ and V † ζ V ζ = P ζ , where P ζ are orthogonal projectors satisfying a completeness relation ∑ d deg ζ =1 P ζ = ✶ . (iii) Use the isometries to project onto the virtual space. If we denote the iMPS local matrix as A (physical index omitted; it may correspond to a block of local tensors if one deals with a multi-site unit cell), the iMPS local tensor for each degenerate ground state is given by ˜ A ζ = V † ζ AV ζ . The full procedure is depicted in Fig. 5.

## III. RESULTS

This section presents our numerical results in two chiral spin liquid models. The first (second) one has the same topological order as the Laughlin (Moore-Read) state of bosons at filling factor 1 / 2 (1) [45, 46]. In both cases, the anyon eigenbasis is constructed and their entanglement spectra are computed [47].

## A. S = 1 2 / Abelian chiral spin liquid

The first model is defined on the square lattice with each site occupied by one spin-1/2 [15, 41, 48]. There are N x ( N y ) sites along the x ( y ) direction. The spin1/2 operators are expressed in terms of fermionic parton operators as

$$S _ { j } ^ { a } = \frac { 1 } { 2 } \sum _ { \alpha, \beta = \uparrow, \downarrow } c _ { j, \alpha } ^ { \dagger } \sigma _ { \alpha, \beta } ^ { a } c _ { j, \beta }, \quad \quad ( 4 7 ) \quad \text{$any$} \\ \text{$true$} \\ \text{$the$}$$

where j is the lattice site index, σ a ( a = x, y, z ) are the Pauli matrices, and c j,α ( c † j,α ) is the fermionic annihilation (creation) operator. A single-occupancy constraint ∑ α = ↑ ↓ , c † j,α c j,α = 1 ensures that the physical states at each site are spin-1/2 states, represented by | ↑ ⟩ j = c † j, ↑ | vac ⟩ and | ↓ ⟩ j = c † j, ↓ | vac , where ⟩ | vac ⟩ is the vacuum of fermionic partons.

The hopping Hamiltonian for the fermionic partons is defined as [48]

$$H = \sum _ { \langle j k \rangle, \alpha } t _ { j k } c _ { j, \alpha } ^ { \dagger } c _ { k, \alpha } + \sum _ { \langle \langle j k \rangle \rangle, \alpha } i \Delta _ { j k } c _ { j, \alpha } ^ { \dagger } c _ { k, \sigma } \,, \quad ( 4 8 )$$

where ⟨ jk ⟩ and ⟨⟨ jk ⟩⟩ denote nearest and next-nearest neighbors, respectively. The hopping parameters have absolute values | t jk | = 1 0 and . | ∆ jk | = 0 5 and their . signs are displayed in Fig. 7(a). The lower band with Chern number 1 for each spin component is fully occupied. This Fermi sea is acted on by Gutzwiller projector P G to generate the desired chiral spin liquid.

To characterize the SU(2) 1 topological order of this model, it is convenient to employ cylinders that is open along the x direction and periodic along the y direction with boundary twist angle θ y . As already demonstrated in previous works, the anyon eigenbasis can be constructed easily when the parton Hamiltonian has exact zero modes localized at the two ends of the cylinders [15, 41]. If N y ∈ 4 N + and θ y = 0 or N y ∈ 4 N + -2 and θ y = π , there are four zero modes on either side of the cylinder (two for each spin component of the partons). The anyon eigenstate in the identity sector is obtained when the spin-up and spin-down zero modes on one end of the cylinder are occupied. The anyon eigenstate in the semion sector is obtained when the spin-up zero mode on one end and spin-down zero mode on the other end are occupied.

̸

One significant advance enabled by our method is that there is no need to manually select the edge modes. Instead, with the iMPS in hand, it is possible to filter out the anyon eigenbasis by analyzing the fixed points of the transfer matrix. When θ y = 0 , π , there is no exact zero mode. Tuning θ y from 0 to π , we numerically observed that only the identity sector is obtained. To get the semion sector, one needs to slightly couple the edge modes on each side. Two sectors get mixed in this case and the filtration scheme works. This procedure has been carried out in various systems and we focus on the case with N x = 64 and N y = 10. The lattice is mapped to a chain with a snake shape to accommodate the MPS. When the system is sufficiently long along the x direction, the MPS tensors at the center can be taken out to serve as unit cells of iMPS. There are 20 sites in each unit cell because the hopping Hamiltonian has period 2 along the x direction. We have constructed the iMPS representation with bond dimension D = 800 for both anyon eigenstates. To compute the entanglement spectrum, the system is divided into left and right halves so the good quantum numbers are K L y (the total momentum along y ) and S L z (the azimuthal spin). The edge physics of the 1 / 2 Laughlin state is captured by the chiral SU(2) 1 Wess-Zumino-Witten (WZW) model. This is a chiral conformal field theory (CFT) with two primary fields corresponding to affine Dynkin labels [0 1] ,

## (a) identity sector

FIG. 6. Λ p level for spin-up and spin-down orbitals in (a) the identity sector and (b) the semion sector of the spin-1/2 Laughlin chiral spin liquid. Only eight single-particle levels are plotted for depiction and those with Λ p = 0 are marked in orange. The Schmidt vectors with the largest and second largest singular value, denoted by λ 0 and λ 1 , respectively, are given and the red crosses represent the occupied virtual mode. In (b), we use arrows to indicate the eight degenerate virtual states at λ 1 (only one of the two modes connected by an arrow is occupied in each virtual state, so we have 2 3 virtual states).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-765a9ec8-0236-4419-b02b-96ccf364673a.jpg)

| Λ p                   |       |                |      |
|-----------------------|-------|----------------|------|
| spin                  |       |                |      |
| edge                  | left  |                |      |
|                       |       | left left left | left |
|                       | left  | left λ         | left |
| (b) semion sector Λ p | λ     |                |      |
|                       | 0     | 1              |      |
| spin                  |       |                |      |
|                       | right |                |      |
| edge                  | left  |                |      |
|                       |       | left right     |      |
|                       | λ 0   | 1              |      |
|                       |       | left right λ   |      |

and [1 0] [49]. , The energy spectrum of this model is well known. As shown in Figs. 7(b) and (c), the counting of entanglement levels agrees with the CFT predictions.

It is helpful to inspect the spin symmetry of the anyon eigenstates. Both the parton Hamiltonian and Gutzwiller projection respect SU(2) spin rotation symmetry. In the Schmidt decomposition of an SU(2)-symmetric MPS, the Schmidt vectors associated with the same virtual block must have the same singular values. For the identity sector, this symmetry can be preserved by the anyon eigenstate as shown in Fig. 6(a). The largest singular value λ 0 is unique when a cut is made at the center, so it belongs to the spin-0 sector. The second-largest singular value λ 1 is four-fold degenerate, which means that the associated Schmidt vectors arise from tensor product of two spin-

1 / 2 and their virtual block is 1 2 ⊗ 1 2 = 0 ⊕ 1. The other virtual blocks can be understood similarly. On the contrary, the SU(2) symmetry is broken done to U(1) in the semion sector. This is already evident in Fig. 7(c) as the entanglement levels with S L z = ± 1 are not identical. If we construct a spin-singlet ansatz for the semion state, the SU(2) symmetry can be restored [15]. The largest singular value λ 0 is two-fold degenerate so its virtual block has spin-1 / 2. The second largest singular value λ 1 corresponds to the virtual block 1 2 ⊗ 1 2 ⊗ 1 2 = 1 2 ⊕ 1 2 ⊕ 3 2 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-4d56e294-09c5-4616-9947-3c7dbfbd7335.jpg)

(c) [0,1] sector

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-2be564ed-4feb-4fcb-bd59-6088cc6b1bcd.jpg)

(b) [1,0] sector

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-b7d98004-9472-481a-8bdc-8b09740ee690.jpg)

ξ

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-23a77f8f-eda7-4baf-80bf-660c2259fd3f.jpg)

(e) [1,1] sector

(f) [0,2] sector

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-58e9c4f0-a9e8-405c-8fbb-731bb5b1d2fb.jpg)

FIG. 7. (a) Schematics of signs of t i,j and ∆ i,j in parton Hamiltonian H . The signs of t i,j are indicated by the ± along the associated link. The signs of ∆ i,j are negative (positive) along (against) the arrows on the colored lines. Entanglement spectrum of (b) the identity sector and (c) the semion sector of the SU (2) 1 chiral spin liquid. Entanglement spectrum of (d) the identity sector, (e) the Ising anyon sector, and (f) the fermion sector of the SU (2) 2 chiral spin liquid.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01155v1-figure-8d0f50b5-ecb8-4c48-a543-8d8fce15f6d9.jpg)

## B. S = 1 non-Abelian chiral spin liquid

terms of fermionic partons as follows:

The (bosonic) Moore-Read state is a canonical example of non-Abelian topological order associated with the chiral SU(2) 2 WZW model [46]. Its three anyon eigenstates are labeled as [2,0] (identity sector), [1,1] (Ising anyon sector), and [0,2] (fermion sector). In the quantum Hall contexts, the bosonic Moore-Read state can be constructed using two copies of the bosonic Laughlin state with proper symmetrization. It has been proposed that the same topological order can be realized in spin-1 models [50-57] as well as the Yao-Lee model [58].

We use the fermionic parton approach to construct a spin-1 chiral spin liquid state which has the same topological order as the bosonic Moore-Read state. To this end, we adopt a two-orbital construction [59] and write the spin-1 operators as

$$S _ { j } ^ { a } = \frac { 1 } { 2 } \sum _ { \tau = 1, 2 } \sum _ { \alpha, \beta = \uparrow, \downarrow } c _ { j, \tau, \alpha } ^ { \dagger } \sigma _ { \alpha, \beta } ^ { a } c _ { j, \tau, \beta }, \quad ( 4 9 ) \quad \text{eige}.$$

where τ = 1 2 stands for two different orbitals , at the same site. The local constraint for restoring the spin-1 Hilbert space is imposed by requiring ⃗ S 2 j = 2, so that each orbital must occupy exactly one fermionic parton (carrying spin-1/2) and the two fermionic partons must form a spin-1 triplet. The three spin-1 states are expressed in

$$\begin{array} { c } \cdot \\ | 1 _ { j } \rangle = c _ { j, 1, \uparrow } ^ { \dagger } c _ { j, 2, \uparrow } ^ { \dagger } | \text{vac} \rangle, \\ | - 1 _ { j } \rangle = c _ { j, 1, \downarrow } ^ { \dagger } c _ { j, 2, \downarrow } ^ { \dagger } | \text{vac} \rangle, \\ | \ 0 _ { j } \rangle = \frac { 1 } { \sqrt { 2 } } ( c _ { j, 1, \downarrow } ^ { \dagger } c _ { j, 2, \uparrow } ^ { \dagger } + c _ { j, 1, \uparrow } ^ { \dagger } c _ { j, 2, \downarrow } ^ { \dagger } ) | \text{vac} \rangle. \quad ( 50 ) \\ \cdot \quad \cdot \cdot \cdot \cdot \cdot \end{array} \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \end{array}$$

Accordingly, the Gutzwiller projector P G is defined to remove states other than those in Eq. (50).

The parton Hamiltonian of this system is identical to Eq. (48), except that the fermionic partons have an additional orbital index. In the ground state, the lower band with Chern number 1 for each spin and orbital component is fully occupied. The Gutzwiller projector is applied on top of the parton Fermi seas to generate a spin1 wave function which is expected to exhibit the SU(2) 2 chiral topological order. However, unlike the SU(2) 1 case in Sec. II A, it is a prior unclear how to define the anyon eigenbasis in the current spin-1 model.

Using our filtration scheme, three anyon eigenstates are constructed using bond dimension D = 1600 and their entanglement spectra are displayed in Fig. 7(d)-(f). The counting of entanglement levels agrees with theoretical predictions based on the chiral SU(2) 2 WZW model. This success is a highly nontrivial corroboration of our method. If the parton hopping Hamiltonian indeed has exact zero modes, we can write down the anyon eigenbasis explicitly. To begin with, all parton modes with

negative energy are occupied to generate the state | Φ . ⟩ The zero modes are created by the operators d † s,τ,α , where s = L, R indicates the left or right edge of the cylinder, τ = 1 2 stands for the two orbitals, and , α = ↑ ↓ , represents spin-up or spin-down partons. The anyon eigenbasis is given by

$$| \psi _ { [ 2, 0 ] } \rangle & = P _ { \text{G} } \, d _ { L, 1, \uparrow } ^ { \dagger } d _ { L, 1, \downarrow } ^ { \dagger } d _ { L, 2, \uparrow } ^ { \dagger } d _ { L, 2, \downarrow } ^ { \dagger } | \Phi \rangle, \quad \text{and} \\ | \psi _ { [ 1, 1 ] } \rangle & = P _ { \text{G} } \, d _ { L, 1, \uparrow } ^ { \dagger } d _ { L, 1, \downarrow } ^ { \dagger } d _ { L, 2, \uparrow } ^ { \dagger } d _ { R, 2, \downarrow } ^ { \dagger } | \Phi \rangle, \quad \text{and} \\ | \psi _ { [ 0, 2 ] } \rangle & = P _ { \text{G} } \, d _ { L, 1, \uparrow } ^ { \dagger } d _ { R, 1, \downarrow } ^ { \dagger } d _ { L, 2, \uparrow } ^ { \dagger } d _ { R, 2, \downarrow } ^ { \dagger } | \Phi \rangle \,. \quad \text{(51)} \quad \text{trith} \\ \text{way}.$$

This construction is not unique. For instance, the state in the identity sector, ψ [2 0] , , can also be expressed as P G d † L, 1 , ↑ d † L, 1 , ↓ d † R, , 2 ↑ d † R, , 2 ↓ | Φ . ⟩

The results in Eq. (51) reveal how to relate the anyon eigenbasis of the bosonic Moore-Read states to that of the bosonic Laughlin states. If we interpret the two different orbitals as different 'copies', the first line of Eq. (51) shows that the identity sector of the Moore-Read state is obtained using two copies of the Laughlin state in the identity sector. If one (both) copy of the Laughlin state is replaced by its semion sector, we would obtain the Ising anyon (fermion) sector of the Moore-Read state. This idea may also be utilized to generate other SU(2) k nonAbelian chiral spin liquid ( k &gt; 2). In particular, the k = 3 case possesses Fibonacci anyons but there are not many studies in the spin liquid context [60].

## IV. SUMMARY AND DISCUSSIONS

In summary, we have demonstrated that the (infinite) MPSrepresentation of a particle-number-conserving FGS can be constructed efficiently using our algorithm. Thanks to the extensive usage of correlation matrices in the present work, computational cost is significantly reduced compared to previous studies so much larger sys-

- [1] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Berlin, 1980).
- [2] V. Bach, E. H. Lieb, and J. P. Solovej, J. Stat. Phys. 76 , 3 (1994).
- [3] S. Bravyi, Quantum Inf. and Comp. 5 , 216 (2005).
- [4] C. V. Kraus and J. I. Cirac, New J. Phys. 12 , 113004 (2010).
- [5] C. Gros, Ann. Phys. 189 , 53 (1989).
- [6] T. Xiang, Density Matrix and Tensor Network Renormalization (Cambridge University Press, 2023).
- [7] J. I. Cirac, D. P´ erez-Garc´ ıa, N. Schuch, and F. Verstraete, Rev. Mod. Phys. 93 , 045003 (2021).
- [8] F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57 , 143 (2008).
- [9] J. I. Cirac and F. Verstraete, J. Phys. A 42 , 504004 (2009).

tem sizes can be reached. For a given MPS that describes a topologically ordered system, the anyon eigenbasis can be filtered out in a straightforward manner by analyzing the fixed points of the transfer matrix. This helps to eliminate the laborious process of adjusting the boundary conditions to search for proper edge mode occupations. The algorithm is applied to study two chiral spin liquids and the entanglement spectra of their anyon eigenbasis agree with the CFT predictions.

It would be very interesting to further test our algorithm. If the iMPSs on the cylinder are twisted in certain ways, the modular matrix could be computed to provide a quite comprehensive characterization [40-42]. Quantum phase transitions of topological states are routinely investigated in the parton framework and numerical verification of many theoretical predictions could be facilitated by our method. When a microscopic model is studied using DMRG, a wisely chosen initial state may substantial speedup the calculation. For a 2D system with topological order, feeding the program with suitable parton wave functions rather than random states could help us to find the anyon eigenbasis more easily.

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

## ACKNOWLEDGMENTS

T.L. thanks Tong-Zhou Zhao for helpful discussions. H.H.T. is grateful to Jan von Delft, Hui-Ke Jin, JhengWei Li, Rong-Yang Sun, Lei Wang, Qi Yang, Yi Zhou for stimulating discussions and fruitful collaborations on related topics. T.L. and T.X. are supported by the National Natural Science Foundation of China under Grant No. 12488201. Y.H.W. is supported by the National Natural Science Foundation of China under Grant No. 12174130.

- [10] U. Schollw¨ck, Ann. Phys. o 326 , 96 (2011).
- [11] R. Or´s, Ann. Phys. u 349 , 117 (2014).
- [12] C.-P. Chou, F. Pollmann, and T.-K. Lee, Phys. Rev. B 86 , 041105 (2012).
- [13] H.-H. Zhao, K. Ido, S. Morita, and M. Imada, Phys. Rev. B 96 , 085103 (2017).
- [14] M. T. Fishman and S. R. White, Phys. Rev. B 92 , 075132 (2015).
- [15] Y.-H. Wu, L. Wang, and H.-H. Tu, Phys. Rev. Lett. 124 , 246401 (2020).
- [16] H.-K. Jin, H.-H. Tu, and Y. Zhou, Phys. Rev. B 101 , 165135 (2020).
- [17] A. M. Aghaei, B. Bauer, K. Shtengel, and R. V. Mishmash, arXiv:2009.12435.
- [18] G. Petrica, B.-X. Zheng, G. K.-L. Chan, and B. K. Clark, Phys. Rev. B 103 , 125161 (2021).

- [19] H.-K. Jin, R.-Y. Sun, Y. Zhou, and H.-H. Tu, Phys. Rev. B 105 , L081101 (2022).
- [20] S. R. White, Phys. Rev. Lett. 69 , 2863 (1992).
- [21] H.-K. Jin, H.-H. Tu, and Y. Zhou, Phys. Rev. B 104 , L020409 (2021).
- [22] J.-Y. Chen, J.-W. Li, P. Nataf, S. Capponi, M. Mambrini, K. Totsuka, H.-H. Tu, A. Weichselbaum, J. von Delft, and D. Poilblanc, Phys. Rev. B 104 , 235104 (2021).
- [23] H.-K. Jin, W. M. H. Natori, and J. Knolle, Phys. Rev. B 107 , L180401 (2023).
- [24] R.-Y. Sun, H.-K. Jin, H.-H. Tu, and Y. Zhou, npj Quantum Mater. 9 , 16 (2024).
- [25] Q. Mortier, N. Schuch, F. Verstraete, and J. Haegeman, Phys. Rev. Lett. 129 , 206401 (2022).
- [26] J.-W. Li, J. von Delft, and H.-H. Tu, Phys. Rev. B 107 , 085148 (2023).
- [27] Q. Yang, X.-Y. Zhang, H.-J. Liao, H.-H. Tu, and L. Wang, Phys. Rev. B 107 , 125128 (2023).
- [28] Y. He, K. Li, Y. Zhang, and H. C. Po, Phys. Rev. Res. 6 , L022016 (2024).
- [29] H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101 , 090603 (2008).
- [30] J. Jordan, R. Or´s, G. Vidal, F. Verstraete, u and J. I. Cirac, Phys. Rev. Lett. 101 , 250602 (2008).
- [31] V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. B 79 , 195119 (2009).
- [32] M. Lubasch, J. I. Cirac, and M.-C. Ba˜uls, Phys. Rev. n B 90 , 064425 (2014).
- [33] P. Corboz, Phys. Rev. B 94 , 035133 (2016).
- [34] L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Phys. Rev. B 94 , 155123 (2016).
- [35] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Phys. Rev. X 9 , 031041 (2019).
- [36] M. Scheb and R. M. Noack, Phys. Rev. B 107 , 165112 (2023).
- [37] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Inf. and Comp. 7 , 401 (2007).
- [38] E. Keski-Vakkuri and X.-G. Wen, Int. J. Mod. Phys. B 07 , 4227 (1993).
- [39] Y. Zhang, T. Grover, A. Turner, M. Oshikawa, and A. Vishwanath, Phys. Rev. B 85 , 235151 (2012).
- [40] L. Cincio and G. Vidal, Phys. Rev. Lett. 110 , 067208 (2013).
- [41] H.-H. Tu, Y. Zhang, and X.-L. Qi, Phys. Rev. B 88 , 195412 (2013).
- [42] M. P. Zaletel, R. S. K. Mong, and F. Pollmann, Phys. Rev. Lett. 110 , 236801 (2013).
- [43] I. Peschel, J. Phys. A 36 , L205 (2003).
- [44] C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81 , 052338 (2010).
- [45] R. B. Laughlin, Phys. Rev. Lett. 50 , 1395 (1983).
- [46] G. Moore and N. Read, Nucl. Phys. B 360 , 362 (1991).
- [47] H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101 , 010504 (2008).
- [48] Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. B 84 , 075128 (2011).
- [49] P. D. Francesco, P. Mathieu, and D. S´n´chal, e e Conformal Field Theory (Springer-Verlag, New York, 1997).
- [50] M. Greiter and R. Thomale, Phys. Rev. Lett. 102 , 207203 (2009).
- [51] M. Greiter, D. F. Schroeter, and R. Thomale, Phys. Rev. B 89 , 165125 (2014).
- [52] P. Lecheminant and A. M. Tsvelik, Phys. Rev. B 95 , 140406 (2017).
- [53] Z.-X. Liu, H.-H. Tu, Y.-H. Wu, R.-Q. He, X.-J. Liu, Y. Zhou, and T.-K. Ng, Phys. Rev. B 97 , 195158 (2018).
- [54] J.-Y. Chen, L. Vanderstraeten, S. Capponi, and D. Poilblanc, Phys. Rev. B 98 , 184409 (2018).
- [55] H.-C. Zhang, Y.-H. Wu, H.-H. Tu, and T. Xiang, Phys. Rev. B 103 , 075130 (2021).
- [56] B. d. z. Jaworowski and A. E. B. Nielsen, Phys. Rev. B 106 , 115131 (2022).
- [57] Y. Huang, W. Zhu, S.-S. Gong, H.-C. Jiang, and D. N. Sheng, Phys. Rev. B 105 , 155104 (2022).
- [58] H. Yao and D.-H. Lee, Phys. Rev. Lett. 107 , 087205 (2011).
- [59] B. S. Shastry, Phys. Rev. B 46 , 8263 (1992).
- [60] W.-W. Luo, Y. Huang, D. N. Sheng, and W. Zhu, Phys. Rev. B 108 , 035130 (2023).
- [61] K. Li, Y.-B. Zhang, and H. C. Po, arXiv:2408.XXXXX.
Tong Liu Ying-Hai Wu Hong-Hao Tu Tao Xiang
Efficient conversion from fermionic Gaussian states to matrix product states
Fermionic Gaussian states are eigenstates of quadratic Hamiltonians and
widely used in quantum many-body problems. We propose a highly efficient
algorithm that converts fermionic Gaussian states to matrix product states. It
can be formulated for finite-size systems without translation invariance, but
becomes particularly appealing when applied to infinite systems with
translation invariance. If the ground states of a topologically ordered system
on infinite cylinders are expressed as matrix product states, then the fixed
points of the transfer matrix can be harnessed to filter out the anyon
eigenbasis, also known as minimally entangled states. This allows for efficient
computation of universal properties such as entanglement spectrum and modular
matrices. The potential of our method is demonstrated by numerical calculations
in two chiral spin liquids that have the same topological orders as the bosonic
Laughlin and Moore-Read states, respectively. The anyon eigenbasis for the
first one has been worked out before and serves as a useful benchmark. The
anyon eigenbasis of the second one is, however, not transparent and its
successful construction provides a nontrivial corroboration of our method.