Finite and Infinite Matrix Product States
for Gutzwiller Projected Mean-Field Wavefunctions

In this section, we show in detail how to convert a Slater determinant into a finite matrix product state. The Schmidt decomposition of a Slater determinant $|SD\rangle$ on $N$ sites bi-partitioned into two regions cut between sites $i$ and $i+1$ will be notated as

$$\ket{SD}=\sum_{\alpha}\lambda^{i;N}_{\alpha}\ket{L_{\alpha}^{i;N}}\ket{R_{\alpha}^{i;N}}$$

(1)

where $|L_{\alpha}\rangle$ and $|R_{\alpha}\rangle$ are the $\alpha$’th left and right Schmidt vectors respectively (with support in their respective subsystem) and $\lambda_{\alpha}$ is the $\alpha$’th Schmidt eigenvalue. Note that, for a Slater determinant, each of the individual left and right Schmidt vectors are also Slater determinants and efficiently computable cheong2004many ; klich2006lower ; knizia2012density ; peschel2012special . Slater determinants are specified by a set of single-particle orbitals and all the Slater determinants in the set of right Schmidt vectors $\{|R_{\alpha}\rangle\}$ are specified by subsets of single particle orbitals from a set of (at most) $N$ single-particle orbitals $\{\phi^{R}_{1}...\phi^{R}_{N}\}$ defined on the (inclusive) sites $[(i+1),\ldots,N]$. There are, at most, $2^{N}$ such subsets. Analogous statements hold for the left Schmidt vectors.

We now describe how to efficiently evaluate the matrix elements of each $A$. We start by noting that $|\sigma_{i+1}|\rangle\otimes|R^{i+1;N}\rangle$ is also a Slater determinant. It is specified by the single particle orbitals

$$\{[0\phi_{a}],[0\phi_{b}],\cdots,[0\phi_{c}],\phi^{i+1}\}$$

(4)

where $[0\phi_{a}]$ is the single particle orbital with coefficients in the lattice basis $[0,\phi_{a}(i+2),\phi_{a}(i+3),\cdots,\phi_{a}(N)]$ and $\phi^{i+1}$ is the single-particle orbital in analogous notation, $[1,0,\cdots,0]$. Eqn. 3 then reduces to the overlap of two Slater determinants of size $(N-i)\times(N-i)$ which can be computed in $O(N^{3})$ time.

While naively each element of $A$ requires such a computation, there is a significant overlap in these different computations which reduce the naive compuational complexity of the tensor computation. There are two steps in computing the overlap of two Slater determinants: evaluating the overlap matrix between all pairs of single particle orbitals that make up the two determinants and computing the determinant of this overlap matrix. All the Slater determinants used in the ket (respectively bra) of eqn. 3 (over different terms in A) come from a subset of single-particle orbitals of the N-orbital set $\{\phi_{1}^{R},..,\phi_{N}^{R}\}$. We can compute the overlap matrix of all these respective single-particle orbitals once per three-tensor $A$ at a cost of $O(N^{3}).$ The entries of $A$ are then determinants of submatrices of this overlap matrix. While naively each determinant also costs $O(N^{3})$ to compute, the submatrices differ only in the bottom $\log_{2}D$ columns and right $\log_{2}D$ rows where $D$ is the bond-dimension of $A$; determinant update formulas can then be used to accelerate this computation letting each determinant be computed in time $O(N^{2}\log_{2}D)$ after an initial $O(N^{3})$ operation to evaluate the inverse of the upper-left $(N-\log_{2}D)\times(N-\log_{2}D)$ block of the overlap matrix. The whole evaluation of each tensor $A$ can be done in $O(N^{3})+O(D^{2}N^{2}\log_{2}D)$ time. This can be further attenuated somewhat by more aggressive use of determinant update formulas harville1998matrix .

Notice that there are significant parts of this algorithm that can be run in parallel. Each three-tensor $A$ can be computed separately. Within each $A,$ the Schmidt decomposition can be partially parallelized; each element of the overlap matrix can be computed in parallel; and, after the initial evaluation of the inverse of the upper-left block of the overlap matrix, each determinant can then be computed in parallel.

The above procedure generates a finite MPS approximation (the accuracy of the representation is given as an input to the algorithm) of any Slater determinant. In this section, we describe how to generate an infinite MPS from the Slater determinant ground state of a gapped mean-field Hamiltonian. This infinite MPS can be described by left $L$ and right $R$ boundary tensors which sandwich the bulk tensor $A$ giving us an infinite matrix product state of the form,

