For $i=z,x$ we thus obtain $|ψ_f⟩=±|ψ_f⟩$, showing that states with $_2$ charge $-1$ are mapped to zero. For $i=1,y$ on the other hand, $σ^x σ^ϕ(i) σ^x = - σ^ϕ(i)$ resulting in $|ψ_f⟩=∓|ψ_f⟩$ and thus a zero projection for $|ψ_s⟩$ with $_2$ charge $+1$.\footnote{Note that in this calculation bosonizing the virtual legs was quite helpful as it allowed to work with multiple, otherwise odd, $\sigma^x$ tensors. Alternatively, one could retain the fermionic degrees of freedom and introduce some consistent form of ordering. This is indispensable for a 2D extension and explained in more detail in Ref.~\cite{Obrien}.} If we would only have used the trivial boundary condition, $i=1$, the MPO would only map from the $_2$ sector with charge $+1$. Analogous considerations on the fermion side show that this $U^1$ would also only target parity even fermion states. We conclude that we need the four $U^i$ to arrive at a complete map between the full spin and fermion Hilbert spaces. In Table \ref{tab:JWtable} we summarize how the different charge sectors are related by these MPO's.

Not only charge sectors are redistributed, also boundary conditions are mapped in a non-trivial way by $U^i$. To see this, consider the translation operator $T_f$ on the fermionic side, defined as
\begin{equation}
T_f a_n T_f^{-1} = a_{n+1}
\end{equation}
for $j<N$ and $T_f a_N T_f^-1 = ±a_1$ for (anti-)periodic boundary conditions so that $T_f^N= ±1$. Diagrammatically $T_f$ can be represented as
\begin{equation}
T_f = \diagram{jw}{20}
\end{equation}
with the so-called twist $σ^z_f=P$ only present in the case of anti-periodic boundary conditions. Analogously, a translation operator $T_s$ can be defined on the spin side with a $σ^x$ twist for anti-periodic boundary conditions. Now suppose we start from a translation-invariant spin state $|ψ_s⟩$ with periodic boundary conditions, i.e., $T_s |ψ_s⟩ = e^i k |ψ_s⟩$ with $k=2π/Nn$ and $n ∈{1,..,N}$. Translating the corresponding fermion state $|ψ_f⟩$ then yields
\begin{equation}
\begin{split}
T_f \ket{\psi_f} &= T_f U^i T_s^{-1} T_s \ket{\psi_s} \\[5pt]
&= e^{ik} \, \diagram{jw}{21}\, .
\end{split}
\end{equation}
Here, the composition of $U^i$ with the translation operators can be related to the original MPO by moving the boundary condition to the new final site. Based on the form of this boundary condition, this procedure can generate twists on the spin and/or fermion side. If $i=1$ for instance, the bosonized boundary condition is $σ^z$ and moving it through the bosonized GHZ and CNOT tensors, we obtain a $σ^z_f$ twist on the fermion side. However, we also reorder the fermion legs and as mentioned before, this generates an additional fermionic parity factor, $(-1)^|f_N|∑_i=1^N-1|f_i|$. Because the total fermion charge is even for $i=1$, this factor reduces to $(-1)^|f_N|$, i.e., yet another $σ^z_f$ on the new first fermion leg. Hence, the twist is removed. Similarly, there is no twist by moving the bosonized BC for $i=z$ but the reordering of the fermion legs introduces a $σ^z_f$. This shows that $|ψ_f⟩$ is not an eigenvector of the periodic $T_f$ then but rather of the anti-periodic version as the latter compensates the twist. We conclude that $U^1$ will map the periodic sector on the spin side to the periodic sector on the fermion side while $U^z$ targets the anti-periodic sector on the fermion side. Similar considerations amount to the twists summarized in Table \ref{tab:JWtable}. Note again that we need all four MPO's to fully map the complete spaces.

%$\tilde{G}^1$ effectively is a GHZ-tensor in the $x$-basis, while $\tilde{G}^0$ and $\tilde{H}$ are those in the $z$-basis. Their symmetry properties can be used to check all known JW mappings.

%In this discussion, we have to be very careful with boundary conditions. If we consider a system of $N$ fermions on a ring, we actually have to consider 4 different cases: even/odd fermion parity and periodic/anti-periodic boundary conditions. Those 4 topological sectors will be mapped to different sectors on the spin side (as is well-known from \eg~the exact solution of the Ising model in a transverse magnetic field). For each of those cases, we will define a different graded MPO $U^i$ that is build

\begin{table}[]
\centering
\[\begin{array}{|c|c|cc|cc|cc|c|} \hline
& \text{Bosonized} & \multicolumn{2}{|c|}{\text{Fermions}} & \multicolumn{2}{|c|}{\text{Spins}} & \multicolumn{2}{|c|}{\text{Dual Spins}} & \text{Unified} \\
i & \text{BC} & \bigotimes_n {\sigma_f^z}_n & \text{Twist} & \bigotimes_n {\sigma^x}_n & \text{Twist} & \bigotimes_n {\sigma^z}_n & \text{Twist} & \text{BC} \\ \hline
1 & \sigma^z & +1 & I & -1 & I & +1 & \sigma^z & I \\
z & I & +1 & \sigma^z_f & +1 & I & +1 & I & I \\
x & \sigma^x & -1 & I & +1 & \sigma^x & -1 & I & \sigma^x \\
y & \sigma^x \sigma^z & -1 & \sigma^z_f & -1 & \sigma^x & -1 & \sigma^z & \sigma^x \\
\hline\end{array}\]
\caption{Summary of how different topological sectors are mapped by $U^i$, the composition of a Jordan-Wigner transformation and a Kramers-Wannier duality (Fermions and Spins columns), and by $U_\text{KW}^{\phi(i)}$, a matched Kramers-Wannier duality (Spins and Dual Spins columns). The mapping of topological sectors by their composition, i.e., the usual Jordan-Wigner transformation, is hence given by the Fermions and Dual Spins columns with the combined bosonic boundary condition (BC) in the final column.}
\label{tab:JWtable}
\end{table}

%Note that the charge on the fermion side is given by the expectation value of the operator $\otimes_{l=1}^N \sigma^z_{f,l}$ and on the spin side by the expectation value of $\otimes_{l=1}^N \sigma^x_{l}$.
%The application of an odd operator on the fermion side results in the switch of boundary conditions. \Eg,
%    \sigma^x_{f,n} U^i_\text{JW} = U^{\xi(i,x)}_\text{JW} \left(\otimes_{l=1}^{n-1}\sigma^x_{l}\right) \sigma^z_{n-1}
%with $\xi(i,x)$ switching the topological sector according to the rule $\sigma^{\xi(i,x)}\simeq \sigma^i\sigma^x$. Similarly, $\sigma^{\xi(i,y)}\simeq \sigma^i\sigma^x$, $\sigma^{\xi(i,z)}= \sigma^i$, and $\sigma^{\xi(i,1)}= \sigma^i$.