$$\ket{\textrm{iMPS}}=\sum_{\sigma}L^{\sigma_{L}}\ldots A^{\sigma_{n-1}}A^{\sigma_{n}}A^{\sigma_{n+1}}\ldots R^{\sigma_{R}}\times\\ \ket{\sigma_{L}\ldots\sigma_{n-1}\sigma_{n}\sigma_{n+1}\ldots\sigma_{R}}$$

(5)

with an arbitrary number of bulk tensors $A$. $L$ and $R$ are tensors which span a fixed number $k$ of sites. Note that any thermodynamic observable can be computed directly in the thermodynamic limit of the Slater Determinant using only the bulk tensor $A$. In addition, we can compute the amplitude for the Slater determinant on any (large enough) system size, by inserting the corresponding number of bulk tensors between the boundary tensors $L$ and $R$ (i.e. to generate the MPS for a $N$ site Slater determinant from the infinite MPS, we therefore use $N-2k$ bulk tensors $A$); see fig.3.

To generate the iMPS, we start off by producing two Slater determinants defined on $2N$ and $2N+1$ sites, where $N$ is sufficiently large such that the entanglement spectrum is constant over cuts in the “bulk” of the wave-functions. For gapped systems, we generically expect the entanglement spectrum over the bulk to be constant; see fig.  S3 in supplement  S5 for an example of this. We then generate the Schmidt decompositions

	$\displaystyle|\Psi^{2N}\rangle$	$\displaystyle=$	$\displaystyle\sum_{\alpha}\lambda^{N;2N}_{\alpha}|L_{\alpha}^{N;2N}\rangle|R_{\alpha}^{N;2N}\rangle$		(6)
	$\displaystyle|\Psi^{2N+1}\rangle$	$\displaystyle=$	$\displaystyle\sum_{\alpha}\lambda^{N;2N+1}_{\alpha}|L_{\alpha}^{N;2N+1}\rangle|R_{\alpha}^{N;2N+1}\rangle$		(7)

Both $|L_{\alpha}^{N;2N}\rangle$ and $|L_{\alpha}^{N;2N+1}\rangle$ are going to be the same up to a gauge freedom. We fix this gauge freedom by defining a unitary

	$\displaystyle C_{\alpha\beta}^{N;2N+1}=$	$\displaystyle 0$	$\displaystyle\textrm{ if }\lambda_{\alpha}\neq\lambda_{\beta}$		(8)
	$\displaystyle=$	$\displaystyle\bra{L^{N;2N}}_{\alpha}\ket{L^{N;2N+1}}_{\beta}$	otherwise

which rotates between Schmidt eigenvectors with the same Schmidt eigenvalue allowing the state on $2N+1$ sites to be defined as

$$\ket{\Psi^{2N+1}}=\sum_{\alpha\gamma}\ket{L^{N;2N}}_{\alpha}\lambda^{N;2N+1}_{\alpha}C^{N;2N+1}_{\alpha\gamma}\ket{R^{N;2N+1}}_{\gamma}$$

(9)

Then the tensor $A$ for the iMPS is

$$A^{\sigma_{N+1}}_{\alpha\beta}=\sum_{\gamma}C^{[2N+1]}_{\alpha\gamma}\bra{\sigma_{N+1}|\bra{R^{N;2N}}_{\beta}}\ket{R^{N,2N+1}}_{\gamma}$$

(10)

As in the finite MPS case, we have that the single-particle orbitals of the Slater determinant $|\sigma_{N+1}\rangle|R^{N;2N}\rangle$ are shifted to the right with an initial zero as their first element. The overlap of this tensor can be computed in exactly the same way as for the finite MPS case. Here, though, we only need to evaluate one tensor $A$ instead of a tensor per site, with the assumption that we are using an iMPS defined by a single tensor (i.e. single site unit cell) $A$. This process can be generalized to multi-site unit cells as well (see supplement S2). Note that by directly applying the finite MPS algorithm to large systems to try to find the bulk tensor $A$ will fail because the gauge freedom available in the tensors will prevent a single identical bulk tensor from being produced at each step.

We numerically validate our algorithms by applying fMPStoSD and iMPStoSD on the ground state of the Su-Schrieffer-Heeger(SSH) model su1979solitons ,

$$H_{SSH}=v\sum_{n}\left(c^{\dagger}_{n,1}c_{n,2}+\text{h.c.}\right)+\\ w\sum_{n}\left(c^{\dagger}_{n+1,1}c_{n,2}+\text{h.c.}\right)$$

(11)