The usual JW transformation is obtained by doing an extra Kramers-Wannier transformation on the spin side. This KW transformation is of the same form as \refeq{eq:JW1f},
\begin{equation}
U^i_\text{KW} = \diagram{jw}{22} \, ,
\end{equation}
\begin{equation}
\sum_{\alpha \beta s} \delta_{\alpha s} \delta_{\beta s} \ket{\alpha} \ket{s} \bra{\beta} = \diagram{jw}{23} \qquad \text{and} \qquad
\sum_{\alpha+\beta=s\mod{2}} \ket{\alpha} \bra{s} \bra{\beta} = \diagram{jw}{24}
\end{equation}
are fully bosonic now but have the same properties as in \refeq{eq:propghz}. One can again work out how charge sectors on the top and bottom side of $U^i_KW$ are mapped onto each other with the $_2$ charge on top, i.e., for the dual spins, being the eigenvalue of $⊗_n σ^z_n$ and the non-trivial dual twists being $σ^z$. The results are summarized in Table \ref{tab:JWtable}. Note that the KW duality essentially has the effect of swapping the charge and twist columns. Furthermore, note that we constructed the table for $U^ϕ(i)_KW$ rather than for $U^i_KW$. Indeed, this matching of the boundary conditions makes sure that when gluing the KW on top of $U^i$, the spin charge sectors are the same and the composite transformation, i.e., the JW, will not reduce to zero. Moreover, one can immediately read off the mapping of the topological sectors for the composite JW MPO. Regarding the tensors in this MPO, we have that
\begin{equation}
\diagram{jw}{25} = \diagram{jw}{26} \quad \Rightarrow \quad \diagram{jw}{27} = \diagram{jw}{28} = \diagram{jw}{29}\, ,
\end{equation}
\begin{equation}
\diagram{jw}{30} = \sum_s \sum_{\alpha+\beta=f\mod{2}} \delta_{sf} \ket{\alpha} |f) \bra{s} \bra{\beta}\, .
\end{equation}
The bond dimension thus reduces to 2. For the combined boundary conditions we obtain
\begin{equation}
\diagram{jw}{31} = \diagram{jw}{33}= \diagram{jw}{32}
\end{equation}
\begin{equation}
\diagram{jw}{34} = \diagram{jw}{36}= \diagram{jw}{35}
\end{equation}
so that
\begin{equation}\label{eq:JW2}
\begin{split}
U_\text{JW}^i &= \diagram{jw}{37}\\[8pt]
&= \sum_{f_1,f_2,\cdots}\rm{Tr}\left(\tilde{G}^{f_1}\tilde{G}^{f_2}\cdots \tilde{G}^{f_N}\sigma^{\chi(i)}\right)|f_1)|f_2)\cdots |f_N)\bra{f_N}\cdots \bra{f_2} \bra{f_1} \, ,
\end{split}
\end{equation}
with $χ(1)=χ(z)=1$ and $χ(x)=χ(y)=x$. As opposed to the combination of JW and KW, exactly the same MPO (i.e., with the same BC) deals with both the periodic and anti-periodic case in each charge sector. On the other hand, putting a $σ^z$ ($σ^y$) BC yields zero in the even (odd) sector. Therefore, we could just as well take sum, returning an MPO with open boundary conditions $|0⟩ ⟨0|$ in the even sector and $|1⟩ ⟨0|$ in the odd. Then the JW transformations, both for the even and odd case, can be unified in a single MPO with boundary term $(|0⟩+ |1⟩)⟨0|$. This is the reason that one typically just talks about one unique JW transformation.

Clearly, the $_2$-graded tensor network language yields a pretty natural framework for formulating and understanding the intricacies of the JW transformation. What is more, it can also be naturally extended to higher dimensions by explicitly retaining the fermionic nature of the virtual legs. Again, we refer to Ref.~\cite{Obrien} for more details. The question whether it is better to use graded tensor networks or to use standard tensor network methods on a JW transformed Hamiltonian also imposes itself now. The answer is certainly not unequivocal, but will depend on the questions that one wishes to ask. For example, the nature of the excitations might change drastically under duality transformations: local excitations can become topological domain-wall like excitations and vice-versa. Also, local symmetries in one representation can become non-local in the other one. %; this happens e.g. with the Hubbard model, where it is not obvious how to preserve all symmetries (particle-hole + SU$(2)$ + U$(1)$) under JW transformations.
Equally so, one has to be very careful in labelling the different entanglement spectra, as duality transformations permute the different labels (as examplified in the table above).

\section{Fermionic PEPS}

Fermionic PEPS (fPEPS) represent the 2D analogue of the fermionic MPS discussed in the previous sections. Here, we explain how one constructs fermionic PEPS within the framework of $_2$-graded Hilbert spaces. Next, contraction of these tensor networks is considered with the aim to calculate expectation values. Both the boundary MPS (also referred to as MPS-MPO) approach \cite{fishman2018, vanderstraeten2022} and CTMRG \cite{nishino1996,nishino1997,orus2009} are considered.

\subsection{General form}
\label{sec:fpeps}

The conceptual construction of a 2D fermionic PEPS strongly resembles that of the 1D fermionic MPS. To each lattice site we associate virtual fermionic degrees of freedom but now in both spatial directions. For a simple square lattice, we thus have a pair of virtual degrees of freedom per lattice site in both the horizontal and the vertical direction. Next, neighboring virtual fermions are brought into a maximally entangled state,
\begin{align}
\ket{\Psi_\mathrm{virt}} &= \prod_{e_h} \left( 1 + a_{e_h,l}\dag
a_{e_h,r}\dag \right) \prod_{e_v} \left( 1 + a_{e_v,d}\dag a_{e_v,u}\dag \right) \ket{\Omega_\mathrm{virt}} =\diagram{peps}{1}\, ,
\end{align}
after which we introduce a linear, parity-preserving map, locally projecting the four virtual modes per site onto a single physical mode,
\begin{equation}
\begin{split}
\mathcal{A}_j = \left(A_j\right)_{\alpha_1 \alpha_2 \sigma \alpha_3 \alpha_4} (a_{j,l})^{\alpha_1} (a_{j,d})^{\alpha_2} (a_j\dag)^{\sigma} (a_{j,r})^{\alpha_3} (a_{j,u})^{\alpha_4} \, .
\end{split}
\end{equation}
As was the case in 1D, each virtual link amounts to a Kronecker delta connecting neighboring virtual indices. However, as the projectors are no longer placed along a line, fermionic parity factors will appear as well. Still, the final result can be reinterpreted as a contraction of local $ℤ_2$-graded tensors,
\begin{equation}
\begin{split}
\ket{\Psi(A)} = \mathcal{C}_\mathrm{virt} \Biggl( \prod_j \, \Biggl(\left(A_j\right)_{\alpha_j\beta_j \gamma_j \varepsilon_j \delta_j} \ket{\alpha_j} \ket{\beta_j} \ket{\gamma_j} \bra{\varepsilon_j} \bra{\delta_j}  \Biggr) \Biggr) ,
\end{split}
\end{equation}
where the fPEPS tensors are defined as
\begin{equation}
\diagram{peps}{2} = \left(A_j\right)_{\alpha \beta \gamma \delta \varepsilon} \ket{\alpha} \ket{\beta} \ket{\gamma} \bra{\varepsilon} \bra{\delta}
\end{equation}
with $|γ⟩$ labeling the physical leg. Assuming translation invariance, the fPEPS construction can thus be summarized as
\begin{equation}
\ket{\Psi(A)} = \diagram{peps}{3} \, .
\end{equation}
As was the case for fMPS, both the physical and virtual Hilbert spaces can be enlarged. Taking $f$ fermions on the physical level for instance amounts to a local Hilbert space of dimension $d=2^f$ whereas for the virtual spaces arbitrary natural bond dimensions $D$ are possible, resulting in a fPEPS object living in a $dD^4$-dimensional, graded Hilbert space. Also note that upon contraction of these objects, fermionic minus signs due to supertraces are inevitable in contrast to what we observed in 1D for fMPS.

Analogous to the 1D case, fermionic PEPS have gauge freedom due to their specific virtual structure, \ie, an operation
\begin{equation}
\diagram{peps}{4} \rightarrow \diagram{peps}{5}
\end{equation}
leaves the fPEPS unchanged but now different gauge transforms can be applied in both directions. However, unlike in 1D, it is not possible to use this gauge transforms to impose canonical forms that significantly simplify the calculation of expectation values and ensure unit normalisation.
% A canonical form, as found in 1D, cannot be derived in two dimensions because simple entanglement cuts cannot be preformed due to the loops in the diagrams.
Just like their bosonic counterparts, fermionic PEPS thus require approximate methods for their contraction.