The model describes spinless fermions on a 1D lattice, with a 2 site unit cell made up of $A$,$B$ sites, with different (real) parameters for intra-cell hopping ($v$) and intercell hopping ($w$). It admits two different quantum ground states, distinct in their topological properties: a trivially gapped phase for $v>w$ and a (symmetry protected) topological gapped ground state, characterized by the presence of two zero-energy edge modes inside the gap, for $v<w$, separated by a quantum critical point at $v=w$.

We will discuss here the trivial ground state. A small subtlety related to choosing the same gapless boundary mode in the Slater determinant wave-functions used for generating the uniform tensor in the iMPS procedure is delegated to the supplement S4. For the finite Slater Determinant, we compare the MPS we generate using two different truncation values against the exact Slater Determinant by comparing all of the amplitudes (see 1st column in fig. 5). For the infinite case, we generate the iMPS and then use the bulk tensor we have computed along with the boundary tensors to compute amplitudes for a much larger system and again compare amplitudes against the exact solution for that much larger system (see 2nd column in fig. 5). In both cases, we find that the amplitudes are in very good agreement for all amplitudes down to the Schmidt eigenvalue cutoff.

The key trick to convert a BdG wave-function into a MPS will be to 1) convert it to a Slater determinant through a particle hole transformation, 2) convert this Slater determinant to a MPS, and 3) then undo the particle-hole transformation in the MPS language.

Consider a BdG Hamiltonian:

$$H_{BdG}=-\sum_{\langle ij\rangle,\sigma}t_{ij}\left(c^{\dagger}_{i\sigma}c_{j\sigma}+\text{h.c.}\right)-\\ \sum_{\langle ij\rangle}\Delta_{ij}\left(c^{\dagger}_{i\uparrow}c^{\dagger}_{j\downarrow}+\text{h.c.}\right)-\mu\sum_{i\sigma}c^{\dagger}_{i\sigma}c_{i\sigma}$$

(12)

Under a canonical particle-hole transformation in the $\downarrow$ sector:

$$\begin{split}f^{\dagger}_{2i-1}=c^{\dagger}_{i\uparrow}\\ f^{\dagger}_{2i}=c_{i\downarrow}\end{split}$$

(13)

the BdG Hamiltonian becomes:

$$H^{ph}_{\text{BdG}}=-\sum_{\langle ij\rangle}t_{ij}\left(f^{\dagger}_{2i-1}f_{2j-1}-f^{\dagger}_{2i}f_{2j}+\text{h.c.}\right)-\\ \sum_{\langle ij\rangle}\Delta_{ij}\left(f^{\dagger}_{2i-1}f_{2j}+\text{h.c.}\right)-\mu\sum_{i\sigma}\left(f^{\dagger}_{2i-1}f_{2i-1}-f^{\dagger}_{2i}f_{2i}\right)$$

(14)

and the new vacuum is $\ket{0^{\text{ph}}}=c^{\dagger}_{1\downarrow}c^{\dagger}_{2\downarrow}\ldots c^{\dagger}_{N\downarrow}\ket{0}$, where $\ket{0}$ is the vacuum of the original theory. The ground state of $H_{\text{BdG}}$ is thus the Slater determinant ground state of $H^{\text{ph}}_{\text{BdG}}$ on top of the new vacuum $\ket{0^{\text{ph}}}$.

Using the results from section II, we convert the Slater determinant ground state of $H^{\text{ph}}_{\text{BdG}}$ into a MPS $\ket{\Psi_{MPS}}=\sum_{\{\sigma\}}A^{[1]\sigma_{1}}\ldots A^{[2N]\sigma_{2N}}\ket{\sigma_{1}\ldots\sigma_{2N}}$ where $\sigma_{2i-1}=\{0,1\}$ ($i\in[1,N])$ indicates the absence/presence of a $\uparrow$ particle and $\sigma_{2i}=\{0,1\}$ indicates the absence/presence of a hole on top of the filled $\downarrow$ Fermi sea at site $i$.

To ‘undo’ the particle-hole transformation, we need to deal with the fact that the $f_{i}$ act on the false vacuum $\ket{0^{ph}}$ (and not the real vacuum) by swapping, for all $i$, the matrices $A^{[2i]1}$ and $A^{[2i]0}$. Moreover, by ordering the fermionic operators by site, and then spin, the matrices $A^{[2i-1]1}$, $A^{[2i]1}$ will pick up factors of $(-1)^{i-1}$. We can now combine these transformations giving us our final MPS for the BdG ground state of the form $|\Phi_{GS}\rangle=\sum_{\{\sigma={0,\uparrow,\downarrow,\uparrow\downarrow\}}}B^{[1]\sigma_{1}}\ldots B^{[N]\sigma_{N}}|\sigma_{1}\ldots\sigma_{N}\rangle$ where

	$\displaystyle B^{[i]0}$	$\displaystyle=(-1)^{i-1}\times A^{[2i-1]0}A^{[2i]1}$
	$\displaystyle B^{[i]\uparrow}$	$\displaystyle=(-1)^{2i-2}\times A^{[2i-1]1}A^{[2i]1}$
	$\displaystyle B^{[i]\downarrow}$	$\displaystyle=(-1)^{0}\times A^{[2i-1]0}A^{[2i]0}$
	$\displaystyle B^{[i]\uparrow\downarrow}$	$\displaystyle=(-1)^{i-1}\times A^{[2i-1]1}A^{[2i]0}$		(15)