\subsection{Boundary MPS}
\label{sec:boundarymps}

One possible approach to contract an fPEPS relies on the determination of boundary MPS. We will first illustrate this for the calculation of the norm of a fPEPS, which can be expressed as
\begin{equation}
\braket{\Psi(A)|\Psi(A)} = \diagram{peps}{8} \, ,
\end{equation}
where the double-layer tensor $O$ is defined by
\begin{equation}
\diagram{peps}{7} = \diagram{peps}{6}
\label{doublelayertensor}
\end{equation}
and where
\begin{equation}
\diagram{peps}{115} = \diagram{peps}{114}
\end{equation}
is just an internally permuted version of the conventional $A$. Again, all the physical legs in $|Ψ(A)⟩$ point outward so that no $P$ tensors appear when taking the norm by contracting with its conjugate. In order to contract this tensor network, we calculate the leading top and bottom eigenvector, also called fixed points, of the horizontal linear transfer matrix,
\begin{equation}
T(O) = \diagram{peps}{9} \, .
\end{equation}
Once computed, both can be contracted, yielding the norm of $|Ψ(A)⟩$.

To obtain the fixed points, we approximate them by a uniform MPS with a certain boundary bond dimension $χ$ and with two ``physical'' legs instead of one. \Ie, we propose
\begin{equation}
\ket{\zeta(M)} = \diagram{peps}{10}
\end{equation}
for the top boundary state. The relevant eigenvalue equation then becomes
\begin{equation}
\diagram{peps}{12} \approx \Lambda_t \diagram{peps}{10} \, .
\label{eq:topfixedpoint}
\end{equation}
Note that this is not equivalent to applying the corresponding operator, $T |ζ(M)⟩ = Λ_t |Ψ(M)⟩$. Indeed, in the latter $P$ tensors are placed on all of the upward pointing legs of $|ζ(M)⟩$. Therefore, we redefine our transfer operator by
\begin{equation}
\begin{split}
\tilde{T}(\tilde{O}) &= \diagram{peps}{14} = \diagram{peps}{13}
\end{split}
\end{equation}
where the circular tensors correspond to $P$. Application of this slightly altered operator to $|ζ(M)⟩$ then amounts to the desired eigenvalue equation. Similarly, the solution of the bottom fixed-point equation could be approximated by an MPS having the structure of the conjugate of $|ζ(M)⟩$, \ie,
\begin{equation}
\bra{\zeta(N)} = \diagram{peps}{15} \, .
\end{equation}
The appropriate eigenvalue equation then becomes $⟨ζ(N)| T̃ = Λ_b ⟨ζ(N)|$, corresponding to
\begin{equation}
\diagram{peps}{16} = \Lambda_b \diagram{peps}{15} \, .
\label{eq:botfp}
\end{equation}
A useful feature of these definitions of the fixed-points is that the bottom eigenvalue equation can be rewritten in the form of the top one. Taking the adjoint of the bottom equation, we obtain $T̃^†|ζ(N)⟩ = Λ_b |ζ(N)⟩$, \ie, the top equation but now for $T̃^†$. Indeed, $|ζ(M)⟩$ is the right eigenvector of $T̃$ (as an operator) while $⟨ζ(N)|$ is the left eigenvector. Writing out the resulting equation in a diagram, we obtain
\begin{equation}
\diagram{peps}{18} = \astt{\Lambda_b} \diagram{peps}{19} \, .
\end{equation}
We conclude that the bottom fixed-point equation can be transformed to a top fixed-point equation as in \refeq{eq:topfixedpoint} but with $M →N$, $Λ_t →Λ_b$ and
\begin{equation}
\begin{split}
\diagram{peps}{7} \rightarrow \diagram{peps}{24} \, .
\end{split}
\label{eq:pepsdagger}
\end{equation}
Alternatively, we could absorb the parities on the vertical legs in the $N$ tensors. This does remove quite some parities from the following equations for the general case of non-hermitian $T̃$ where $M ≉N$. However, for hermitian $T̃$, one would have to add these $P$ tensors again instead of just approximating the $N$ by $M$ tensors (see \refeq{eq:doublelayerherm}). Therefore, we do not absorb the parities here. With $Λ_t$ being the dominant eigenvalue of $T̃$ and $Λ_b$ of $T̃^†$, we also immediately have that $Λ_t ≈Λ_b$ and denote both with $Λ$.  To solve the top fixed-point equations for $M$ and $N$, the VUMPS algorithm can be applied again but now for MPS tensors with two physical legs. This introduces some slight alterations to the method that we work out in \refapp{app:vumpsmpo2}.

%While this seems to have done the trick, writing out the resulting equation in a diagram,
%	\diagram{peps}{18} = \astt{\Lambda_b} \diagram{peps}{19} \, ,
%shows that this is still not formally identical to \refeq{eq:topfixedpoint}. Indeed, even though the $P$ tensors and the twisted legs can be absorbed in the MPO tensor, the virtual arrows of this operator are still pointing in the opposite way. Therefore, we flip them by using an arbitrary flipper, where
%	\begin{split}
%		\diagram{peps}{20} &= \diagram{peps}{21}\\
%		\diagram{peps}{22} &= \diagram{peps}{23}
%	\end{split}\, .
%	\label{eq:flips}
%We then conclude that the bottom fixed-point equation can be transformed to a top fixed-point equation as in \refeq{eq:topfixedpoint} but with $M \rightarrow N$, $\Lambda_t \rightarrow \astt{\Lambda_b}$ and
%	\begin{split}
%		\diagram{peps}{7} \rightarrow \diagram{peps}{24} \, .
%	\end{split}
%	\label{eq:pepsdagger}

After computing the fixed-points in their center-site gauge with the VUMPS algorithm, the norm and local expectation values of fermionic PEPS can easily be computed. For instance, sequential application of the eigenvalue equation to the overlap between the boundary states yields
\begin{equation}
\begin{split}
\braket{\Psi(A)|\Psi(A)} &= \Lambda^{|\ZZ|} \braket{\zeta(N)|\zeta(M)} = \Lambda^{|\ZZ|} \diagram{peps}{25} \, .
\end{split}
\end{equation}
Indeed, with every application of \refeq{eq:topfixedpoint} or \refeq{eq:botfp}, an additional $T(O)$ layer appears, finally yielding the expression for the norm up to a (possibly infinite) factor. Essentially, this factor is equal to $Λ^N_v$ with $N_v$ the infinite number of vertical layers. Furthermore, \refapp{app:vumpsmpo2} learns that $Λ= λ^N_h$ with $N_h$ the (infinite) number of sites in the horizontal direction and $λ$ the dominant eigenvalue obtained in the determination of the environments (\refeq{eq:app:envvumpsmpo2}). We thus find that $⟨Ψ(A)|Ψ(A)|=⟩ λ^N_h N_v ⟨ζ(N)|ζ(M)|$⟩. For the remaining MPS overlap, we determine left and right environments,

\begin{equation}
\diagram{peps}{26} = \mu \diagram{peps}{27} \quad \text{and} \quad \diagram{peps}{28} = \mu \diagram{peps}{29}
\end{equation}

that we normalize as

\begin{equation}
\diagram{peps}{30} = 1
\end{equation}

with $C_M$ ($C_N$) the canonical $C$ tensor for the top (bottom) boundary MPS so that for the norm $\braket{\Psi(A)|\Psi(A)} = \lambda^{N_h N_v} \mu^{N_h}$. This shows that normalizing an arbitrary PEPS is not as straightforward as for MPS. However, one typically does not require a normalized state for the calculation of local expectation values. Indeed, consider a strictly local operator $W$, then

\begin{equation}
\braket{\Psi(A)|W|\Psi(A)} = \diagram{peps}{31}
\end{equation}