This approach works both for the finite and infinite MPS as we just used our (i)MPS $\rightarrow$ Slater determinant approach as a subroutine. For the infinite MPS it produces a unit cell of size 2 as every other $B$ differs by a sign.

As a check of our algorithm, we consider the ground state of the BdG Hamiltonian of the form:

$$H_{BdG}=-\sum_{\langle ij\rangle,\sigma}t_{ij}\left(c^{\dagger}_{i\sigma}c_{j\sigma}+\text{h.c.}\right)-\\ \sum_{\langle ij\rangle}\Delta_{ij}\left(c^{\dagger}_{i\uparrow}c^{\dagger}_{j\downarrow}+c^{\dagger}_{j\uparrow}c^{\dagger}_{i\downarrow}+\text{h.c.}\right)-\mu\sum_{i\sigma}c^{\dagger}_{i\sigma}c_{i\sigma}$$

(16)

and compare the amplitudes of the exact ground state with the MPS generated (see 3rd column in fig. 5).

We will show how to generate a MPS representa-tion of the Pfaffian ground state of the Kitaev p+ip chain:

$$H_{c}=-t\sum_{n}\left(c^{\dagger}_{n}c_{n+1}+\text{h.c.}\right)+\Delta\sum_{n}\left(c^{\dagger}_{n}c^{\dagger}_{n+1}+\text{h.c.}\right)\\ +\mu\sum_{n}c^{\dagger}_{n}c_{n}$$

(17)

We first consider $H^{\text{ext}}=H_{c}\bigoplus H_{d}$ whose ground state is given by the tensor product of two identical pfaffians:

$$\ket{GS}=\sum_{\sigma_{c},\sigma_{d}}Pf(M_{\sigma_{c}})\ket{\sigma_{c}}\bigotimes Pf(M_{\sigma_{d}})\ket{\sigma_{d}}$$

(18)

where $M$ is a $N\times N$ matrix built from parameters of the model and $M_{\sigma_{c}}$ is a submatrix of $M$ obtained by selecting indices as given by the $\sigma_{c}$ configuration.

Using the local canonical transformation of fermions,

$$\begin{split}c^{\dagger}_{n}=\frac{\left(\bar{c}^{\dagger}_{n,\uparrow}+\bar{c}^{\dagger}_{n,\downarrow}\right)}{\sqrt{2}}\\ d^{\dagger}_{n}=i\frac{\left(\bar{c}^{\dagger}_{n,\uparrow}-\bar{c}^{\dagger}_{n,\downarrow}\right)}{\sqrt{2}}\end{split}$$

(19)

converts $H^{\text{ext}}$ to

$$H_{BdG}=-t\sum_{n,\sigma}\left(\bar{c}^{\dagger}_{n,\sigma}\bar{c}_{n+1,\sigma}+\text{h.c.}\right)\\ +\Delta\sum_{n}\left(\bar{c}^{\dagger}_{n\uparrow}\bar{c}^{\dagger}_{n+1\downarrow}+\bar{c}^{\dagger}_{n+1,\downarrow}\bar{c}^{\dagger}_{n,\uparrow}+\text{h.c.}\right)\\ +\mu\sum_{n\sigma}\bar{c}^{\dagger}_{n,\sigma}\bar{c}_{n,\sigma}$$

(20)

The transformation leaves the vacuum unchanged. Given a BdG Hamiltonian, we can obtain the ground-state as an MPS as done in section III.1.

$$\ket{GS}=\sum_{\sigma}\ldots A^{2i-1\sigma_{2i-1}}A^{2i\sigma_{2i}}\ldots\ket{\sigma_{1}\sigma_{2}\ldots\sigma_{2N}}$$

(21)

where $A^{2i-1\sigma_{2i-1}}$ is the tensor on site $i$ for the $\uparrow$ local physical sector and $A^{2i\sigma_{2i}}$ is the tensor on site $i$ for the $\downarrow$ local physical sector.