\begin{equation}
\diagram{peps}{32} = \diagram{peps}{33}\, .
\end{equation}

Using the fixed points, this reduces to

\begin{equation}
\braket{\Psi(A)|W|\Psi(A)} =\lambda^{N_h(N_v-1)} \diagram{peps}{34}
\end{equation}

where environments can be defined as

\begin{equation}
\diagram{peps}{35} = \mu \diagram{peps}{36} \quad \text{ and } \quad \diagram{peps}{37} = \mu \diagram{peps}{38}
\label{eq:environments}
\end{equation}

\begin{equation}
\diagram{peps}{39} = 1
\end{equation}

so that

\begin{equation}
\braket{\Psi(A)|W|\Psi(A)} =\lambda^{N_h(N_v-1)} \mu^{N_h-1} \diagram{peps}{40} \, .
\end{equation}

However, as we did not explicitly normalize the PEPS, we are rather interested in

\begin{equation}
\frac{\braket{\Psi(A)|W|\Psi(A)}}{\braket{\Psi(A)|\Psi(A)}}
\end{equation}

where the norm can also be expressed as

\begin{equation}
\begin{split}
\braket{\Psi(A)|\Psi(A)} &= \lambda^{N_h(N_v-1)} \diagram{peps}{41} = \lambda^{N_h(N_v-1)} \mu^{N_h}
\end{split}
\end{equation}

so that

\begin{equation}
\frac{\braket{\Psi(A)|W|\Psi(A)}}{\braket{\Psi(A)|\Psi(A)}} = \frac{1}{\mu} \diagram{peps}{40} \, .
\end{equation}

In this way, all the possibly infinite prefactors have been eliminated. We conclude that calculating local expectation values poses no problem for fermionic PEPS. Similarly, two-point correlation functions in the horizontal and vertical direction can be calculated efficiently in this way (see sec:2d).

The contraction of PEPS becomes particularly manageable when the linear transfer operator is hermitian. Indeed, when $\tilde{T}=\tilde{T}^\dagger$, the bottom fixed point is equal to the adjoint of the top fixed point and both have the same real eigenvalue. As a result only a single VUMPS run has to be performed to obtain both. Furthermore, the norm can be expressed as

\begin{equation}
\begin{split}
\braket{\Psi(A)|\Psi(A)} &= \lambda^{N_h N_v} \braket{\Psi(M)|\Psi(M)}\\
&= \lambda^{N_h N_v} \diagram{peps}{42} = \lambda^{N_h N_v}\, .
\end{split}
\end{equation}

with $\lambda \in \RR$ and positive so that it can be regarded as a norm density. Hence, dividing each $A$ tensor with $\sqrt{\lambda}$ immediately yields a normalized PEPS. To impose the hermiticity on an arbitrary fPEPS, one could use eq:pepsdagger. When

\begin{equation} \label{eq:doublelayerherm}
\begin{split}
\diagram{peps}{7} = \diagram{peps}{24} \, ,
\end{split}
\end{equation}

$\tilde{T}$ will be hermitian. Writing out this equality in terms of the PEPS tensor and its conjugate, we get

\begin{equation}
\begin{split}
\diagram{peps}{6} = \diagram{peps}{113} \, .
\end{split}
\end{equation}

Note that on the left-hand side of the equation, the upper left legs belong to a single tensor whereas this is not the case on the right-hand side. Therefore, we introduce the following unitary operations in the latter,

\begin{equation}
\begin{split}
\diagram{peps}{113} \quad \rightarrow \quad \diagram{peps}{112} \, .
\end{split}
\end{equation}

Indeed, making a chain, the swaps and flippers are undone so that the sufficient condition remains satisfied. Furthermore, we can express it on the level the PEPS tensor as the combination of

\begin{equation} \label{eq:pepsherm}
\begin{split}
\diagram{peps}{4} = \diagram{peps}{43}
\end{split}
\end{equation}

\begin{equation} \label{eq:pepsherm}
\begin{split}
\diagram{peps}{115} = \diagram{peps}{45} \, .
\end{split}
\end{equation}

Rewriting $A^\ast$ in terms of $\thiccbar{A}$, one can easily show that both equations are equivalent. Imposing eq:pepsherm, we are thus guaranteed to have a hermitian $\tilde{T}$.
By additionally integrating gauge transforms, one turns this sufficient condition into a necessary one.

Using the techniques above, one could for instance calculate the energy density of a fPEPS for a certain Hamiltonian. Combining this with automatic differentiation <cit.>, the gradient of the energy density can be computed, thus allowing a variational optimization of the states as previously demonstrated for CTMRG <cit.>.

§.§ Corner Transfer Matrix Renormalization Group

The corner transfer matrix renormalization group algorithm <cit.> represents a complementary alternative to the boundary MPS approach for contracting PEPS. At its core is the assumption that the double layer of a PEPS and its conjugate,

\begin{equation}
\braket{\Psi(A)|\Psi(A)} = \scaleddiagram{peps}{46}{0.3} \, ,
\end{equation}

can be represented as

\begin{equation}
\diagram{peps}{47} \, ,
\label{eq:ctmrg}
\end{equation}

where the corner transfer matrices (CTMs) $\{C_1,C_2,C_3,C_4\}$ approximate the infinite corners of the network while the half rows and columns are approximated by $\{T_1,T_2,T_3,T_4\}$. The virtual legs regulating the size of these transfer matrices correspond to spaces of dimension $\chi$. As a result, a local expectation value can be computed by performing the same contraction as in eq:ctmrg but with $O$ replaced by $O_W$. To determine the $C$ and $T$ tensors, the CTMRG algorithm in its most general form consists of three steps respectively applied in the left, right, upward and downward direction:

* insertion: the network from eq:ctmrg is enlarged, here in the left direction (, a “left move”) by inserting a new column containing $T_1$, $O$ and $T_3$,

\begin{equation}
\diagram{peps}{48} \, .
\end{equation}

* absorption: the inserted tensors are contracted with their neighboring $C$ and $T$ yielding the bigger

\begin{equation}
\diagram{peps}{49} = \diagram{peps}{50}\, ,
\end{equation}

\begin{equation}
\diagram{peps}{51} = \diagram{peps}{52}
\end{equation}

\begin{equation}
\diagram{peps}{53} = \diagram{peps}{54}\, .
\end{equation}

* renormalization: $\tilde{C}_1$, $\tilde{T}_4$ and $\tilde{C}_4$ are truncated back to a virtual space of the original size by applying the projectors

\begin{equation}
\diagram{peps}{55} \qquad \text{and} \qquad \diagram{peps}{56} \qquad \text{with} \qquad
\diagram{peps}{57} \approx \diagram{peps}{58}\, ,
\label{eq:projectorcondition}
\end{equation}

yielding the updated transfer matrices

\begin{equation}
\diagram{peps}{59} = \diagram{peps}{60}\,, \quad
\diagram{peps}{61} = \diagram{peps}{62} \, , \quad
\diagram{peps}{63} = \diagram{peps}{64} \, .
\end{equation}

Cycling through the four directions while repeating these steps and normalizing the transfer matrices, the CTMRG algorithm will converge towards a fixed point, invariant under new insertions. As a result, the norm of the PEPS as expressed in eq:ctmrg can equally well be expressed as

\begin{equation}
\diagram{peps}{65} \quad , \quad
\diagram{peps}{66} \quad \text{or} \quad
\diagram{peps}{67}\, .
\end{equation}

Multiplying the former and eq:ctmrg and dividing by the latter two then yields a quantity whose difference between consecutive iterations can be used as an error measure.

Determination of a proper choice for the projector in the renormalization step is crucial for the stability and efficiency of the algorithm. Many different methods for obtaining these projectors have been proposed. Here, we focus on the one introduced in Ref. <cit.> and reformulated in Ref. <cit.>, a recurring choice in many state of the art CTMRG codes. For the left move it is constructed as follows.

* First, define the half-system transfer matrices $U$ and $D$,

\begin{equation}
\diagram{peps}{68} = \, \diagram{peps}{69},
\end{equation}

and perform an SVD of their left-side contraction,

\begin{equation}
\diagram{peps}{70} \approx \, \diagram{peps}{71},
\end{equation}

where the approximate equality indicates that we only keep the highest singular values to truncate back to the desired bond dimension.

* Construct the projectors

\begin{equation}
\qquad \diagram{peps}{72} = \, \diagram{peps}{73}
\end{equation}

\begin{equation}
\qquad \diagram{peps}{74} = \, \diagram{peps}{75} \, ,
\end{equation}

where $\Sigma_L^+$ is the pseudoinverse of $\Sigma_L$ and where the circular tensors are parity tensors. Indeed, these are necessary to make use of the left (right) isometry of $U_L$ ($V_L$) so that eq:projectorcondition holds. The only alteration to the CTMRG algorithm for fermions thus consists of adding $P$ tensors to the SVD decompositions in the renormalization step.

§ APPLICATIONS IN 2D

To benchmark the fermionic TNS methods in 2D, we reproduce correlation functions of Gaussian fermionic TNS (GfTNS) <cit.> approximating the ground state of a chiral model. Other studies have found that in chiral phases, variationally optimized PEPS tend to reproduce the exponentially decaying correlation functions up to a certain bulk correlation length, regulated by the PEPS bond dimension $D$, after which a, possibly polynomial, long-range tail takes over <cit.>. Here, we intend to verify this behavior with GfTNS for the p-wave superconductor. Subsequently, generic fPEPS algorithms are used to reproduce the Gaussian correlation functions.

§.§ Long-range correlations in chiral superconductors

Consider the p-wave superconductor,

\begin{align}
H = &- t\sum_\mathbf{n} \left( a_\mathbf{n}^\dagger a_{\mathbf{n}_\rightarrow}+ a_\mathbf{n}^\dagger a_{\mathbf{n}_\uparrow} +h.c.\right) -\mu \sum_\mathbf{n} a_\mathbf{n}^\dagger a_{\mathbf{n}} \nonumber - \Delta \sum_\mathbf{n} \left(a_\mathbf{n}^\dagger a_{\mathbf{n}_\rightarrow}^\dagger +i a_\mathbf{n}^\dagger a_{\mathbf{n}_\uparrow}^\dagger + h.c.\right) \, ,
\label{eq:peps:pwave}
\end{align}

where $\mathbf{n}_{\rightarrow(\uparrow)}$ corresponds to the nearest neighbor of site $\mathbf{n}$ on the right (upper) side in a square lattice. More specifically, we examine the point where $\Delta/t=1$ and $\mu/t = 2$, , a gapped and chiral topological superconductor with Chern number $C=1$. As this model is free it can be solved by a Bogoliubov transformation. The quadratic nature also justifies the use of GfTNS to efficiently approximate the ground state,  by minimizing their energy density with a gradient descent method. Indeed, for GfTNS these local observables can be calculated analytically using correlation matrices rather than via approximate environments, allowing for a simpler and faster optimization <cit.>. However, the restriction to GfTNS also comes at a cost. Calculating local observables for a 2D translation invariant GfTNS in the thermodynamic limit for instance requires the computation of continuous $\mathbf{k}$-space integrals. In more than one spatial dimension this can no longer be done analytically and therefore one has to work with finite lattices containing $N_s$ sites and (anti-)periodic boundary conditions, thus replacing the integrals with finite sums. More importantly, the no-go theorem of Dubail and Read applies, stating that non-singular GfTNS cannot realize chiral features and therefore have a trivial Chern number <cit.>. As explained in Ref. <cit.>, optimizing GfTNS while starting from a random and thus non-singular initial guess, these states will never become truly singular since this requires a strict fine-tuning of their parameters. The approximate GfTNS will therefore never become truly chiral, and ultimately exponential decay will take over again after the polynomial decay at intermediate distances. Indeed, a more detailed examination shows that the Berry curvature generated all through the Brillouin zone is compensated in a very small region around $\mathbf{k}=0$. In real space, this results in a polynomial tail for the correlations up to a corresponding (large) length scale after which the final exponential decay takes over. The compensating region in $\mathbf{k}$-space was found to shrink by increasing the bond dimension as well as the finite system size $N_s$ used during the GfTNS optimization. Indeed, a higher $N_s$ samples more points close to the zone center so that the cost function tries to shrink this region further. Increasing $N_s$ and $D$, the length scale for the onset of the exponential tail hence increases. All of these features can be confirmed in fig:tsc for $||\Gamma_{(0,0),(x,0)}||$, the Frobenius norm of the real-space correlation matrix between sites $(0,0)$ and $(x,0)$,

\begin{equation}
\Gamma^{kl}_{(0,0),(x,0)} = \frac{i}{2} \braket{\Psi|\left[c^k_{(0,0)},\, c^l_{(x,0)}\right]|\Psi}
\end{equation}

with $c^{1}_\mathbf{n} = a^\dagger_\mathbf{n} + a_\mathbf{n}$ and $c^{2}_\mathbf{n} = -i( a^\dagger_\mathbf{n} - a_\mathbf{n})$ the Majorana operators. At elevated $D$ and $N_s$ the exact, exponentially decaying region is followed by a polynomial regime. However, fig:tsc also shows that this polynomial regime is followed again by a second exponential decay. When $D$ and $N_s$ are too low, the polynomial region shrinks or completely disappears. If $D=2$, for instance, the two exponential regimes immediately follow each other, independent of $N_s$. For $D=4$ on the other hand, $N_s=50^2$ results in a limited polynomial regime quickly followed by the second exponential decay while $N_s=1000^2$ results in a significantly longer polynomial tail. If $D=8$, the second exponential region is never even probed for separations $x$ ranging up to $500$ sites. The notion that chiral features in PEPS require long-range tails in their correlation functions is thus sustained and even explained in the sense that it is a consequence of the no-go theorem. However, for approximate, non-fine-tuned GfTNS, this polynomial region also has a finite length scale after which a second exponential regime sets in.

Semi-logarithmic (top panels) and double logarithmic (bottom panels) plots of the real-space correlations, more specifically $||\Gamma_{(0,0),(0,x)}||$, the Frobenius norm of the real-space correlation matrix between sites $(0,0)$ and $(x,0)$) for the chiral p-wave superconductor at $\Delta/t=1$ and $\mu/t = 2$. Exact results obtained by means of a Bogoliubov transformation are compared to those for GfTNS (circular markers) evaluated on $3000\times3000$ lattice but optimized on both small ($N_s=50^2$, left panels) and large lattices ($N_s=1000^2$, right panels) and with bond dimensions $D=2,4,8$. For the latter, polynomial tails are only resolved at elevated $D$ and $N_s$. These Gaussian results were reproduced with generic fPEPS methods and the $\times$ ($+$) sign markers correspond to values obtained with boundary MPS (CTMRG) where $\chi = 80, 200, 400$ to contract the infinite-size fPEPS with $D=2,4,8$ obtained from the optimized GfTNS tensors.

The peculiar correlation functions described above represent an interesting test case for the contraction methods outlined in the previous sections. Indeed, after casting the approximate GfTNS in a generic fPEPS format as described in app:gftns, these methods should reproduce the Gaussian results. Normalized correlations in the horizontal direction like $\braket{a^\dagger_{(0,0)} a^\ph_{(x,0)}}$ and $⟨a_(0,0) a_(x,0)|$⟩ can be calculated using the boundary MPS by contracting

\begin{equation}
\diagram{peps}{76} \, ,
\end{equation}