Notice that the canonical transformation given in eqn. 19 mixes the $\uparrow$, $\downarrow$ physical sectors on site $i$. Hence we can obtain the MPS for the GS in the $c,d$ space by choosing the on-site tensor in the following way:

$$\ket{GS}=\sum_{\sigma}C^{1\sigma_{1}}C^{2\sigma_{2}}\ldots C^{N\sigma_{N}}\ket{\sigma_{1}\sigma_{2}\ldots\sigma_{N}}$$

(22)

with $\sigma=\{0,c,d,cd\}$ where:

$$\begin{split}C^{n,0}=A^{2n-1,0}A^{2n,0}\\ C^{n,c}=\frac{A^{2n-1,1}A^{2n,0}+A^{2n-1,1}A^{2n,0}}{\sqrt{2}}\\ C^{n,d}=i\frac{A^{2n-1,1}A^{2n,0}-A^{2n-1,1}A^{2n,0}}{\sqrt{2}}\\ C^{n,cd}=iA^{2n-1,1}A^{2n,1}\end{split}$$

(23)

By projecting out the $d$ particles in the $\ket{GS}$ wavefunction, we obtain a Pfaffian wavefunction in the $c$-particle sector: $\ket{GS}=\text{constant}\times\sum_{\sigma_{c}}Pf(M_{\sigma_{c}})\ket{\sigma_{c}}$. At the level of the MPS this projection is realized by eliminating the sectors $\sigma=\{d,cd\}$.

$$\ket{\text{Pf}}=\sum_{\sigma=\{0,c\}}C^{1\sigma_{1}}C^{2\sigma_{2}}\ldots C^{N\sigma_{N}}\ket{\sigma_{1}\sigma_{2}\ldots\sigma_{N}}$$

(24)

While we focused in the previous section on a specific example, here we consider a generic quadratic Hamiltonian $H=\sum_{n,m}C^{\dagger}_{n}h_{n,m}C_{m}$ with $g$ species of fermions per unit cell where the vector $C^{\dagger}_{n}=\left(c^{\dagger}_{n,1},c_{n,1},c^{\dagger}_{n,2},c_{n,2}\ldots c^{\dagger}_{n,g},c_{n,g}\right)$.

We form an extended Hamiltonian $H^{\text{ext}}$ which is a sum of two copies of $H$:

$$H^{ext}=\sum h^{\text{hop}}_{i\alpha,j\beta}\left(c^{\dagger}_{i,\alpha}c_{j,\beta}+d^{\dagger}_{i,\alpha}d_{j,\beta}+\text{h.c.}\right)\\ +\sum h^{\text{pair}}_{i\alpha,j\beta}\left(c^{\dagger}_{i,\alpha}c^{\dagger}_{j,\beta}+d^{\dagger}_{i,\alpha}d^{\dagger}_{j,\beta}+\text{h.c.}\right)$$

(25)

where $\alpha,\beta\in(1,\ldots,g)$. As before, its ground state $\ket{\text{GS}}^{\text{ext}}=\sum_{\sigma_{c},\sigma_{d}}Pf(M_{\sigma_{c}})\ket{\sigma_{c}}\bigotimes Pf(M_{\sigma_{d}})\ket{\sigma_{d}}$ is a tensor product of two identical Pfaffian wavefunctions. We then obtain the Pfaffian ground state of $H$ by projecting out all the $d$ sectors. Under the following linear canonical transformation

$$\begin{split}\bar{c}^{\dagger}_{n,\alpha,\uparrow}=\frac{c^{\dagger}_{n,\alpha}+id^{\dagger}_{n,\alpha}}{\sqrt{2}}\\ \bar{c}^{\dagger}_{n,\alpha,\downarrow}=\frac{c^{\dagger}_{n,\alpha}-id^{\dagger}_{n,\alpha}}{\sqrt{2}}\end{split}$$

(26)

$H^{\text{ext}}$ becomes a BdG-like hamiltonian when expressed in terms of $c_{\uparrow}$ and $c_{\downarrow}$:

$$\begin{split}H^{ext}_{\text{BdG}}=\sum h^{\text{hop}}_{i\alpha,j\beta}\left(\bar{c}^{\dagger}_{i,\alpha,\uparrow}\bar{c}_{j,\beta,\uparrow}+\bar{c}^{\dagger}_{i,\alpha,\downarrow}\bar{c}_{j,\beta,\downarrow}+\text{h.c.}\right)+\\ +\sum h^{\text{pair}}_{i\alpha,j\beta}\left(\bar{c}^{\dagger}_{i,\alpha,\uparrow}\bar{c}^{\dagger}_{j,\beta,\downarrow}+\bar{c}^{\dagger}_{i,\alpha,\downarrow}\bar{c}^{\dagger}_{j,\beta,\uparrow}+\text{h.c.}\right)\end{split}$$