\begin{equation}
\diagram{peps}{77} = \diagram{peps}{78} \qquad \text{and} \qquad \diagram{peps}{79} = \diagram{peps}{80}
\end{equation}

contain the physical creation and annihilation operators as in eq:kitaevmpo,

\begin{equation}
\begin{split}
\diagram{peps}{81} &= \ket{1} \bra{0} (1| \qquad
\diagram{peps}{82} = \ket{0} \bra{1} (1| \\
\diagram{peps}{83} &= |1) \ket{1} \bra{0} \qquad
\diagram{peps}{84} = |1) \ket{0} \bra{1} \, .
\end{split}
\end{equation}

For the normalized correlations, we thus obtain

\begin{equation}
\braket{a^{(\dagger)}_{(0,0)} a^\ph_{(x,0)}} = \frac{1}{\mu^{x+1}} \diagram{peps}{85}
\label{eq:corrsmpsmpo}
\end{equation}

\begin{equation}
\diagram{peps}{86} = \diagram{peps}{87}
\end{equation}

and $\mu$, $G^l$ and $G^r$ as in eq:environments.

Similarly, the converged CTMs and half row/column transfer matrices yield the unnormalized correlations

\begin{equation}
\braket{a^{(\dagger)}_{(0,0)} a^\ph_{(x,0)}} = \diagram{peps}{88}
\end{equation}

\begin{equation}
\diagram{peps}{89} = \diagram{peps}{90}
\end{equation}

the horizontal transfer matrix. An efficient and stable way to calculate these correlations  proceeds by setting

\begin{equation}
\diagram{peps}{91} = \diagram{peps}{92} \quad , \quad \diagram{peps}{93} = \diagram{peps}{94}
\end{equation}

\begin{equation}
\diagram{peps}{95} = \mu \diagram{peps}{96} \, .
\end{equation}

Thus assuming that $G^l$ $(G^r)$ approximates the dominant left (right) eigenvector upon contraction with $E_h$ for the converged CTMRG environments, the normalized correlations can be expressed as

\begin{equation}
\braket{a^{(\dagger)}_{(0,0)} a^\ph_{(x,0)}} = \frac{1}{\mu^{x+1}} \diagram{peps}{97} \, .
\end{equation}

Form the $\braket{a^\dagger_{(0,0)} a^\ph_{(x,0)}}$ and $⟨a_(0,0) a_(x,0)|$⟩ correlation functions we compute $||\Gamma_{(0,0),(0,x)}||$. Comparing the results to the Gaussian correlation functions in fig:tsc confirms that both generic contraction methods yield a satisfactory reproduction and this with manageable environment bond dimensions $\chi = 80, 200, 400$ for $D=2,4,8$.

§ CONCLUSIONS

We have shown that fermions can be directly built into tensor network states without a need to keep track of leg crossings explicitly. Indeed, constructing fermionic tensors as elements of graded Hilbert spaces, we derived a methodology that takes into account their anti-commuting nature within the tensor structure itself, as well as in other tensor operations such as contractions and decompositions. The resulting formalism can be diagrammatically represented, with fermionic tensors having a specific orientation for each leg, in line with common practices to incorporate global symmetries. As a result, it can be readily integrated in standard tensor network toolboxes with symmetries.

Once built up, we have used the formalism to construct fermionic versions of MPS and PEPS, and discussed standard algorithms like DMRG, tangent-space methods and CTMRG. Highlighting the differences with their bosonic counterparts in the form of additional parity tensors, we have fermionized important elements of the traditionally bosonic TNS toolbox. We have benchmarked these algorithms against interesting test cases like the Kitaev chain where we resolve the gapless edge modes, the 1D Hubbard model exhibiting spin-charge separation and the chiral p-wave superconductor in 2D with characteristic long-range tails in the TNS approximations of real-space correlations.

Though this work establishes the main rationale behind a consistent fermionic formalism, it only serves as a starting point. Numerous concepts and algorithms within the extensive tensor network literature have yet to be explored. RG-like procedures that alter the lattice geometry, such as the tensor renormalization group (TRG), and their connection to spin structures, as a condition for consistent spinful fermionic theories in the continuum, remain unexplored.
The $\ZZ_2$ grading underlying our formalism can also be extended to $\ZZ_p$, effectively incorporating parafermions <cit.>. This approach has already been followed to explore potential gapped phases but one could again aim to extend the TNS toolbox based on this more complicated grading <cit.>.
Another direction for future research involves a detailed comparison of the advantages of working directly with fermions versus employing a Jordan-Wigner transformation first <cit.>. Especially in higher dimensions, the interplay between more involved global symmetries and periodic as well as more exotic boundary conditions remains to be explored. Nonetheless, the graded formalism presented here appears to be the most natural language to address these questions.

§ ACKNOWLEDGEMENTS

We acknowledge useful discussions with Maarten Van Damme and Laurens Lootens.
This research is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program Grant agreements No. 101076597 - SIESS (NB), the Research Foundation Flanders, EOS (grant No. 40007526), IBOF (grant No. IBOF23/064), and BOF-GOA (grant No. BOF23/GOA/021). LB is supported by a PhD fellowship 11H7223N from the Research Foundation Flanders. LV is supported by the FNRS (Belgium).

§ FERMIONIC TANGENT SPACE METHODS

§.§ Tangent-space projector for fermionic MPS

We derive an expression for $P_{\ket{\Psi(A)}}$, the tangent space projector in and orthogonal to a point $\ket{\Psi(A)}$ on the uniform fMPS manifold. Therefore, consider an arbitrary translation invariant state $\ket{\chi}$. Its projection onto the tangent space, $\ket{\Phi(B(X);A)}=P_{\ket{\Psi(A)}} \ket{\chi}$, can be parametrized by a tensor $X$ as in eq:cantangentleft where

\begin{align}
X &= \min_X \left\lVert \ket{\chi} - \ket{\Phi(B(X);A)} \right\rVert^2\\
&= \min_X \left( \braket{\Phi(B(X);A) | \Phi(B(X);A) } - \braket{\Phi(B(X);A) | \chi } - \braket{\chi |\Phi(B(X);A) } \right) \nonumber.
\end{align}

Differentiating the complex conjugate of the entries of $X$, this minimum is characterized by

\begin{equation}
0 = \frac{\partial}{\partial \bar{X}} \left\lVert \ket{\chi} - \ket{\Phi(B(X);A)} \right\rVert^2 \, ,
\end{equation}

\begin{equation}
\diagram{appendix}{1}-\diagram{appendix}{2} = 0
\end{equation}

so that

\begin{equation}
\ket{\Phi(B(X);A)} = \sum_j \, \diagram{appendix}{3} \, .
\end{equation}

Using the completeness relation from eq:comptangent, we can rewrite this as

\begin{equation}
\ket{\Phi(B(X);A)} = \sum_j \quad \diagram{appendix}{4} - \diagram{appendix}{5} \, ,
\end{equation}

yielding the final result for the tangent space projector

\begin{equation}
P_{\ket{\Psi(A)}} = \sum_j \quad \diagram{appendix}{6} - \diagram{appendix}{7} \, .
\end{equation}

§.§ Vumps for fermionic MPOs

Consider a 1D fermionic Hamiltonian represented as an MPO,

\begin{equation}
H = T(O) = \diagram{appendix}{8} \, .
\end{equation}

Its ground state minimizes the real-valued Rayleigh quotient,

\begin{equation}
\frac{\braket{\Psi | H | \Psi}}{\braket{\Psi | \Psi}} \, ,
\label{eq:app:eigfunc}
\end{equation}

and so does its optimal MPS approximation $\ket{\Psi(\tilde{A})}$ when restricting $\ket{\Psi}$ to the MPS manifold. Put differently, the derivative of eq:app:eigfunc to the complex conjugate of the MPS parameters equals zero, ,