(27)

We can then solve for the MPS representation of the groundstate of the above BdG-like Hamiltonian using the methods described in section  III.1.

We obtain the Slater determinant ground state $\bra{\Psi_{ext}}$ by diagonalizing the particle-hole transformed $H^{\text{ext}}_{\text{BdG}}$ and then computing its MPS representation. Each unit cell $M$ is described by $2g$ tensors $A^{[M,p]\sigma}$ with $\sigma=0,1$ signifying the absence/presence of a particle of type $p\in[0,1,\ldots 2g-1]$. A particle of type $p=2k$ corresponds to the flavor $k,\uparrow$; a particle of type $p=2k+1$ corresponds to the flavor $k,\downarrow$. From the above MPS (which is in $\bar{c}_{\uparrow}$, $\bar{c}_{\downarrow}$ local physical space) we construct the MPS tensors in $c$, $d$ space. In particular, the matrices describing the absence/presence of a particle of type $c_{i}$ on site $M$ are given by:

	$\displaystyle B^{[M,i]0}$	$\displaystyle=A^{[M,2i-1][0]}A^{[M,2i][1]}$		(28)
	$\displaystyle B^{[M,i]1}$	$\displaystyle=\frac{(-1)^{g(M-1)+i-1}}{\sqrt{2}}\times$
		$\displaystyle\left[A^{[M,2i-1][1]}A^{[M,2i][1]}+A^{[M,2i-1]0}A^{[M,2i][0]}\right]$

where $(-1)^{g(M-1)+i-1}$ takes care of fermionic ordering. and the MPS representation of $\ket{\Psi_{GS}}$ defined on $Ng$ sites is given by:

$$\ket{\Psi_{GS}}=\sum_{\{\sigma\}}\left(B^{[1,1]\sigma_{1_{1}}}B^{[1,2]\sigma_{1_{2}}}\ldots B^{[1,g]\sigma_{1_{g}}}\right)\ldots\times\\ \left(B^{[N,1]\sigma_{N_{1}}}B^{[N,2]\sigma_{N_{2}}}\ldots B^{[N,g]\sigma_{N_{g}}}\right)\times\\ \ket{\left(\sigma_{1_{1}}\sigma_{1_{2}}\ldots\sigma_{1_{g}}\right)\ldots\left(\sigma_{N_{1}}\sigma_{N_{2}}\ldots\sigma_{N_{g}}\right)}$$

(29)

By suitably contracting tensors we can obtain a $N$-tensor MPS representation with physical dimension $2^{g}$: $\ket{Pf}=\sum_{\sigma}C^{1\sigma_{1}}C^{2\sigma_{2}}\ldots C^{N\sigma_{N}}\ket{\sigma_{1}\sigma_{2}\ldots\sigma_{N}}$. For instance, the tensor corresponding to the presence of particles of type $t_{1},t_{2},\ldots t_{s}$ on site $M$ is

$$C^{[M](t_{1},t_{2},\ldots,t_{s})}=B^{[M,1][0]}\ldots\times\\ B^{[M,t_{1}][1]}B^{[M,t_{1}+1][0]}\ldots B^{[M,t_{s}][1]}\ldots B^{[M,g][0]}$$

(30)

In this section we describe how to obtain the MPS representation of a wavefunction

$$\ket{\psi_{1/n}}=\sum_{r_{1},r_{2},\cdots r_{n}}\braket{r_{1},r_{2},\ldots r_{n}}{\psi}^{n}|r_{1},r_{2},\ldots,r_{n}\rangle$$

(31)

where $\ket{\psi}$ is a Slater determinant. We will use as an example n = 3. Products of other mean-field wave-functions can be obtained similarly.

We extend our N-site system to a $3N$-site system for which we label the sites as: $\{1_{1}1_{2}1_{3}2_{1}2_{2}2_{3}\ldots N_{1}N_{2}N_{3}\}$. We then write a single Slater determinant (by padding and interlacing the orbitals to keep the above ordering) of the form $\ket{\psi}\bigotimes\ket{\psi}\bigotimes\ket{\psi}$ for which we then convert into an MPS given by

$$\ket{MPS}=\sum_{{i_{p}}}\left(B^{[1_{1}]i_{1_{1}}}B^{[1_{2}]i_{1_{2}}}B^{[1_{3}]i_{1_{3}}}\right)\ldots\times\\ \left(B^{[N_{1}]i_{N_{1}}}B^{[N_{2}]i_{N_{2}}}B^{[N_{3}]i_{N_{3}}}\right)\times\\ \ket{i_{1_{1}}i_{1_{2}}i_{1_{3}}\ldots i_{N_{1}}i_{N_{2}}i_{N_{3}}}$$

(32)

Projecting on the sector $i_{n_{1}}=i_{n_{2}}=i_{n_{3}}$ gives us the desired results of

$$\ket{\psi_{1/3}}=A^{[1]i_{1}}A^{[2]i_{2}}\ldots A^{[N]i_{N}}\ket{i_{1}i_{2}\ldots i_{N}}$$

(33)

where we define

$$A^{[n]i_{n}}=B^{[n_{1}]i_{n_{1}}}B^{[n_{2}]i_{n_{2}}}B^{[n_{3}]i_{n_{3}}}$$

(34)

The relevant mean-field Hamiltonian that arises from the parton construction of the bilinear-biquadratic S=1 model liu2010fermionic ,liu2014gutzwiller is

$$H_{mf}=-J\chi\sum_{i,\alpha=-1,0,1}\left[c^{\dagger}_{i,\alpha}c_{i,\alpha}+\text{h.c.}\right]\\ +(J-K)\Delta\sum_{i,j}\left[c^{\dagger}_{i,-1}c^{\dagger}_{j,1}-c^{\dagger}_{0i}c^{\dagger}_{0j}+c^{\dagger}_{1i}c^{\dagger}_{-1j}+\text{h.c.}\right]\\ +\lambda\sum_{i,\alpha}c^{\dagger}_{i\alpha}c_{i\alpha}$$

(37)

where the $c^{\dagger}_{-1}$, $c^{\dagger}_{0}$ and $c^{\dagger}_{1}$ are the on-site fermion parton flavors corresponding to $S_{z}=-1,0,1$.

This could be converted into a MPS by treating it as a general Pfaffian and then applying the techniques in section  III.2. In models such as this, though, where the mean-field Hamiltonian in the parton basis has a tensor sum structure where one or more of the Hilbert subspaces can be treated with a simpler mean-field (i.e. with a SD or BdG ground state) it makes computational sense to obtain the MPS representation in each sector and then “glue” the two MPS together; we exemplify this approach here.

Introducing a Nambu spinor, in the $k$-basis the Hamiltonian is block-diagonal:

	$$\displaystyle H^{k}_{mf}=\frac{1}{2}\begin{bmatrix}c^{\dagger}_{k,1}&c_{-k,-1}&c^{\dagger}_{k,0}&c_{-k,0}\end{bmatrix}\times$$
	$$\displaystyle\begin{bmatrix}\chi_{k}&\Delta_{k}&0&0\\ \Delta_{k}^{*}&-\chi_{k}&0&0\\ 0&0&\chi_{k}&-\Delta_{k}\\ 0&0&-\Delta_{k}^{*}&-\chi_{k}\end{bmatrix}\begin{bmatrix}c_{k,1}\\ c^{\dagger}_{-k,-1}\\ c_{k,0}\\ c^{\dagger}_{-k,0}\end{bmatrix}$$		(38)

The one-body hamiltonian is a tensor sum of BdG-like Hamiltonian $H_{BdG}$ and a $p+ip$ Hamiltonian $H_{p+ip}$. The mean field ground state is (we consider anti-periodic boundary conditions and an even number of sites):

$$\ket{\Psi_{GS}}=\Pi_{0<k<2\pi}\left(u_{k}+v_{k}c^{\dagger}_{k,1}c^{\dagger}_{-k,-1}\right)\times\\ \Pi_{0<q<\pi}\left(u_{q}-v_{q}c^{\dagger}_{q,0}c^{\dagger}_{-q,0}\right)\ket{0}$$

(39)

where $u_{k}$ and $v_{k}$ are given in terms of the parameters of the Hamiltonian.

By performing a particle-hole transformation in the $S_{z}=\{\uparrow,\downarrow\}$ sector (see section III.1) and the Pfaffian artificial extension in the $S_{z}=0$ sector (see section III.2),

	$\displaystyle f^{\dagger}_{1,k}$	$\displaystyle=c^{\dagger}_{k,1}$		(40)
	$\displaystyle f^{\dagger}_{2,k}$	$\displaystyle=c_{-k,-1}$
	$\displaystyle f^{\dagger}_{3,k}$	$\displaystyle=c^{\dagger}_{k,0}$
	$\displaystyle f^{\dagger}_{4,k}$	$\displaystyle=c_{-k,0}$