\begin{equation}
\label{eq:vumpsderiv}
\frac{\partial}{\partial \bar{A}} \left( \bra{\Psi(A)} \right) \left( H - \frac{\braket{\Psi(A) | H | \Psi(A)}}{\braket{\Psi(A) | \Psi(A)}} \right)  \ket{\Psi(A)} = 0
\end{equation}

when $A = \tilde{A}$. This expresses that the tangent space projection of the original eigenvalue equation equals zero in the optimal MPS, ,

\begin{equation}
P_{\ket{\Psi(\tilde{A})}} \left( H- E \right) \ket{\Psi(\tilde{A})} = 0 \, .
\end{equation}

Using the tangent space projector derived in app:tanspaceproj, this condition gives rise to the fixed-point equations, ${A^c}' \propto A^c$ and $C' \propto C$, where we dropped the tilde and

\begin{equation}
\diagram{appendix}{11} = \diagram{appendix}{12}\, \quad , \quad
\diagram{appendix}{13} = \diagram{appendix}{14} \, .
\end{equation}

The proportionality factor in the fixed-point equations is equal to the (infinite) eigenvalue $E$. Invoking the MPO structure of $H$ then leads to the set of eigenvalue and linear problems Eqs. (<ref>), (<ref>), (<ref>) and (<ref>), as detailed in sec:fvumps, that can be iteratively solved to arrive at the fixed-point solution (eq:vumpsfp).

Note that in our derivation of the VUMPS algorithm (eq:vumpsderiv) we have made use of the fact that $H$ is hermitian. However, all of the techniques involved, based on the determination of dominant eigenvalues (typically involving Krylov methods), also apply when searching the dominant eigenvector for a more general, non-hermitian MPO. However, the method is no longer variational then. In this regard, one could wonder if there exists a local way to express that a fermionic MPO is hermitian. , is there a condition on $O$ under which we know that $T$ is hermitian? Building on sec:tensorproperties, we know that hermiticity of such an operator comes down to

\begin{equation}
T = \diagram{appendix}{8} = \diagram{appendix}{33} = T^\dagger\, .
\end{equation}

Hence, one could infer a sufficient condition by relating the local tensors on both sides of the equality. However, as the virtual arrows point in opposite directions, we need to flip them. Furthermore, both MPOs might be related via a gauge transform. Hence, if there exists some arbitrary flipper satisfying

\begin{equation}
\diagram{appendix}{42} = \diagram{appendix}{43} \,
\end{equation}

such that the relation

\begin{equation}
\diagram{appendix}{34} = \diagram{appendix}{44} \, ,
\end{equation}

is satisfied, then $T(O)$ will be hermitian.

§.§ Vumps for a double-layer transfer matrix

This type of VUMPS algorithm is typically encountered when computing PEPS expectation values as the contraction of a 2D bra and ket layer with possibly a local operator in between. Hence, we will determine the leading eigenvector of the transfer operator

\begin{equation}
T(O) = \diagram{appendix}{62}
\end{equation}

where $O$ is a double-layer tensor as in doublelayertensor. For the fixed-point, we propose an MPS Ansatz,

\begin{equation}
\ket{\Psi(M)} = \diagram{appendix}{63} \, .
\end{equation}

The optimal $M$ tensors are then determined by projecting the eigenvalue equation for $T(O)$ onto the MPS tangent space, yielding

\begin{equation}
P_{\ket{\Psi(M)}} \left( T(O)- \Lambda \right) \ket{\Psi(M)} = 0 \, .
\end{equation}

To do so, we have to determine the tangent space projector $P_{\ket{\Psi(M)}}$ in and orthogonal to an MPS with two physical legs. Note therefore that all the concepts mentioned in app:tanspaceproj can readily be extended to multiple physical legs, with $P$ tensors on the upward physical legs as the only exception. In particular, we can bring the MPS in the mixed gauge

\begin{equation}
\ket{\Psi(M)} = \diagram{appendix}{64}
\end{equation}

\begin{equation}
\diagram{appendix}{65} = \diagram{appendix}{66} = \diagram{appendix}{67}
\end{equation}

\begin{equation}
\diagram{appendix}{68} = \diagram{appendix}{69} \qquad \text{and} \qquad \diagram{appendix}{70} = \diagram{appendix}{71} \, .
\end{equation}

Here, the circular tensors in the physical contractions again indicate a parity tensor. A tangent vector for the $\ket{\Psi(M)}$ states can always be expressed as

\begin{equation}
\ket{\Phi(B;M)} = \sum_j \diagram{appendix}{72} \, .
\end{equation}

When orthogonal to the MPS, we can parametrize the $B$ tensor as

\begin{equation}
\diagram{appendix}{73} \qquad \text{or} \qquad \diagram{appendix}{74} \, ,
\end{equation}

, the left, respectively, right canonical form where $V^l$ spans the right null-space of $\thiccbar M^l$ while $V^r$ describes the left null-space of $\thiccbar M^r$ so that

\begin{equation}
\diagram{appendix}{75} =0\, , \qquad \text{respectively,} \qquad \diagram{appendix}{76} =0
\end{equation}

and the $X$ tensor thus essentially contains all the parametric freedom in $\ket{\Phi(B;M)}$. Furthermore, both  $V^l$ and $V^r$ are normalized in the same way as their corresponding MPS tensor such that

\begin{equation}
\diagram{appendix}{77} + \diagram{appendix}{78}  = \diagram{appendix}{79} \quad \text{and} \quad
\diagram{appendix}{80} + \diagram{appendix}{81}  = \diagram{appendix}{82}
\end{equation}

due to completeness. In both gauges the overlap between two tangent vectors reduces to

\begin{equation}
\braket{\Phi(B_2;M)|\Phi(B_1;M)} = |\mathbb{Z}| \diagram{appendix}{83} = |\mathbb{Z}| \diagram{appendix}{84} \, .
\end{equation}

The tangent space projection in and orthogonal to $\ket{\Psi(M)}$ of a general translation-invariant state $\ket{\chi}$ is the tangent vector parametrized by $X$ for which

\begin{equation}
X = \min_X \left\lVert \ket{\chi} - \ket{\Phi(B(X);M)} \right\rVert^2 \, .
\end{equation}

Again differentiating with respect to the complex conjugate of the entries of $X$, this minimum is characterized by

\begin{equation}
\diagram{appendix}{85} - \diagram{appendix}{86} =0 \,.
\end{equation}

so that

\begin{equation}
\ket{\Phi(B(X);M)}=P_{\ket{\Psi(M)}} \ket{\chi} = \sum_j \diagram{appendix}{87} \, .
\end{equation}

Using the aforementioned completeness relations, we obtain

\begin{equation}
P_{\ket{\Psi(M)}} = \sum_j \diagram{appendix}{88} - \diagram{appendix}{89}\, .
\end{equation}

for the “two-leg” tangent space projector as the $P$ tensors appear due to the operator application. The MPS approximation for the fixed-point of $T(O)$ thus satisfies the typical VUMPS fixed-point equations, ${M^c}' \propto M^c$ and $C' \propto C$, where now

\begin{equation}
\diagram{appendix}{90} = \diagram{appendix}{91}
\quad \text{and} \quad
\diagram{appendix}{92} = \diagram{appendix}{93}\, .
\end{equation}

, $M^c$ and $C$ are the solutions to an eigenvalue problem which we can solve by determining the environments $G^l$ and $G^r$, defined by

\begin{equation}
\diagram{appendix}{94} = \lambda \diagram{appendix}{95} \qquad \text{and} \qquad \diagram{appendix}{96} = \lambda \diagram{appendix}{97} \, .
\label{eq:app:envvumpsmpo2}
\end{equation}

The VUMPS fixed-point equations then reduce to