where $\{f_{k,\alpha},f_{q,\beta}\}=\delta_{\alpha\beta}\delta_{kq}$, giving us

$$\ket{\Psi_{GS}^{ext}}=\Pi_{k}\left(u_{k}f^{\dagger}_{k,2}+v_{k}f^{\dagger}_{k,1}\right)\\ \Pi_{0<q<\Pi}\left(u_{q}f^{\dagger}_{q,4}-v_{q}f^{\dagger}_{k,3}\right)\ket{\text{vac}}$$

(41)

with $\ket{\text{vac}}=\Pi_{0<k<2\pi}c_{k\downarrow}\ket{0}\Pi_{0<q<\pi}c_{k,0}\ket{0}$. Since $\ket{\Psi^{ext}_{GS}}$ is a Slater Determinant, we can obtain the MPS representation using the methods in section II. We now “undo” the transformation (see again section III.1 and section III.2) and write the MPS in the following form:

$$\ket{MPS}=\sum_{\{i\}}\left(C^{[1_{\uparrow}]i_{1\uparrow}}C^{[1_{\downarrow}]i_{1\downarrow}}C^{[1,\rightarrow]i_{1,\rightarrow}}\right)\times\\ \ldots\left(C^{[N_{\uparrow}]i_{N\uparrow}}C^{[N_{\downarrow}]i_{N\downarrow}}C^{[N,\rightarrow]i_{N_{\rightarrow}}}\right)\\ \ket{\left(i_{1_{\uparrow}}i_{1_{\downarrow}}i_{1_{\rightarrow}}\right)\ldots\left(i_{N_{\uparrow}}i_{N_{\downarrow}}i_{N_{\rightarrow}}\right)}$$

(42)

where $i_{n,\alpha}\in\{0,1\}$ indicates the absence/presence of a particle of type $\alpha$ on site $n$.

To obtain the MPS with onsite tensors $A^{[n]\sigma_{n}},\sigma_{n}\in\{\uparrow,\downarrow,\rightarrow,\uparrow\rightarrow,\downarrow\rightarrow,\uparrow\downarrow,\uparrow\downarrow\rightarrow\}$, we ”glue” together appropriate sectors. For example

$$A^{[n]\uparrow}=C^{[n\uparrow]1}C^{[n\downarrow]0}C^{[n\rightarrow]0}$$

(43)

Gutzwiller projection is realized by summing only over the one-particle per site physical indices $\sigma_{n}\in\{\uparrow,\downarrow,\rightarrow\}$:

$$P_{G}\ket{\Psi_{GS}}=\sum_{\sigma_{n}\in\{\uparrow,\downarrow,\rightarrow\}}A^{[1]\sigma_{1}}\ldots A^{[N]\sigma_{N}}\ket{\sigma_{1}\ldots\sigma_{N}}$$

(44)

In this section, we will discuss orthogonalizing our iMPS states. This includes a brief overview of the standard iMPS orthogonalization as well as a detailed description of how we address the degeneracies that appear when Gutzwiller projecting slave-fermion mean-field states onto degenerate ground state manifolds.

The orthogonalization procedure for a typical iMPS is standard (see ref. orus2008infinite, and supplement  S6 for more details). The method relies on obtaining the leading right/left eigenvectors of the transfer matrix operator $E=\sum_{\sigma}A^{\sigma}\otimes A^{\sigma*}$ where $\sigma$ runs over the on-site physical index. $E$ admits the following decomposition:

$$E=\sum_{i}\lambda_{i}\ket{R}_{i}\bra{L}_{i}$$

(45)

where $\ket{L}_{i}$ and $\ket{R}_{i}$ are left/right eigenvectors of $E$ and $\braket{L_{i}}{R_{i}}=0$ for $\lambda_{i}$ non-degenerate. In the infinite limit only the leading left/right eigenvectors of $E$ survive. If the dominant eigenvalue is non-degenerate, the transfer matrix is given by

$$\lim_{N\rightarrow\infty}E^{N}=\ket{R}\bra{L}$$

(46)

where $\ket{R(L)}$ are by definition the eigenvectors corresponding to the dominant eigenvalue. Thus, the dominant left and right eigenvectors correspond to a pure state. The Implicitly Restarted Arnoldi Method can be efficiently used for this purpose by noting that $\left(\sum_{\sigma}A^{\sigma}\otimes A^{\sigma*}\right)\mathrm{vec}(v)=\sum_{\sigma}\mathrm{vec}\left(A^{\sigma*}vA^{T}\right)$, where the $\mathrm{vec}(v)$ operation takes the square matrix $v$ and stacks the columns together. The entanglement spectrum and observables are then easily obtained.