\begin{equation}
\diagram{appendix}{90} = \diagram{appendix}{98} = \lambda \diagram{appendix}{99}
\end{equation}

\begin{equation}
\diagram{appendix}{92} = \diagram{appendix}{100} = \diagram{appendix}{101}\, .
\end{equation}

when the environments are normalized as

\begin{equation}
\diagram{appendix}{102} = 1\, .
\end{equation}

§ REFORMULATION OF GFTNS AS GENERIC FTNS

To define Gaussian fermionic tensor network states, the authors of Ref.<cit.> explore two approaches: the Kraus-Schuch (KS) method, which maximally exploits the Gaussian property and constructs GfTNS by applying a local Gaussian map to to a Gaussian state of maximally entangled pairs, and the Gu-Verstraete-Wen (GVW) construction, which primarily relies on the tensor network structure and utilizes Grassmann variables to contract local Gaussian tensors. Under straightforward constraints, these are equivalent and can be transformed into each other.

The GfTNS in this work (Sec. <ref>) are obtained by minimizing their energy density within the KS formalism. Subsequently, they are cast in their respective GVW form. As this definition resembles the $\ZZ_2$-graded formalism the most, we will use it as our starting point to transform GfTNS into generic fTNS. Therefore, consider an arbitrary GVW Gaussian fermionic tensor at site $\mathbf{n}$, defined as

\begin{equation}\label{eq:Tx}
\hat{T}_\mathbf{n} = \exp\left(\frac{1}{2}\chi^T_\mathbf{n} A \chi^\ph_\mathbf{n}\right) = \, \diagram{appendix}{103},
\end{equation}

where the column vector $\chi_\mathbf{n} \equiv \left(a^{\dagger}_{\mathbf{n}},\theta^{h}_{\mathbf{n}},\bar{\theta}^{h}_{\mathbf{n}},\theta^{v}_{\mathbf{n}},\bar{\theta}^{v}_{\mathbf{n}}\right)^T$ collects the $N$ physical creation operators $a^\dagger_{\mathbf{n},i}$ ($i=1,...,N$) as well as $4M$ virtual Grassmann variables $\theta^h_{\mathbf{n},\alpha}, \bar{\theta}^h_{\mathbf{n},\beta}, \theta^v_{\mathbf{n},\lambda}$ and $\bar{\theta}^v_{\mathbf{n},\sigma}$ ($\alpha, \beta, \lambda, \sigma =1,...,M$), all assigned to site $\mathbf{n}$. The latter are mutually anti-commuting and square to zero. Moreover, we make them anti-commute with the physical creation operators. The matrix $A \in \CC^{(N+4M)\times(N+4M)}$, hence anti-symmetric, contains the variational parameters determined by the energy minimization and is assumed to be independent of $\mathbf{n}$, i.e., we are restricting ourselves to translation invariant states. With this definition of Gaussian tensors in place, we can define an (unnormalized) contracted GfTNS via the Berezin integral,

\begin{equation}\label{contraction}
|\psi\rangle = \int [D\theta] \int [D\bar{\theta}] \prod_{\mathbf{n}} e^{\bar{\theta}_{\mathbf{n}}^{h^T} \theta_{\mathbf{n}+\mathbf{e}_x}^h } e^{\bar{\theta}_{\mathbf{n}}^{v^T} \theta_{\mathbf{n}+\mathbf{e}_y}^v } \hat{T}_{\mathbf{n}}|0\rangle\,,
\end{equation}

where $|0\rangle$ is the physical Fock vacuum and $\mathbf{e}_{x/y}$ are unit vectors along the $x/y$-direction.

To understand how the Grassmann variables can incorporate the correct fermionic features on the virtual level, we show that there exists a straightforward isomorphism between $\mathcal{G}(M)$, the Grassmann algebra of $M$ Grassmann variables and the Clifford algebra of $M$ fermionic creation and annihilation operators, represented on a $\ZZ_2$-graded Hilbert space of dimension $2^M$. Therefore, consider the case of a single Grassmann number $\theta$. To every monomial in $\theta$ we associate a basis state of the super vector space as follows,

\begin{equation}
\theta^n \cong |n\rangle\, ,\quad n\in \{0,1\} \, .
\label{eq:grassmannisomorphism}
\end{equation}

The dual space is isomorphic to polynomials of another Grassmann number $\bar{\theta}$,

\begin{equation}
\bar{\theta}^n \cong \langle n|\,,\quad n\in\{0,1\} \, .
\end{equation}

The evaluation/contraction map is then given by the following Berezin integral,

\begin{equation}
\mathcal{C}: \langle n|\otimes |m\rangle \cong \bar{\theta}^n\theta^m\rightarrow \int \mathrm{d}\theta \int\mathrm{d}\bar{\theta} \,e^{\bar{\theta}\theta} \bar{\theta}^n\theta^m = \langle n|m\rangle = \delta_{nm} \, .
\end{equation}

This mapping of monomials of Grassmann numbers to basis states of a super vector space generalizes straightforwardly to the case with more than one Grassmann number. Keeping this in mind, we can interpret the definition of the local tensor in eq:Tx as a coherent state that couples the $N$ physical fermions to $4M$ virtual fermions. Moreover, we can regard it as generic fermionic tensor of the form

\begin{equation}
\diagram{appendix}{104} \, .
\label{eq:gtensor}
\end{equation}

Indeed, $\theta$ Grassmann variables correspond to standard $\ZZ_2$-graded Hilbert spaces, while $\bar{\theta}$ variables yield a dual space, i.e., this distinction is equivalent to placing arrows on the tensor legs. The Berezin integral performs the contraction of the local tensors by means of the $e^{\bar{\theta}_\mathbf{n}\theta_{\mathbf{n}+\e_{x/y}}}$ factors. Note that the resulting composite bond dimension of the GfTNS is indeed $D= 2^M$ in each direction.

The relation between the virtual Grassmann variables in the GVW construction and $\ZZ_2$-graded Hilbert spaces allows us to rephrase the Gaussian states as generic fTNS and thus to determine their tensor entries in eq:gtensor based on the parametric $A$ matrix in eq:Tx. Indeed, also replacing the physical creation operators with Grassmann variables in $\chi$, we can expand the exponential in the resulting $T$ tensor as

\begin{equation}\label{expanded}
T = \sum_{\{n_i\}=0}^1 T_{n_1,n_2,\dots,n_N} \theta_1^{n_1}\theta_2^{n_2}\dots \theta_N^{n_N}\,,
\end{equation}

where for simplicity of notation, we have denoted all the different Grassmann variables as $\theta_i$, where the the index $i$ enumerates the Grassmann variables in the order they appear in $\chi$. Given a bit string of $n_i\,'s$, the corresponding tensor element can be obtained as follows: take the expansion in Eq. (<ref>), take all Grassmann variables for which $n_i=0$ and put them equal to zero by hand, and integrate over those Grassmann variables for which $n_i=1$. We now apply the same procedure to the exponential expression in eq:Tx. This gives

\begin{equation}
T_{\{n_i\}} = \prod_{i:n_i=1}\int \mathrm{d}\theta_i \exp\left(\frac{1}{2}\chi_{\{n_i\}}^T A[n_i] \chi_{\{n_i\}} \right)\,,
\end{equation}

where the product over integrals is descending order, and $\chi_{\{n_i\}}$ is the vector of Grassmann variables for which $n_i=1$, and $A[n_i]$ is the corresponding submatrix of $A$ where only the rows and columns corresponding to the Grassmann variables with $n_i=1$ are kept. The result of the integral is

\begin{equation}
T_{\{n_i\}} = \text{Pf}(A[n_i]) \, .
\end{equation}

The generic tensor entries can thus be obtained by calculating Pfaffians of submatrices of the parametric matrix $A$.