# Fermionic tensor networks for higher-order topological insulators from charge pumping

Anna Hackenbroich, B. Andrei Bernevig, Norbert Schuch, Nicolas Regnault

► To cite this version:

Anna Hackenbroich, B. Andrei Bernevig, Norbert Schuch, Nicolas Regnault. Fermionic tensor networks for higher-order topological insulators from charge pumping. Physical Review B, 2020, 101 (11), pp.115134. 10.1103/PhysRevB.101.115134. hal-02526733

HAL Id: hal-02526733

https://hal.science/hal-02526733v1

Submitted on 2 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fermionic tensor networks for higher-order topological insulators from charge pumping

Anna Hackenbroich<sup>1,2</sup> B. Andrei Bernevig<sup>3,4</sup> Norbert Schuch<sup>1,2</sup> and Nicolas Regnault<sup>3,5</sup>

<sup>1</sup>Max-Planck-Institute of Quantum Optics, Hans-Kopfermann-Str. 1, 85748 Garching, Germany

<sup>2</sup>Munich Center for Quantum Science and Technology, Schellingstraße 4, 80799 München, Germany

<sup>3</sup>Joseph Henry Laboratories and Department of Physics, Princeton University, Princeton, New Jersey 08544, USA

<sup>4</sup>Max-Planck-Institute of Microstructure Physics, 06120 Halle, Germany

<sup>5</sup>Laboratoire de Physique de l'École normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 24 December 2019; revised manuscript received 19 February 2020; accepted 19 February 2020; published 19 March 2020)

We apply the charge-pumping argument to fermionic tensor network representations of *d*-dimensional topological insulators (TIs) to obtain tensor network states (TNSs) for (*d* + 1)-dimensional TIs. We exemplify the method by constructing a two-dimensional projected entangled pair state (PEPS) for a Chern insulator starting from a matrix product state (MPS) in *d* = 1 describing pumping in the Su-Schrieffer-Heeger (SSH) model. In extending the argument to second-order TIs, we build a three-dimensional TNS for a chiral hinge TI from a PEPS in *d* = 2 for the obstructed atomic insulator (OAI) of the quadrupole model. The (*d* + 1)-dimensional TNSs obtained in this way have a constant bond dimension inherited from the *d*-dimensional TNSs in all but one spatial direction, making them candidates for numerical applications. From the *d*-dimensional models, we identify gapped next-nearest-neighbor Hamiltonians interpolating between the trivial and OAI phases of the fully dimerized SSH and quadrupole models, whose ground states are given by an MPS and a PEPS with a constant bond dimension equal to 2, respectively.

DOI: 10.1103/PhysRevB.101.115134

## I. INTRODUCTION

Higher-order TIs [1–4] have recently been introduced as a new class of symmetry-protected topological systems generalizing the framework of TIs with surface states [5]. According to one definition, a TI of order *n* in *d* spatial dimensions has topological boundary modes at the (*d* − *n*)-dimensional intersection of *n* crystal faces. In this terminology, conventional TIs such as Chern insulators [6] are of order *n* = 1 with protected boundary modes at (*d* − 1)-dimensional edges. On the other hand, second-order TIs protected by mirror or rotation symmetries have zero-dimensional corner states in *d* = 2 dimensions (in which case they may exemplify obstructed atomic limits [7] with local Wannier states) and one-dimensional chiral or helical hinge states in *d* = 3 dimensions (and bulk bands without a local Wannier description). Prototypical examples of second-order topological phases include the two-dimensional quadrupole model of Ref. [1], a natural extension [8] of the Su-Schrieffer-Heeger (SSH) model [9], and the three-dimensional chiral hinge insulator of Ref. [2], both of which have been experimentally observed in either materials [10] or mechanical [11], acoustic [12,13], photonic [14–17], and electrical [18–20] systems.

The bulk-boundary correspondence states that the topological properties of a system are reflected in the excitations at its physical boundary. For instance, a two-dimensional Chern insulator is characterized by an integer number of gapless chiral edge modes [21]. Similarly, second-order TIs in two dimensions possess gapless corner modes at the intersection of two edges compatible with the crystal symmetry [1]. In three dimensions, one type of second order TI is characterized by the presence or absence of a chiral hinge mode [2]. For strong TIs [22] as well as chiral topological phases, the universal features in the boundary energy spectrum are encoded in its entanglement spectrum (ES) [23] characterizing the bulk entanglement properties. The convenience of this bulk characterization makes the ES an important tool for the numerical analysis of topological phases. Recently, it was also observed, as expected, that the ES of higher-order TIs displays characteristic (*d* − *n*)-dimensional boundary states as long as the entanglement cut preserves the crystal symmetries of the phase [24–26].

A natural platform for the bulk-boundary correspondence is provided by tensor network states (TNSs) [27] in which the entanglement between physical particles is mediated through virtual particles hosting the lower-dimensional boundary theory [28]. Physical and virtual particles are related by a local tensor whose bond dimension determines the maximal amount of entanglement in the state [29,30]. The structure of the local tensor encodes the topological properties of the quantum state, and an analysis of this relation has led to valuable insight for the systematic understanding of one-dimensional symmetry-protected topological (SPT) phases [31] and intrinsically ordered topological phases [32]. Free-fermion systems

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.
can be described using Gaussian fermionic TNSs (GfTNSs) that are defined purely in terms of two-point correlation functions [33]. Starting from such Gaussian TNSs, interacting states can be constructed, whose topological properties derive from the features of the initial Gaussian model [34].

Many nonchiral topologically ordered phases in two dimensions [35,36] have a simple representation in terms of projected entangled pair states (PEPSs) [37] with a finite and constant bond dimension [38,39]. However, Gaussian PEPSs with chiral topological properties do not adequately describe the bulk of a gapped phase since their correlation functions decay algebraically [40,41]. Indeed, the existence of gapped Gaussian PEPSs with finite bond dimension for SPT phases of free fermions is forbidden in dimensions $d > 1$ [41,42]. The no-go theorem states that only Wannierizable phases corresponding to the product of one- and zero-dimensional systems can be exactly represented as a Gaussian TNS.

This no-go theorem does not prevent TNSs from being a useful numerical description at finite system size for chiral topological phases like the fractional quantum Hall effect [43]. Hence, numerically efficient TNSs for such phases are valuable, especially in two and three dimensions. One method for the construction of TIs in $d+1$ dimensions from $d$-dimensional TIs is given by charge pumping [44]. For instance, a two-dimensional Chern insulator is obtained from a charge-pumping interpolation in the one-dimensional SSH model when the periodic time direction is identified with the momentum in the second spatial direction [45]. For second-order TIs, dipole pumping in the quadrupole model defines a three-dimensional chiral hinge insulator [2,3]. In both cases, the zero-dimensional boundary modes of the $d$-dimensional model give rise to one-dimensional chiral boundary modes in the $(d+1)$-dimensional system.

In this paper, we exemplify how the charge-pumping argument applied to $d$-dimensional GfTNSs with constant bond dimension yields gapped GfTNSs for TIs in $d+1$ dimensions. By construction, the $(d+1)$-dimensional TNS has a constant finite bond dimension in a hybrid coordinate system with $d$ real-space axes and one momentum axis in the additional dimension. In order to obtain a real-space tensor network for the $(d+1)$-dimensional state, we apply to the hybrid TNS an inverse Fourier transform (FT) in the direction $d+1$. As a result, the bond dimension of the real-space local tensor in this direction generically grows with the system size due to the nonlocality of the FT, whereas it is identical to the finite bond dimension of the $d$-dimensional TNS in the other $d$ directions. We apply this construction both to a matrix product state (MPS) [46] for the SSH model in order to obtain a PEPS for a Chern insulator, and to a Gaussian PEPS with finite bond dimension for the topological quadrupole model in order to obtain a three-dimensional GfTNS for the second-order chiral hinge state TI of Refs. [2,3]. Therefore, our approach provides us with a gapped TNS with one-dimensional chiral boundary states and a constant finite bond dimension in all but one of the spatial directions. This representation is therefore potentially useful for tensor network algorithms.

This paper is structured as follows: In Sec. II we begin with a brief overview of charge pumping in the SSH model and fermionic TNSs, and continue by studying an MPS for the SSH ground state along the charge-pumping interpolation.

In the following Sec. III, we employ this MPS to construct a two-dimensional PEPS describing a Chern insulator. In Sec. IV we extend our method to the two-dimensional second-order quadrupole insulator and the construction of a three-dimensional PEPS for the chiral hinge state higher-order TI. Finally, we summarize our results and discuss remaining open questions in Sec. V.

## II. FERMIONIC MPS FOR CHARGE PUMPING IN THE SSH MODEL

In this section we introduce a fermionic MPS that describes a charge-pumping cycle in the SSH model corresponding to a Chern insulator with Chern number $C=1$ in two spatial dimensions. We begin by briefly reviewing the SSH model, the charge-pumping argument, and its relation to Chern insulators in Sec. II A. We continue with a short introduction to fermionic MPSs in Sec. II B. In Sec. II C, we construct the MPS for the SSH model and study its parent Hamiltonian.

### A. Chern insulator from charge pumping in the SSH model

The SSH model describes a one-dimensional chain of spinless fermions with two orbitals denoted A and B per unit cell [9]. We consider the model at half-filling where the number of particles is equal to the number $N_x$ of unit cells. For open boundary conditions, the Hamiltonian reads as

$$H_{\text{SSH}} = t^{(0)} \sum_{x=0}^{N_x-1} (\hat{a}_{A,x}^{\dagger} \hat{a}_{B,x} + \text{H.c.}) + t^{(1)} \sum_{x=1}^{N_x-1} (\hat{a}_{B,x-1}^{\dagger} \hat{a}_{A,x} + \text{H.c.}), \quad (1)$$

where we use the notation $\hat{a}_{A,x}$ and $\hat{a}_{B,x}$ for the fermionic annihilation operators of the orbitals A and B in unit cell $x$ with $x=0, \dots, N_x-1$, respectively. Here, $t^{(0)}$ denotes the hopping amplitude between A and B orbitals within the same unit cell, and $t^{(1)}$ the hopping amplitude between sites on neighboring unit cells [see Fig. 1(a)]. For $|t^{(1)}| > |t^{(0)}|$, the hopping between unit cells dominates the hopping within unit cells and the SSH model is in a phase topologically different from the case $|t^{(0)}| > |t^{(1)}|$. This phase is protected by spatial inversion symmetry and characterized by fermionic edge modes [45] and a “filling anomaly” [8]. It is called obstructed atomic insulator (OAI) [7]. When the intracell hopping $t^{(0)}$ vanishes in the deep OAI phase, the system is dimerized since its bulk splits into decoupled two-site blocks. In this case, the edge excitations are created by the mode operators $\hat{a}_{A,0}^{\dagger}$ and $\hat{a}_{B,N_x-1}^{\dagger}$ and have exactly zero energy in this specific model, but can generally be moved in energy. On the other hand, for $|t^{(0)}| > |t^{(1)}|$ the SSH model is trivial, with a dimerized point at $t^{(1)}=0$ and no state at zero energy.

A gapped interpolation between the trivial and OAI phases of the SSH model can be constructed by adding to the Hamiltonian $H_{\text{SSH}}$ the term

$$\sum_{x=0}^{N_x-1} (\mu_A \hat{a}_{A,x}^{\dagger} \hat{a}_{A,x} + \mu_B \hat{a}_{B,x}^{\dagger} \hat{a}_{B,x}) \quad (2)$$

introducing a staggered chemical potential $\mu_A = -\mu_B = \mu$ which breaks the inversion symmetry as shown in Fig. 1(a).
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Extended SSH model with a staggered chemical potential $\mu$, nearest-neighbor hoppings $t^{(0)}$ and $t^{(1)}$ within and between unit cells, and next-nearest-neighbor hoppings $t_A^{(2)}$ and $t_B^{(2)}$ on the A and B sublattices, respectively. The unit cells consisting of one site of each sublattice are marked by gray rectangles. (b) Action on the local MPS tensors of the U(1) symmetry ensuring that the MPS of Eq. (11) lies at half filling of the chain.</figcaption></figure>

In Ref. [3], the authors consider the periodic time-dependent Hamiltonian $H_{\text{pump}}(t)$ with time $t \in (-\pi, \pi]$ defined by the couplings

$$ (\mu, t^{(1)}, t^{(0)}) = \begin{cases} (\cos(t), 0, -\sin(t)), & -\pi < t \le 0 \\ (\cos(t), \sin(t), 0), & 0 < t \le \pi. \end{cases} \qquad (3) $$

At $t = -\pi$ (equivalent to $t = \pi$), the system is in an atomic phase with only the A orbitals occupied. For $-\pi < t < 0$, the coupling between A and B sites in the same unit cell is nonzero and the charge is transferred from left to right by the changing chemical potential until only the B orbitals are occupied at $t = 0$. At $t = -\pi/2$, the staggered chemical potential vanishes and the system is in the trivial dimerized phase of the SSH model. For $0 < t < \pi$, the intra-unit-cell coupling vanishes whereas the hopping between different unit cells is nonzero and the charge is transferred from left to right such that at $t = \pi$ the system returns to the state with all A orbitals occupied. At $t = \pi/2$, the chemical potential vanishes and the system is in the OAI dimerized phase of the SSH model. The ground state of the pumping Hamiltonian $H_{\text{pump}}(t)$ is continuous as a function of time $t$ if the chain has periodic and antiperiodic boundary conditions for $N_x$ odd and even, respectively.

The charge-pumping interpolation of Eq. (3) can be used to generate a lattice model in one dimension higher with the topology of a Chern insulator [3]. Indeed, $H_{\text{pump}}(t)$ corresponds to the time-dependent Bloch Hamiltonian

$$ \mathcal{H}_{\text{pump}}(k_x, t) = [t^{(0)}(t) + t^{(1)}(t) \cos k_x] \sigma_1 + t^{(1)}(t) \sin k_x \times \sigma_2 + \mu(t) \sigma_3, \quad (4) $$

where $\sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$, $\sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}$, and $\sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$ are the Pauli matrices and $k_x = \frac{\pi}{N_x}(2j - N_x) \in [-\pi, \pi]$ for $0 \le j \le N_x - 1$ the lattice momentum. Since the interpolation is cyclic, the time $t \in [-\pi, \pi]$ may be interpreted as the lattice

momentum $k_y$ of a second spatial direction y and Eq. (4) as the Bloch Hamiltonian of a two-dimensional system with closed boundaries in both directions. This Hamiltonian has a Chern number $C = 1$ due to the charge transport between unit cells induced by the interpolation for $t \in [0, \pi] [3]$.

We note that the cycle of Eq. (3) can be deformed to a smooth charge-pumping interpolation, without changing the topology or breaking the dimerization. In this cycle, the couplings $\mu(t)$, $t^{(1)}(t)$, and $t^{(0)}(t)$ and hence the Bloch Hamiltonian are smooth functions of the time $t$, i.e., they are infinitely often continuously differentiable with respect to $t$. As explained in the previous paragraph, a two-dimensional Chern insulator model is obtained by identifying the time $t$ and the momentum $k_y$. Due to the smoothness of the cycle, the real-space representation of the Chern insulator has couplings which decay faster than any polynomial [47]. The smooth cycle is obtained by smoothing out the nonanalyticities of Eq. (3) at $t = 0, \pi$ using smooth functions that interpolate between 0 and 1. This can be done in such a way that, like in Eq. (3), at each time $t$ the smooth couplings satisfy $\mu(t)^2 + t^{(1)}(t)^2 + t^{(0)}(t)^2 = 1$ [48].

## B. Fermionic MPSs for a half-filled lattice

Similarly to bosonic MPSs for spin chains, fermionic MPSs describing chains of electrons are obtained by associating virtual particles to each physical particle which mediate the entanglement between different physical constituents [33]. In the case of fermionic tensor networks, the virtual particles obey fermionic statistics. A fermionic MPS with $f$ physical complex fermionic modes per site and $\xi$ virtual complex fermionic modes per lattice site and nearest-neighbor bond has physical dimension $d_p = 2^f$ and bond dimension $D = 2^\xi$. The state is fully characterized by the set of local maps which relate the virtual and physical particles associated to one lattice site. With respect to orthonormal bases $\{|i\rangle\}$ with $i = 0, \dots, d_p - 1$ for the physical Hilbert space and $\{|l\rangle\}$, $\{|r\rangle\}$ with $l, r = 0, \dots, D - 1$ for the left and right virtual spaces on one site, respectively, the local maps are represented by $d_p$ local matrices $A_{lr}^i$ of dimension $D \times D$ [49].

In order to fix the global parity of the state, we restrict ourselves to parity-even local tensors which preserve the fermion parity between the physical and virtual layers [49]. We choose the orthonormal basis for the physical Hilbert space such that all basis states $|i\rangle$ have either even fermion parity $|i| = 0$ or odd fermion parity $|i| = 1$, and similarly for the virtual Hilbert spaces. A local tensor is parity even if for all nonvanishing elements $A_{lr}^i \neq 0$, the total fermion parity $|i| + |l| + |r| = 0$ is even, such that

$$ (-1)^{|i|} A_{lr}^i = (-1)^{|l|} A_{lr}^i (-1)^{|r|} \qquad (5) $$

for all configurations of $i, l$, and $r$ (where no summation is implied). If the MPS is constructed from tensors $A[j]$ at each site $j$ with this property, the state on a closed chain with $N_x$ sites is given by

$$ |\psi\rangle = \sum_{i_0, \dots, i_{N_x-1}} \sum_l (-1)^{\epsilon^{|l|}} \times (A[0]^{i_0} \dots A[N_x-1]^{i_{N_x-1}})_{ll} |i_0, \dots, i_{N_x-1}\rangle. \quad (6) $$
Here, $|l|$ is the parity of the virtual basis state on the link between the first and last site, and $\epsilon = 1$ or 0 corresponding to periodic and antiperiodic boundary conditions for the many-body state, respectively [49].

Below, we construct a fermionic MPS with parity-even local tensors for the ground state of the SSH model along the dimerized charge-pumping interpolation. Due to the parity symmetry of Eq. (5), this MPS necessarily has an even number of physical particles on a closed chain where all virtual bonds are contracted. On the other hand, the ground state of the SSH model is half-filled such that the number of particles is odd if the number of unit cells is odd. Therefore, to construct the parity-even SSH MPS we use a many-body basis built by acting with creation operators on the state $|\Omega\rangle$ that contains $N_x$ physical particles. To that end, we define new mode operators $a_{A,x}$ and $a_{B,x}$ for the physical fermions by performing a particle-hole transformation on all $B$ orbitals of the SSH chain while leaving the $A$ orbitals unaltered,

$$a_{A,x} \equiv \hat{a}_{A,x}, \qquad (7a)$$

$$a_{B,x} \equiv \hat{a}_{B,x}^{\dagger}. \qquad (7b)$$

The state $|\Omega\rangle$ is the vacuum of the new operators, given in terms of the original vacuum state $|\hat{\Omega}\rangle$ with $\hat{a}_{A,x}|\hat{\Omega}\rangle = \hat{a}_{B,x}|\hat{\Omega}\rangle = 0$ as

$$|\Omega\rangle = \prod_{x} \hat{a}_{B,x}^{\dagger} |\hat{\Omega}\rangle. \qquad (8)$$

It satisfies $a_{A,x}|\Omega\rangle = a_{B,x}|\Omega\rangle = 0$. Therefore, the new vacuum contains $N_x$ of the original fermions and thus is half-filled. The MPS of the SSH model along the interpolation is then defined with respect to the Fock states constructed from the modes $a_{A,x}^{\dagger}$ and $a_{B,x}^{\dagger}$ acting on the new vacuum $|\Omega\rangle$. This particle-hole transformation is motivated by our desire to write the MPS in terms of parity-even local tensors, which will in turn enable us to express the state as a GfTNS and thereby compute a Bloch parent Hamiltonian analytically.

### C. MPS for the SSH model

In this section, we study a fermionic MPS for the SSH model. In Sec. II C 1, we define the state in terms of its local tensors and identify the U(1) symmetry leading to a conserved number of particles. In Sec. II C 2, we derive a parent Hamiltonian for the MPS, allowing us to conclude in Sec. II C 3 that the MPS describes the SSH model along the charge-pumping cycle of Eq. (3).

#### 1. Local tensors and U(1) symmetry

In order to describe the ground state of the SSH model along the dimerized charge-pumping interpolation, we consider a fermionic MPS with physical dimension $d_p = 2$ corresponding to one fermionic mode per site and bond dimension $D_{SSH} = 2$. This is the minimal bond dimension for an MPS along the charge-pumping cycle since the ES of an open SSH chain in the OAI and trivial dimerized phases has two degenerate levels with respect to all cuts between and within unit cells, respectively.

The fermionic MPS is translation invariant with a unit cell of two sites. It is therefore fully specified by the local

matrices $A_{Ir}^i$ and $B_{Ir}^i$ with $i, l, r \in \{0, 1\}$ for sites on the $A$ and $B$ sublattices, respectively. In terms of the local tensors the physical state on a closed chain is given as

$$|\psi\rangle = \sum_{\substack{i_{A,0}, \dots, i_{A,N_x-1} \\ i_{B,0}, \dots, i_{B,N_x-1}}} |i_{A,0}, i_{B,0}, i_{A,1}, \dots, i_{B,N_x-1}\rangle \times \sum_l (-1)^{\epsilon^{|l|}} \times (A^{i_{A,0}} B^{i_{B,0}} A^{i_{A,1}} \dots B^{i_{B,N_x-1}})_{ll}, \qquad (9)$$

where $\epsilon = 1$ or 0 for periodic and antiperiodic boundary conditions. The physical many-body basis state is

$$|i_{A,0}, i_{B,0}, i_{A,1}, \dots, i_{B,N_x-1}\rangle = (a_{A,0}^{\dagger})^{i_{A,0}} (a_{B,0}^{\dagger})^{i_{B,0}} (a_{A,1}^{\dagger})^{i_{A,1}} \dots (a_{B,N_x-1}^{\dagger})^{i_{B,N_x-1}} |\Omega\rangle \qquad (10)$$

with the vacuum $|\Omega\rangle$ from Eq. (8). Guided by the dimerized limits to be discussed below, for the local MPS matrices we make the ansatz

$$A^0 = \begin{pmatrix} \gamma & 0 \\ 0 & 0 \end{pmatrix}, \quad A^1 = \begin{pmatrix} 0 & \beta \\ -\alpha & 0 \end{pmatrix}, \qquad (11a)$$

$$B^0 = \begin{pmatrix} \gamma & 0 \\ 0 & 0 \end{pmatrix}, \quad B^1 = \begin{pmatrix} 0 & -\alpha \\ -\beta & 0 \end{pmatrix}, \qquad (11b)$$

depending on parameters $\alpha, \beta, \gamma \in \mathbb{R}$. Note that the normalized quantum state defined by these local MPS matrices depends only on the two quotients $\alpha/\gamma$ and $\beta/\gamma$. Nonetheless, we choose to work with the parametrization of Eq. (11) since the case $\gamma = 0$ can be treated more conveniently without divergences.

In order to motivate the ansatz of Eq. (11) for the local matrices, let us see how the MPS from Eq. (9) can describe the ground state of the SSH model both in the trivial and the OAI dimerized phases with appropriate choices for the parameters $\alpha, \beta$, and $\gamma$. In the trivial dimerized phase of the SSH model, each unit cell decouples from the rest of the system and is in the state $(a_A^\dagger a_B^\dagger - 1)|\Omega\rangle$. In order to ensure the absence of entanglement on the bonds between unit cells, the blocked MPS matrices $(A^{i_A} B^{i_B})_{Ir}$ for one unit cell should be nonzero only if $l = r = 0$. Then, the restriction of the MPS to the unit cell, $\sum_{i_A, i_B} (A^{i_A} B^{i_B})_{00} |i_A, i_B\rangle$, is proportional to $(a_A^\dagger a_B^\dagger - 1)|\Omega\rangle$ if the blocked MPS matrices have only two nonzero entries $(A^0 B^0)_{00} = -(A^1 B^1)_{00}$. Therefore, the MPS represents the trivial dimerized phase of the SSH model if the parameters are chosen as $\alpha = 0$ and $|\beta| = |\gamma| \neq 0$. Similarly, in the OAI dimerized phase the chain splits into decoupled plaquettes composed of two adjacent sites from different unit cells. The state on each of the plaquettes is given by the superposition $(a_A^\dagger a_B^\dagger + 1)|\Omega\rangle$ such that the blocked MPS matrices for the plaquette should have only two nonzero entries $(B^0 A^0)_{00} = (B^1 A^1)_{00}$. This is achieved if $\beta = 0$ and $|\alpha| = |\gamma| \neq 0$. The ansatz of Eq. (11) is chosen such as to allow an interpolation between these two cases.

The MPS of Eq. (9) with the parametrization of Eq. (11) has a U(1) symmetry which ensures that the physical state lies at half-filling of the chain. The U(1) rotation acts on a single complex fermion with the matrix $U(\varphi) = \begin{pmatrix} 1 & 0 \\ 0 & e^{i\varphi} \end{pmatrix}$. The local tensors for the two sublattices are invariant under a combination of U(1) rotations of the physical and virtual legs
as sketched in Figs. 1(b) and 1(c) (see for instance Ref. [50] for MPSs with physical symmetries):

$$A_{lr}^i = \sum_j \sum_{l'r'} U(\varphi)_{ij} U(\varphi)_{l'r'}^\dagger A_{l'r'}^j U(\varphi)_{r'r}^\dagger, \quad (12a)$$

$$B_{lr}^i = \sum_j \sum_{l'r'} U(\varphi)_{ij}^\dagger U(\varphi)_{ll'} B_{l'r'}^j U(\varphi)_{r'r}. \quad (12b)$$

If we choose $\varphi = \pi$, these identities correspond to the parity symmetry of Eq. (5), showing that the local MPS tensors are parity even. For each nearest-neighbor bond, one virtual leg transforms with $U(\varphi)$ and the other with $U(\varphi)^\dagger$ in Eq. (12), such that the two U(1) rotations cancel if the bond is contracted. Therefore, the state on a chain with closed boundaries after the contraction of all virtual bonds is invariant under the physical part of the U(1) rotations of Eq. (12), given by staggered transformations $U(\varphi)$ and $U(\varphi)^\dagger$ on A and B orbitals, respectively. Invariance under this global symmetry implies that

$$0 = \left\langle \sum_x [a_{A,x}^\dagger a_{A,x} - a_{B,x}^\dagger a_{B,x}] \right\rangle \\ = \left\langle \sum_x [\hat{a}_{A,x}^\dagger \hat{a}_{A,x} + \hat{a}_{B,x}^\dagger \hat{a}_{B,x}] \right\rangle - N_x, \quad (13)$$

forcing the number of particles measured in terms of the original physical modes to be equal to the number of unit cells as required for the SSH ground state. In Eq. (12), the physical legs on the A and B sublattices transform as particles and holes, respectively, as expected due to the particle-hole transformation of Eq. (7). In addition, Eq. (12) implies that the virtual legs on the A and B sublattices transform as holelike and particlelike degrees of freedom (DOFs), respectively.

In order to gain a better understanding of the parameters $\alpha$, $\beta$, and $\gamma$ as well as of the systems described by the MPS from Eq. (9), it is helpful to consider a parent Hamiltonian for which it is the exact ground state. This Hamiltonian can be computed directly in terms of its Bloch representation once the MPS of Eq. (9) is expressed as a Gaussian fermionic TNS.

## 2. Parent Hamiltonian

Since the charge-pumping Hamiltonian of Eq. (3) describes noninteracting fermions, its ground state is a fermionic Gaussian state. It is thus fully characterized by its covariance matrix (CM) (see Appendix A for a summary of our conventions) [51]. As we review in Appendix B, the tensor network formalism may be used to construct Gaussian fermionic TNSs which are ground states of free-fermion Hamiltonians [33]. The CM in Fourier space of a translationally invariant Gaussian TNS is given by a simple expression which can often be evaluated analytically. Then, any positive function $\epsilon(\mathbf{k}) > 0$ on the Brillouin zone gives rise to a parent Hamiltonian with dispersion relation $\epsilon(\mathbf{k})$, whose Bloch representation is obtained by multiplying the CM by $\epsilon(\mathbf{k})$.

For all values of the parameters $\alpha$, $\beta$, and $\gamma$, the MPS of Eq. (11) corresponds to a Gaussian fermionic TNS whose Fourier CM is computed analytically in Appendix E. We show that this MPS with bond dimension $D=2$ is the

unique ground state of a longer-range SSH-like model with a staggered chemical potential $\mu$ and next-nearest-neighbor hoppings $t_A^{(2)}$ and $t_B^{(2)}$ on the A and B sublattices, respectively, as sketched in Fig. 1(a). The coupling constants of the parent Hamiltonian $H_{MPS}$ depend on the parameters of the MPS as

$$\mu = \frac{\gamma^4 - \alpha^4 - \beta^4}{\alpha^4 + \beta^4 + \gamma^4}, \quad (14a)$$

$$t^{(0)} = \frac{2\beta^2\gamma^2}{\alpha^4 + \beta^4 + \gamma^4}, \quad (14b)$$

$$t^{(1)} = \frac{2\alpha^2\gamma^2}{\alpha^4 + \beta^4 + \gamma^4}, \quad (14c)$$

$$t_A^{(2)} = -t_B^{(2)} = \frac{-2\alpha^2\beta^2}{\alpha^4 + \beta^4 + \gamma^4}, \quad (14d)$$

where $t^{(0)}$ and $t^{(1)}$ denote the hopping amplitudes between nearest-neighbor sites on the same and adjacent unit cells, respectively.

Depending on the parameter values, $H_{MPS}$ describes different phases of the fermionic chain. If any two out of the three parameters $\alpha$, $\beta$, $\gamma$ vanish, the system is in an atomic state without entanglement between different sites. Indeed, if only $\gamma$ is nonzero, all hopping constants vanish and the chemical potential is $\mu = +1$ such that we obtain the state with all B sites occupied. On the other hand, if either only $\alpha$ or only $\beta$ is nonzero, all hopping constants vanish and the chemical potential is $\mu = -1$ such that the MPS is equal to the state with all A sites occupied. This implies that the atomic state with occupied A orbitals is obtained from two distinct parameter configurations $(\alpha, \beta, \gamma) = (1, 0, 0)$ and $(\alpha, \beta, \gamma) = (0, 1, 0)$ whose corresponding local MPS matrices are related by a virtual unitary gauge transformation with representation matrix $\sigma_1$.

On the other hand, if only $\gamma$ and $\beta$ are nonzero, the system is in a dimerized phase where the next-nearest-neighbor hopping as well as the nearest-neighbor hopping between unit cells vanishes, whereas the nearest-neighbor hopping within unit cells is finite. Unless $|\gamma| = |\beta|$, inversion symmetry is broken by the nonzero chemical potential. For $|\gamma| = |\beta|$, we obtain the trivial dimerized phase of the SSH model. Similarly, if only $\gamma$ and $\alpha$ are nonzero, all couplings vanish except for the nearest-neighbor hopping between unit cells and the chemical potential. For $|\gamma| = |\alpha|$, inversion symmetry is restored and we obtain the OAI dimerized phase of the SSH model. Finally, if all three parameters are nonzero, both nearest-neighbor hopping constants $t^{(0)}$ and $t^{(1)}$ are nonzero and there is an additional next-nearest-neighbor hopping $t_A^{(2)} = -t_B^{(2)}$ which is odd under inversion.

## 3. Charge-pumping interpolation

We now use the MPS of Eq. (11) to describe charge pumping by considering the evolving state obtained from time-dependent parameters $\alpha(t)$, $\beta(t)$, and $\gamma(t)$ with time $t \in (-\pi, \pi]$. In Eq. (14) we saw that the MPS can represent charge pumping with a dimerized nearest-neighbor as well as a longer-range nondimerized Hamiltonian. For simplicity, we focus on the dimerized case. For the MPS given by the
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a), (b) Coupling constants of $H_{MPS}$ and (c), (d) single-particle ES for the MPS of Eq. (11) along topologically nontrivial and trivial interpolations. In (a) and (c), along the charge-pumping interpolation of Eq. (15) corresponding to a Chern insulator with Chern number $C = 1$. In (b) and (d), along the topologically trivial interpolation of Eq. (17). In both (a) and (b), the vanishing next-nearest-neighbor coupling $t_A^{(2)} = t_B^{(2)} = 0$ of $H_{MPS}$ is not shown. The ES in (c) and (d) was computed on half of a periodic chain with $N_x = 10$ unit cells.</figcaption></figure>

parametrization

$$ \phi_{\text{pump}} : (\alpha(t), \beta(t), \gamma(t)) = \begin{cases} (0, \sqrt{|\sin t/2|}, \sqrt{|\cos t/2|}), & -\pi < t \le 0 \\ (\sqrt{|\sin t/2|}, 0, \sqrt{|\cos t/2|}), & 0 < t \le \pi \end{cases} \quad (15) $$

the parent Hamiltonian $H_{MPS}$ is exactly equal to the time-dependent Hamiltonian $H_{pump}(t)$ of the dimerized charge-pumping cycle of Eq. (3). The evolution of the couplings in $H_{MPS}$ along the parametrization is shown in Fig. 2(a). Hence, the MPS of Eq. (11) can describe the ground state of the pumping model of Eq. (3) with the topology of a Chern insulator for all times $t \in (-\pi, \pi]$.

The ground state of the charge-pumping cycle $H_{pump}(t)$ and its correlation functions are continuous along the entire interpolation. However, the local MPS tensors parametrized by $\phi_{pump}$ are discontinuous at the point $t = \pm \pi$ along the cycle: $\alpha \to 0$ and $\beta \to 1$ as $t \to -\pi$ from above, whereas $\alpha \to 1$ and $\beta \to 0$ as $t \to \pi$ from below. In fact, at $t = \pm \pi$ the system is in the atomic state with all A orbitals occupied. As discussed in the paragraph beneath Eq. (14), there are two distinct configurations for the MPS parameters corresponding to this state which are attained for $t = \pm \pi$.

For the MPS we are considering, the discontinuity in the parametrization $\phi_{pump}$ is related to the chiral edge mode of the Chern insulator defined by the interpolation. Indeed, charge pumping in the SSH model with open boundaries corresponds to a Chern insulator on a cylinder with periodic and open boundaries in the vertical and horizontal directions, respectively. In the topological phase with Chern number $C = 1$, each edge of the cylinder hosts a one-dimensional chiral mode [21]. These modes are reflected in the ES of the SSH chain

along the charge-pumping cycle. For instance, the single-particle ES [52] of half the chain with periodic boundaries computed from the representation of the state as a Gaussian fermionic TNS along $\phi_{pump}$ is shown in Fig. 2(c). It has both a left- and a right-moving mode which are localized at the two virtual edges and which are degenerate in the SSH topological phase at $t = \pi/2$.

Using the representation of the MPS in terms of local tensors, we can also compute the many-body ES of half of an infinite chain with open boundaries. Indeed, the many-body ES is isometric to the spectrum of the logarithm of the normalized and symmetrized left and right fixed points of the MPS transfer matrix [28]. In the case of the MPS from Eq. (11), the fixed point is a matrix of dimension $2 \times 2$ with eigenvalues

$$ \frac{1}{2} \pm \frac{-\alpha^4 + \beta^4 + \gamma^4}{2\sqrt{(-\alpha^4 + \beta^4 + \gamma^4)^2 + 4\alpha^2\beta^2\gamma^4}} \quad (16) $$

Therefore, the many-body ES has two nontrivial levels that are related by normalization and which describe the spinless fermion at the single virtual boundary of the half-infinite chain. If and only if $\alpha^4 = \beta^4 + \gamma^4$, the two eigenvalues from Eq. (16) are degenerate. This corresponds to a crossing of the Fermi level by the edge fermion. If the interpolation describes a Chern insulator with Chern number $C = 1$, a degeneracy of the eigenvalues from Eq. (16) should therefore occur exactly once along the cycle. For the parametrization $\phi_{pump}$ this happens in the SSH topological phase at $t = \pi/2$ with $\alpha = \gamma$ and $\beta = 0$. Indeed, $\alpha^4 < \beta^4 + \gamma^4$ for $t < \pi/2$ and $\alpha^4 > \beta^4 + \gamma^4$ for $t > \pi/2$. Hence, the nontrivial Chern number implies that for MPS of the form of Eq. (11) the parametrization has to be discontinuous along the cycle in order to combine the two parts of the interpolation where $\alpha^4 \gtrsim \beta^4 + \gamma^4$ with only a single degenerate point.

We emphasize that this discontinuity in the MPS interpolation is required by the topology of the charge-pumping cycle, even if the pumping cycle itself is infinitely often continuously differentiable. For example, the smooth deformation of the topological pumping cycle, whose existence is discussed in Sec. II A, has an MPS ground state of the form of Eq. (11). Due to the nontrivial topology of the cycle, this MPS has a discontinuous interpolation despite the smoothness of its Bloch Hamiltonian.

The MPS of Eq. (11) can also describe cyclic interpolations corresponding to topologically trivial two-dimensional models. For instance, the MPS with parametrization

$$ \phi_{\text{triv}} : (\alpha(t), \beta(t), \gamma(t)) = \left( \frac{1}{2} e^{-(\tan \frac{|t|}{2})^{-2}}, 0, 1 \right) \quad (17) $$

for $t \in (-\pi, \pi]$ corresponds to the ground state of a model with finite chemical potential and nearest-neighbor hopping $t^{(1)}$ between unit cells, whereas the hopping within unit cells and the next-nearest-neighbor hopping vanish [see Fig. 2(b)]. From the single-particle ES in Fig. 2(d) we see that the conduction and valence bands in this system are not connected by the edge modes such that the Chern number is zero. Indeed, both the MPS parametrization $\phi_{triv}$ and the
corresponding Bloch parent Hamiltonian are smooth along the cycle.

### III. CHERN INSULATOR PEPS FROM SSH CHARGE PUMPING

In this section we discuss how charge pumping can be used to construct tensor networks in $d+1$ dimensions starting from $d$-dimensional TNSs. Specifically, we show how the MPS from Sec. II for the SSH model along a charge-pumping cycle leads to a two-dimensional PEPS for a Chern insulator. In Sec. III B we define a hybrid real-momentum space PEPS with finite bond dimension for the Chern insulator. In Sec. III C, we perform an inverse FT in the direction $d+1$ required to transform the state to a fully real-space representation. In Sec. III D we study how the resulting state can be expressed as a real-space PEPS, whose bond dimension is investigated in Sec. III E.

#### A. ($d+1$)-dimensional TIs from charge pumping

Charge pumping provides a systematic method to obtain a TI in $d+1$ spatial dimensions from a $d$-dimensional TI. Indeed, if the Bloch Hamiltonian of the $d$-dimensional model is smooth along the charge-pumping interpolation as a function of the cyclic time $t \in (-\pi, \pi]$, $t$ can be identified with the momentum $k_{d+1}$ in the ($d+1$)st direction of the ($d+1$)-dimensional system. The time-dependent Hamiltonian of the $d$-dimensional model then gives the Bloch Hamiltonian of the ($d+1$)-dimensional system as a function of $k_{d+1}$. For instance, charge pumping in the SSH model with $d=1$ defines a two-dimensional Chern insulator with Chern number $C=1$ (cf. Sec. II A).

We discretize the ($d+1$)st dimension with a finite number $N_{d+1}$ of lattice sites. The discrete momentum values in $(-\pi, \pi]$ are $k_{d+1}^{(j)} = \frac{\pi}{N_{d+1}}(2j - N_{d+1})$ for even $N_{d+1}$, and $k_{d+1}^{(j)} = \frac{\pi}{N_{d+1}}(2j - N_{d+1} + 1)$ for odd $N_{d+1}$. They are identified with discrete times $t^{(j)} = k_{d+1}^{(j)}$ for $j=0, \dots, N_{d+1}-1$.

Let us express the identification of time $t$ and momentum $k_{d+1}$ as a relation between the mode operators of the hybrid and real-space systems. The $d$-dimensional system has annihilation operators $\hat{a}_{\tau,x}(t^{(j)})$ for the physical fermion on the orbital $\tau=1, \dots, f$ on the unit cell $\mathbf{x} \in \mathbb{Z}^d$, where $f$ is the number of orbitals per unit cell. They depend on the discretized time value $t^{(j)}$ along the pumping interpolation at which the $d$-dimensional model is evaluated. For example, in the SSH model we have mode operators $\hat{a}_{A,x}(t^{(j)})$, $\hat{a}_{B,x}(t^{(j)})$ with $x=0, \dots, N_x-1$.

The ($d+1$)-dimensional model obtained from charge pumping has the same number of orbitals as the $d$-dimensional model it is derived from. For instance, the Chern insulator constructed from the SSH model also has sublattices A and B. Therefore, the physical creation operators of the ($d+1$)-dimensional system in real space are $\hat{a}_{\tau,(x,x_{d+1})}$, where $x_{d+1}=0, \dots, N_{d+1}-1$ is the real-space coordinate in the ($d+1$)st direction and $\tau=1, \dots, f$. The charge-pumping construction of the ($d+1$)-dimensional model requires periodic boundary conditions in the direction $d+1$. We may therefore consider the FT of the mode operators in the

direction $d+1$, while keeping the real-space coordinate $\mathbf{x}$ in the other $d$ dimensions, given by

$$ \hat{a}_{\tau, (\mathbf{x}, k_{d+1}^{(j)})} = \sum_{x_{d+1}=0}^{N_{d+1}-1} \mathcal{F}_{k_{d+1}^{(j)} x_{d+1}}^{\tau, (\mathbf{x}, x_{d+1})} \quad (18a) $$

with

$$ \mathcal{F}_{k_{d+1}^{(j)} x_{d+1}} = \frac{e^{-ik_{d+1}^{(j)} x_{d+1}}}{\sqrt{N_{d+1}}}. \quad (18b) $$

In terms of these mode operators, the identification of time $t$ and momentum $k_{d+1}$ is the expressed by the identity

$$ \hat{a}_{\tau, \mathbf{x}}(t^{(j)}) = \hat{a}_{\tau, (\mathbf{x}, k_{d+1}^{(j)})} \quad (19) $$

for all $d$-dimensional unit cells $\mathbf{x}$ and sublattices $\tau=1, \dots, f$. For example, in the case of the SSH charge pumping we have $\hat{a}_{\tau,x}(t^{(j)}) = \hat{a}_{\tau,(x,k_{d+1}^{(j)})}$ for $\tau=A, B$ and $x=0, \dots, N_x-1$.

The Bloch Hamiltonian of the ($d+1$)-dimensional model is the time-dependent Hamiltonian of the $d$-dimensional system. Hence, with the identification of Eq. (19), the ($d+1$)-dimensional ground state is given by the direct product

$$ |\psi_{d+1}\rangle = \bigotimes_{j=0}^{N_{d+1}-1} |\psi_d(t^{(j)})\rangle \quad (20) $$

of the many-body ground states $|\psi_d(t)\rangle$ of the $d$-dimensional model evaluated at the $N_{d+1}$ discrete times along the interpolation. From Eq. (19) it is clear that this defines the ground state with respect to hybrid ($d+1$)-dimensional real-momentum space coordinates $(\mathbf{x}, k_{d+1}^{(j)})$.

In order to obtain the state in terms of ($d+1$)-dimensional real-space coordinates $(\mathbf{x}, x_{d+1})$, we need to apply the inverse of the FT of Eq. (18) to the hybrid state of Eq. (20). If the $d$-dimensional pumping Bloch Hamiltonian is smooth as a function of time, the real-space correlation functions are guaranteed to decay faster than any polynomial [47].

#### B. Hybrid real-momentum-space Chern PEPS

We now specialize this construction, which we explained above in terms of generic free-fermionic TIs, to $d$-dimensional TIs which are described by a GfTNS at all times along their charge-pumping cycle. We will thereby obtain a GfTNS for the ($d+1$)-dimensional TI.

For the remainder of this section and for pedagogical purposes, we will focus on the one-dimensional case (i.e., $d=1$) and the corresponding notation. We will mostly rely on charge pumping in the SSH model which is described by the MPS with bond dimension $D_{SSH}=2$ from Eq. (9). Hence, we will obtain a Gaussian fermionic PEPS for the Chern insulator with Chern number $C=1$. Extensions to other models and higher dimensions are straightforward.

In case the $d$-dimensional ground state along the pumping cycle is given as a GfTNS, Eq. (20) allows to define the ($d+1$)-dimensional ground state as a hybrid TNS, where the first $d$ dimensions correspond to real space and the ($d+1$)st dimension is expressed in momentum space. Indeed, the local tensor of the hybrid TNS at the position $(\mathbf{x}, k_{d+1}^{(j)})$ is given
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Hybrid real-momentum-space PEPS for $N_y = 3$ sites in the vertical direction obtained by stacking rows of the MPS of Eq. (11) evaluated at different times $t^{(j)}$ along the charge-pumping interpolation. Physical and virtual legs transforming as particles and holes under the U(1) symmetry of Eq. (12) are marked in red and blue, respectively. Since the vertical bond dimension is trivial, i.e., equal to one, it is omitted in this sketch. (b) Inverse FT $\tilde{\mathcal{F}}$ acting separately on the physical and horizontal virtual legs of one column of A sites of the hybrid PEPS. The result defines the column tensor $A^{\text{col}}$ for the real-space PEPS. (c) Decomposition of the column tensor $A^{\text{col}}$ into the contraction of local PEPS tensors $A^{2D}$ for the two-dimensional state. The vertical virtual legs of $A^{2D}$ marked in red have a bond dimension $D_{y,A}$ which generally grows with the system size $N_y$ due to the nonlocality of the inverse FT. The inverse FT and decomposition of a column of B sites is analogous.</figcaption></figure>

by the local tensor of the $d$-dimensional model at the $d$-dimensional real-space position $\mathbf{x}$ and time $t^{(j)}$.

Thus, the virtual DOFs of the hybrid TNS in the first $d$ directions are identical to those of the $d$-dimensional model. In other words, the identification of time and momentum from Eq. (19) holds not just for the physical mode operators, but also for the virtual mode operators in the first $d$ dimensions.

On the other hand, due to the direct product in the direction of $k_{d+1}$ in Eq. (20), the hybrid TNS does not need virtual legs in the direction $d+1$, and we say that the bond dimension in this direction is equal to one (implying that the contraction of this direction corresponds trivially to a product).

The hybrid Chern PEPS obtained from the SSH charge pumping MPS of Eq. (9) is sketched in Fig. 3(a) in the hybrid coordinate system where the horizontal axis describes $\mathbf{x}$ and the vertical axis corresponds to $k_y \equiv k_{d+1}$. The local tensor for a site on the A sublattice at the position $(\mathbf{x}, k_y^{(j)})$ is given by the SSH local tensor $A_{lr}^i(t^{(j)})$ at the time $t^{(j)}$ along the cycle, and similarly for the B sublattice. In the horizontal direction, the two-dimensional hybrid state hence inherits both the translation invariance and the constant finite bond dimension $D_{SSH}$ of the MPS. As discussed in the previous paragraph, virtual legs in the vertical direction are not needed for the hybrid Chern PEPS and are hence not shown in Fig. 3(a).

### C. Inverse Fourier transform

As discussed below Eq. (20), the hybrid state obtained from charge pumping is mapped to a $(d+1)$-dimensional real-space coordinate system by applying an inverse FT in the direction $d+1$ to the fermionic mode operators. For GfTNSs, which have virtual in addition to physical DOFs, the inverse FT should be applied to both the physical mode operators and the virtual fermionic mode operators in the first $d$ directions. For the hybrid Chern PEPS, we therefore apply the inverse FT in the vertical direction to the physical and the horizontal virtual legs. The extension of the FT to the horizontal virtual

modes amounts to a virtual basis change and does not alter the physical state, but ensures that the real-space PEPS is invariant under vertical translations $y \mapsto y+1$ of its physical and virtual legs (see also Appendix F).

For the hybrid Chern PEPS, the correct definition of the inverse FT in the vertical direction $y$ entails a subtlety. Indeed, recall that the SSH charge-pumping MPS from Eq. (9) is defined with respect to mode operators $a_{\tau,x}(t^{(j)})$ related to the $\hat{a}_{\tau,x}(t^{(j)})$ used in Eq. (19) by the particle-hole transformation of Eq. (7). For the Chern PEPS, we perform an analogous particle-hole transformation in two-dimensional real space and define new mode operators

$$a_{A,(x,y)} = \hat{a}_{A,(x,y)}, \qquad (21a)$$

$$a_{B,(x,y)} = \hat{a}_{B,(x,y)}^{\dagger}. \qquad (21b)$$

The operators $a_{A,(x,y)}$, $a_{A,(x,y)}^{\dagger}$, $a_{B,(x,y)}$, $a_{B,(x,y)}^{\dagger}$ span the basis in which we want to express the real-space Chern PEPS.

The SSH particle-hole transformation of Eq. (7) acts within the set of operators at time $t^{(j)}$. Since we identify time and momentum, this is equivalent to a particle-hole transformation acting on the modes of the Chern PEPS in momentum space rather than real space like in Eq. (21). Due to the antiunitarity of the particle-hole transformation, the inverse FT relating the modes $a_{B,x}(t^{(j)})$ and $a_{B,(x,y)}$ on the B sublattice therefore contains an additional complex conjugation. Hence, the physical mode operators $a_{\tau,x}(t^{(j)})$ of the SSH models, providing the basis for the hybrid Chern PEPS, and the mode operators $a_{\tau,(x,y)}$ for the real-space Chern PEPS are related as

$$a_{\tau,(x,y)} = \sum_{j=0}^{N_y-1} \tilde{\mathcal{F}}_{\tau,\tau;y,t^{(j)}} a_{\tau,x}(t^{(j)}) \qquad (22a)$$

with the inverse vertical FT

$$\tilde{\mathcal{F}}_{\tau, \tau'; y, t^{(j)}} = \frac{\delta_{\tau, \tau'}}{\sqrt{N_y}} e^{i \eta_{\tau, y} t^{(j)}}. \qquad (22b)$$
Here, $\eta_A = 1$ and $\eta_B = -1$ for the physical modes which are particlelike on the A sublattice and holelike on the B sublattice.

The inverse FT of the horizontal virtual legs of the hybrid PEPS is analogous to Eq. (22) for the physical modes. Here, the particlelike or holelike character of the virtual modes is determined by their transformation under the U(1) symmetry of Eq. (12). Specifically, the left and right virtual legs on the A sublattice transform as holes, i.e., $\eta_{L,A} = \eta_{R,A} = -1$, whereas the left and right virtual legs on the B sublattice transform as particles such that $\eta_{L,B} = \eta_{R,B} = 1$.

Let us now study how the inverse vertical FT $\tilde{\mathcal{F}}$ acts on the hybrid Chern PEPS. Due to the translation invariance in the horizontal direction, it is enough to consider one column of the hybrid state given by the sites $(x, k_y^{(j)})$ for a fixed horizontal position $x$ and $j = 0, \dots, N_y - 1$ [see Fig. 3(a)]. Under $\tilde{\mathcal{F}}$, the hybrid column is mapped to one column of the two-dimensional real-space state, given by the sites $(x, y)$ with $y = 0, \dots, N_y - 1$.

Since the inverse FT $\tilde{\mathcal{F}}$ is nonlocal, for a generic interpolation there are long-range correlations in the real-space column $\{(x, y)\}_{0\le y \le N_y - 1}$ of the PEPS. In the tensor network language, the real-space column is therefore described by a single tensor obtained from the application of $\tilde{\mathcal{F}}$ to one column of the hybrid PEPS [in contrast to the hybrid column $\{(x, k_y^{(j)})\}_{0\le j \le N_y - 1}$ which is described by a product of individual tensors for each $k_y^{(j)}$, signifying the absence of correlations in the vertical direction]. For a real-space column of A sites, we denote this tensor by $A^{\text{col}}$ as shown in Fig. 3(b). Similarly, we define a tensor $B^{\text{col}}$ for a real-space column of B sites. The tensors $A^{\text{col}}$ and $B^{\text{col}}$ have $N_y$ physical legs of dimension $d_p = 2$ and $N_y$ left and right virtual legs of dimension $D_{\text{SSH}} = 2$, which describe the physical and horizontal virtual DOFs of the real-space PEPS.

The SSH pumping MPS is a Gaussian MPS for free fermions. Thus, the column tensors $A^{\text{col}}$ and $B^{\text{col}}$ are also Gaussian states, which can be described by their covariance matrices (CMs). These CMs are computed in Appendix F for the general $d$-dimensional case. There we show that the CMs of the real-space columns $A^{\text{col}}$ and $B^{\text{col}}$, defined by the application of the inverse FT of Eq. (22) to the physical and horizontal virtual legs of the hybrid columns, are given by the inverse FT of the time-dependent CMs describing the local SSH tensors $A_{lr}^i(t^{(j)})$ and $B_{lr}^i(t^{(j)})$ along the pumping cycle. This result fully characterizes the column tensors $A^{\text{col}}$ and $B^{\text{col}}$.

## D. Translation-invariant real-space PEPS

From the previous subsections, we are now ready to discuss how the two-dimensional real-space state can be expressed as a TNS with local tensors for each site. For that purpose, the column tensors $A^{\text{col}}$ and $B^{\text{col}}$ have to be decomposed into a column of local PEPS tensors $A^{\text{2D}}$ and $B^{\text{2D}}$ with both horizontal virtual legs of dimension $D_{\text{SSH}}$ and vertical virtual legs with dimension $D_{y,A}$ and $D_{y,B}$, respectively. We require the local PEPS tensors to be identical on all sites of the same sublattice as shown in Fig. 3(c) in order to make the invariance under real-space translations explicit.

The decomposition of the column into PEPS tensors can be achieved using tools developed for one-dimensional MPSs [53,54]. To that end, we define one-dimensional pure states $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$, whose many-body wave functions are the tensor elements of $A^{\text{col}}$ and $B^{\text{col}}$, respectively. Therefore, $|\psi_{1D}(A^{\text{col}})\rangle$ effectively describes a chain of length $N_y$, where the local Hilbert space at each site of the fictitious chain is the tensor product of the physical and horizontal virtual Hilbert spaces of $A^{\text{col}}$ at the corresponding physical site (and similarly for $|\psi_{1D}(B^{\text{col}})\rangle$ and $|B^{\text{col}}\rangle$). Pictorially, $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$ are obtained by merging at each site the physical and horizontal virtual legs of $A^{\text{col}}$ and $B^{\text{col}}$ into an effective physical index of dimension $\tilde{d}_p = d_p D_{\text{SSH}}^2$ per site.

Our goal is now to express $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$ as translation-invariant MPSs with periodic boundary conditions. Then, after unfolding the effective physical index into the physical and horizontal virtual indices of the PEPS, the local tensors of these MPSs will define the PEPS tensors $A^{\text{2D}}$ and $B^{\text{2D}}$, respectively. Similarly, the MPS bond dimensions $D_{y,A}$ and $D_{y,B}$ are equal to the vertical bond dimensions of the PEPS local tensors $A^{\text{2D}}$ and $B^{\text{2D}}$, respectively.

Numerically, the translation-invariant MPSs for $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$ can be obtained as follows. Multiple steps are necessary since the common and stable MPS algorithms rely on the presence of open boundary conditions and do not directly find a periodic translation-invariant MPS for a translation-invariant state. Instead, in a first step the Gaussian states $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$ can be decomposed into non-translation-invariant Gaussian MPSs with open boundary conditions by performing successive Schmidt decompositions of the system [53,55]. We are not aware of any method to transform a Gaussian open-boundary MPS into a Gaussian translation-invariant MPS. We therefore interpret the Gaussian MPSs for $|\psi_{1D}(A^{\text{col}})\rangle$ and $|\psi_{1D}(B^{\text{col}})\rangle$ as regular fermionic open-boundary MPSs by choosing physical and virtual basis states and computing the local MPS tensors for each site. Finally, these open-boundary MPSs are transformed into translation-invariant and generically non-Gaussian MPSs with periodic boundary conditions following the standard procedure described in Ref. [54]. This approach leads to the bond dimension of the translation-invariant MPS being given by the sum of the bond dimensions of the open-boundary MPS for each site. Therefore, the bond dimension grows at least linearly with the system size [54].

## E. Vertical PEPS bond dimension

Above, we showed that the inverse FT of the hybrid PEPS can be decomposed into a $(d+1)$-dimensional real-space TNS, where the bond dimension in the first $d$ directions is equal to that of the original $d$-dimensional state. The core question is the scaling of the bond dimension in the $(d+1)^d$ real-space dimension with respect to the system size in this direction.

### 1. Lower bound from ES

The bond dimension of a tensor network is intimately related to the entanglement between its physical particles: For
<figure><img src="image_6.png" /><figcaption>FIG. 4. Schmidt decomposition of the one-dimensional effective state $|\psi_{1D}\rangle$ for the PEPS column with physical dimension $\tilde{d}_p$ into the subsystem $A_L$ of the first $L$ sites and its complement. When $|\psi_{1D}\rangle$ is written as a translation-invariant MPS with bond dimension $D_y$, the virtual boundary of the subsystem $A_L$, marked in blue, crosses two bonds of the MPS.</figcaption></figure>

a TNS describing a physical system with a subsystem $A$, the total dimension of all virtual legs at the boundary $\partial A$ can be no smaller than the rank of the state’s Schmidt decomposition into the DOFs of $A$ and its complement. This constraint allows us to infer a lower bound for the vertical bond dimension of the real-space PEPS from the ES of the one-dimensional column states.

We consider a generic one-dimensional effective column state denoted $|\psi_{1D}\rangle$ with an effective physical dimension $\tilde{d}_p$. Below, we will choose either $|\psi_{1D}\rangle = |\psi_{1D}(A^{\text{col}})\rangle$ or $|\psi_{1D}\rangle = |\psi_{1D}(B^{\text{col}})\rangle$ with $\tilde{d}_p = d_p D_{\text{SSH}}^2$. We assume that $|\psi_{1D}\rangle$ is expressed as a translation-invariant MPS with a bond dimension $D_y$ corresponding to the vertical PEPS bond dimension (see Fig. 4). Such a representation can be obtained for instance using the procedure described in the previous subsection. We now perform a Schmidt decomposition of the pure state $|\psi_{1D}\rangle$ with respect to the subsystem $A_L$ of the first $L$ sites $0 \le y \le L-1$, and denote the rank of the decomposition by $r_L(|\psi_{1D}\rangle)$. As sketched in Fig. 4, the cut between $A_L$ and its complement crosses exactly two virtual bonds of the MPS, namely, those between sites $(L-1, L)$ and $(N_y-1, 0)$, each of dimension $D_y$. Since the virtual dimension can be no smaller than the Schmidt rank for any subsystem size $L$, we obtain the bound

$$D_y^2 \ge \max_{1 \le L \le N_y-1} r_L(|\psi_{1D}\rangle). \quad (23)$$

The Schmidt rank $r_L(|\psi_{1D}\rangle)$ may be obtained by counting the nontrivial levels in the single-particle ES of the Gaussian column state $|\psi_{1D}\rangle$. This ES can be computed numerically based on the explicit result for the CM of $|\psi_{1D}\rangle$ given in Appendix F. Since the state $|\psi_{1D}\rangle$ with physical dimension $\tilde{d}_p$ has $\log_2(\tilde{d}_p)$ free-fermionic modes per lattice site, its single-particle ES with respect to the subsystem $A_L$ consists of levels

$0 \le |\lambda_i(L)| \le 1$ with $1 \le i \le L \log_2(\tilde{d}_p)$ [see Eq. (A5) for a definition of the single-particle ES and its levels $|\lambda_j|$. For an improved numerical stability, we actually compute $\mu_i(L) = \sqrt{1 - \lambda_i(L)^2}$ that may be conveniently extracted from the CM [see Eq. (A6)]. Here, values $\mu_i(L) = 1$ and $\mu_i(L) = 0$ correspond to maximally entangled and nonentangled modes, respectively. The many-body Schmidt rank of the state $|\psi_{1D}\rangle$ is therefore given by the exponential

$$r_L(|\psi_{1D}\rangle) = 2^{\#\{\mu_i(L) > 0\}} \quad (24)$$

of the number of entangled modes in the single-particle ES with $\mu_i(L)$ bigger than zero.

## 2. Exponential growth

The number of entangled modes in the single-particle ES of a column $A^{\text{col}}$ of unit cells of the real-space PEPS is displayed in Fig. 5(a) for two different cyclic interpolations of the MPS from Eq. (11): on one hand, the parametrization $\phi_{\text{pump}}$ of Eq. (15) giving rise to a two-dimensional Chern insulator, and on the other hand the parametrization $\phi_{\text{triv}}$ of Eq. (17) corresponding to a topologically trivial two-dimensional state. The corresponding data of the column $B^{\text{col}}$ is identical.

For $\phi_{\text{pump}}$ the number of entangled modes is equal to

$$\#\{\mu_i(L) > 0\} = \min\{3L, N_y\} \quad (25)$$

when $L \le N_y/2$ (the spectra for $L$ and $N_y - L$ are identical). In Appendix G we show that this is the maximal number of entangled modes which is compatible with the global U(1) symmetry of $A^{\text{col}}$ inherited from the SSH model MPS. Here, the factor 3 in Eq. (25) is a consequence of the column tensor having three fermions (one physical fermion and two virtual fermions) per site. The validity of Eq. (25) can be seen in Fig. 5(a) for the smallest system size $N_y = 28$. For larger system sizes, the number of entangled modes shown in Fig. 5(a) is lower than Eq. (25) because of the finite numerical resolution, as shown by the ES in Fig. 5(b). We have checked that for the Chern PEPS obtained from the MPS corresponding to the smooth charge-pumping cycle described in Sec. II A, the number of entangled modes is also given by Eq. (25).

For $\phi_{\text{triv}}$, the number of entangled modes is given by $2L$ when $L \le N_y/2$. This is lower than Eq. (25) for $\phi_{\text{pump}}$ except at $L = N_y/2$, where both numbers agree. As discussed in Appendix G, this reduction of entangled modes is due to the decoupling of all right virtual legs and is not related to topology. Indeed, for the trivial cycle, the parameter $\beta(t) = 0$ vanishes for the entire duration of the interpolation $\phi_{\text{triv}}$. Any small perturbation of $\phi_{\text{triv}}$ by a nonzero $\beta$ increases the number of entangled modes to Eq. (25). Moreover, these additional decoupled modes disappear when blocking the two columns $A^{\text{col}}$ and $B^{\text{col}}$, i.e., by contracting the bonds connecting them. Hence, the column for AB unit cells has the same Schmidt rank for both $\phi_{\text{pump}}$ and $\phi_{\text{triv}}$.

The number of entangled modes from Eq. (25) corresponds to a maximal Schmidt rank

$$\max_{1 \le L \le N_y-1} r_L(|\psi_{1D}(A^{\text{col}})\rangle) = 2^{N_y} \quad (26)$$
<figure><img src="image_7.png" /><figcaption>FIG. 5. Single-particle ES of a column Acol of the two-dimensional real-space TNSs defined by two different parametrizations for the MPS of Eq. (11): the charge-pumping interpolation $\phi_{\text{pump}}$ of Eq. (15) leading to a Chern insulator, and the parametrization $\phi_{\text{triv}}$ of Eq. (17) leading to a topologically trivial two-dimensional state. (a) Number of modes in the single-particle ES with a finite entanglement corresponding to a value $\mu_i > \delta$ with $\delta = 10^{-13}$ the numerical accuracy for different system sizes $N_y$ as a function of the subsystem size $L$. (b) Single-particle ES for a column with $N_y = 256$ sites with respect to the subsystem $A_L$ with $L = 20$. The double degeneracy in the ES corresponding to $\phi_{\text{triv}}$ is due to the decoupling of all right virtual legs since $\beta(t) = 0$ along $\phi_{\text{triv}}$.</figcaption></figure>

for the state $|\psi_{1D}(A^{\text{col}})\rangle$, and similarly for $|\psi_{1D}(B^{\text{col}})\rangle$. The bound of Eq. (23) therefore implies that the vertical bond dimensions $D_{y,A}$ and $D_{y,B}$ grow exponentially as a function of $N_y$. This is a generic feature of our construction and not related to the topology of the state: It originates in the nonlocality of the inverse FT which couples all states in the exponentially growing Hilbert space of one column. Indeed, a generic quantum state requires an exponentially growing bond dimension to be represented *exactly* as a TNS.

Despite the faster decay of the single-particle entanglement energies for the topologically trivial state than for the Chern insulator in Fig. 5(b), this does not allow us to make any statement about differences in the growth of the vertical bond dimension required for an *approximative* PEPS for the two systems.

## IV. FERMIONIC PEPS FOR TWO-DIMENSIONAL HIGHER-ORDER TI

In this section, we study a Gaussian fermionic PEPS for the topological quadrupole model from Ref. [1], which is reviewed in Sec. IV A. The PEPS is defined in Sec. IV B, where we also provide its parent Hamiltonian. In Sec. IV C we discuss a three-dimensional PEPS with chiral hinge states obtained from a dipole-pumping interpolation of the quadrupole model [2,3].

### A. Second-order quadrupole insulator

A two-dimensional second-order topological phase with a quantized bulk quadrupole moment was recently proposed theoretically [1] and has subsequently been observed experimentally in mechanical [11], photonic [14–17], and electrical [18,19] systems. This phase is described by a microscopic free-fermionic model with one spinless fermionic mode per

site and a unit cell of $2 \times 2$ sites depicted in Fig. 6(a). We consider the system at half-filling where only the lowest two bands are occupied. With open boundary conditions, the nearest-neighbor Hamiltonian is given by

$$H_{\text{Quad}} = \sum_{\mathbf{x}} \sum_{j=0,1} \left[ t_{\mathbf{x}}^{(j)} (\hat{a}_{1,\mathbf{x}}^{\dagger} \hat{a}_{3,\mathbf{x}+j\hat{\mathbf{x}}} + \hat{a}_{4,\mathbf{x}}^{\dagger} \hat{a}_{2,\mathbf{x}+j\hat{\mathbf{x}}}) + t_{\mathbf{y}}^{(j)} (\hat{a}_{1,\mathbf{x}}^{\dagger} \hat{a}_{4,\mathbf{x}+j\hat{\mathbf{y}}} - \hat{a}_{3,\mathbf{x}}^{\dagger} \hat{a}_{2,\mathbf{x}+j\hat{\mathbf{y}}}) + \text{H.c.} \right], \quad (27)$$

where $\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ denote the unit vectors in the horizontal and vertical direction, respectively. The positions of the unit cells are $\mathbf{x} = x\hat{\mathbf{x}} + y\hat{\mathbf{y}}$ with $0 \le x \le N_x - 1$ and $0 \le y \le N_y - 1$ on a lattice with $N_x$ and $N_y$ unit cells in the horizontal and vertical directions. For $\tau = 1, \dots, 4$, $\hat{a}_{t,\mathbf{x}}^{\dagger}$ denotes the creation operator for a fermion in orbital $\tau$ and unit cell at position $\mathbf{x}$. The Hamiltonian of Eq. (27) contains four nearest-neighbor hopping amplitudes $t_x^{(0)}, t_y^{(0)}, t_x^{(1)},$ and $t_y^{(1)}$, where the subscripts $x$ and $y$ refer to hoppings in the horizontal and vertical directions, and the superscripts (0) and (1) indicate hopping between sites on the same and adjacent unit cells, respectively. The signs of the hopping amplitudes ensure that there is a flux $\pi$ through every plaquette of the square lattice.

The quadrupole model of Eq. (27) is in a second-order topological OAI phase protected by the horizontal and vertical mirror symmetries $M_x$ and $M_y$ if $t_x^{(0)}/t_x^{(1)} \in [-11]$ and $t_y^{(0)}/t_y^{(1)} \in [-11]$ [1,56]. In this phase, corners host gapless protected states at the intersection of two gapped edges. Moreover, the system possesses a quantized bulk quadrupole moment, edge dipole moment, and corner charge [3]. When the hopping amplitudes $t_x^{(0)}$ and $t_y^{(0)}$ within unit cells vanish in the OAI phase, the model is in a dimerized OAI phase where every site in the bulk is contained in a decoupled plaquette shifted from the unit cell by one site in the horizontal and vertical directions. Hence, the corner modes are fully localized on the four corner sites of a rectangular patch. On
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Sketch of the quadrupole model from Ref. [3]. A unit cell (marked with a green square) contains four sites. The horizontal and vertical nearest-neighbor hoppings $t_x^{(0)}$ and $t_y^{(0)}$ within unit cells are sketched in black, whereas the nearest-neighbor hoppings $t_x^{(1)}$ and $t_y^{(1)}$ between unit cells are sketched in red. The couplings corresponding to dashed bonds carry a negative sign to ensure a flux $\pi$ through every plaquette. Lattice sites marked in blue and red have a chemical potential $\mu$ and $-\mu$, respectively. (b) Sketch of a unit cell of the quadrupole PEPS from Eq. (29). For each site 1, 2, 3, or 4 marked by a dashed blue square, the blue (red) circles denote its physical (virtual) fermions. Each edge connecting a physical and virtual fermion is labeled with the amplitude of the local tensor when these two fermions are in state $|1\rangle$ and all other virtual fermions of the same site are in state $|0\rangle$. (c) Next-to-nearest-neighbor hopping terms of Eq. (32) in the parent Hamiltonian $H_{PEPS}$.</figcaption></figure>

the other hand, for $|t_x^{(0)}/t_x^{(1)}| > 1$ or $|t_y^{(0)}/t_y^{(1)}| > 1$ the model of Eq. (27) is in a trivial phase with gapped edges and corners. If the hoppings $t_x^{(1)}$ and $t_y^{(1)}$ between adjacent unit cells vanish, the system is in a trivial dimerized phase where each unit cell decouples from the rest of the system.

Similarly to the charge-pumping cycle discussed in Sec. II A, a dipole-pumping cycle interpolating between the trivial and OAI dimerized phases of the quadrupole model can be defined by adding the chemical potential $\mu$ to all sites on the sublattices 1 and 2 and $-\mu$ to all sites on the sublattices 3 and 4 [3]. Therefore, the staggering pattern breaks the symmetries protecting the topological phase but preserves the $C_2$ rotation symmetry [see Fig. 6(a)]. The dipole-pumping cycle is obtained from $C_4$-symmetric hopping amplitudes $t^{(0)} = t_x^{(0)} = t_y^{(0)}$ and $t^{(1)} = t_x^{(1)} = t_y^{(1)}$ and a chemical potential $\mu$ evolving according to the same cyclic interpolation of Eq. (3) as for the SSH charge pump.

The properties of the dipole-pumping cycle for the quadrupole model are analogous to those of the charge-pumping cycle for the SSH model. In particular, the system remains in a dimerized state throughout the interpolation since at all times $t \in [-\pi, \pi]$ only one of the two hopping amplitudes $t^{(0)}(t)$ and $t^{(1)}(t)$ is nonzero. At $t = \pm \pi$ and $t = 0$, the system is in an atomic state with only the orbitals 1 and 2 occupied at $t = \pm \pi$ and only the orbitals 3 and 4 occupied at $t = 0$. For $-\pi < t < 0$, the hopping within unit cells is nonzero, whereas the hopping between unit cells is nonzero for $0 < t < \pi$. In both cases, charge gets transferred by the changing chemical potential. At $t = -\pi/2$ and $\pi/2$, the chemical potential vanishes such that the mirror symmetries are restored and the Hamiltonian corresponds to the trivial and OAI dimerized phase of the quadrupole model, respectively.

In the same manner as charge pumping relates the SSH model to a Chern insulator, dipole pumping induces a model with chiral hinge states from the quadrupole model [2,3]. The chiral hinge model is a three-dimensional second-order

topological insulator whose one-dimensional protected boundary modes occur at the intersection of a pair of two-dimensional faces. The topology of the hinge model obeys a $\mathbb{Z}_2$ classification protected by the product $M_x M_y T$ of the horizontal and vertical mirror symmetries and time reversal $T$ [2].

## B. PEPS for the quadrupole model

In this section, we provide a Gaussian fermionic PEPS for the ground state of the quadrupole model. After giving the details of the construction in Sec. IV B 1, we compute its parent Hamiltonian in Sec. IV B 2 and discuss its ES in Sec. IV B 3.

### 1. Construction

In order to construct a Gaussian PEPS for the ground state of the quadrupole model, we use similar ideas as in the SSH model MPS derived in Sec. II. In analogy to the particle-hole transformation of Eq. (7) on the B sublattice of the SSH chain, we define new physical modes $a_{\tau,x}$ by performing a particle-hole transformation on the sublattices 3 and 4 while leaving the sublattices 1 and 2 unaffected,

$$a_{1,x} = \hat{a}_{1,x}, \quad a_{2,x} = \hat{a}_{2,x}, \qquad (28a)$$

$$a_{3,x} = \hat{a}_{3,x}^{\dagger}, \quad a_{4,x} = \hat{a}_{4,x}^{\dagger}. \qquad (28b)$$

The vacuum $|\Omega\rangle$ of the new modes, satisfying $a_{\tau,x}|\Omega\rangle = 0$ for $\tau = 1, \dots, 4$, contains exactly $2N_x N_y$ physical particles as required for the quadrupole model at half-filling. Thus, the particle-hole transformation allows us to express the quadrupole PEPS in terms of a separate and parity-even local tensor for each lattice site. Parity evenness is required in order to ensure that the state is independent of the order of contractions in the network [57].
The quadrupole PEPS has bond dimension $D_{Quad} = 2$ and is constructed from four local tensors $A^{[\tau]i}_{\tau \text{lurd}}$ for sites on the four sublattices $\tau = 1, \dots, 4$, respectively. Here, $i \in \{0, 1\}$ corresponds to the physical index and $l, u, r, d \in \{0, 1\}$ to the left, up, right, and down virtual indices, respectively. The physical basis states $|0\rangle = |\Omega\rangle$ and $|1\rangle = a_{\tau,x}^{\dagger}|\Omega\rangle$ are obtained from the transformed mode operators of Eq. (28).

Due to the close relation between the SSH and the quadrupole models, we are guided in our ansatz for the local tensors of the quadrupole PEPS by the MPS tensors of Eq. (11). This can be most easily seen from Fig. 6(b). We obtain the quadrupole model by coupling neighboring sites both horizontally and vertically according to the pattern of an SSH chain. In the PEPS tensors, the couplings of the horizontal SSH chains are transmitted by the left and right virtual fermions associated to parameters $\alpha_x, \beta_x$ analogous to Eq. (11). Similarly, the vertical SSH chains are implemented using the top and bottom virtual fermions associated to the parameters $\alpha_y, \beta_y$. Using the analogy to Eq. (11), we therefore obtain local PEPS tensors whose nonvanishing elements are given in terms of five real parameters $\gamma, \alpha_x, \alpha_y, \beta_x$, and $\beta_y$ as

$$A^{[1]0}_{0000} = A^{[2]0}_{0000} = A^{[3]0}_{0000} = A^{[4]0}_{0000} = \gamma, \quad (29a)$$

$$-A^{[1]1}_{1000} = A^{[2]1}_{1001} = A^{[3]1}_{1001} = A^{[4]1}_{1000} = \beta_x, \quad (29b)$$

$$A^{[1]1}_{0001} = A^{[2]1}_{0100} = A^{[3]1}_{0001} = A^{[4]1}_{0100} = \beta_y, \quad (29c)$$

$$A^{[1]1}_{0010} = -A^{[2]1}_{1000} = A^{[3]1}_{1000} = A^{[4]1}_{0010} = \alpha_x, \quad (29d)$$

$$-A^{[1]1}_{0100} = -A^{[2]1}_{0001} = A^{[3]1}_{0100} = A^{[4]1}_{0000} = \alpha_y. \quad (29e)$$

All the other tensor elements are equal to zero. The phases in Eq. (29) were chosen such as to ensure that there is a flux $\pi$ through every plaquette. As sketched in Fig. 6(b), the parameters $\beta_x$ and $\beta_y$ represent the coupling of the physical leg to the virtual legs corresponding to the horizontal and vertical bonds pointing into the unit cell, respectively. Similarly, $\alpha_x$ and $\alpha_y$ control the coupling of the physical leg to the virtual legs pointing out of the unit cell.

The PEPS of Eq. (29) has a global U(1) symmetry analogous to Eq. (12) for the SSH MPS. Indeed, the local tensors are invariant under a combination of U(1) rotations of the physical and virtual legs given by

$$A^{[\tau]i}_{\text{lurd}} = \sum_j \sum_{l'u'r'd'} U(\varphi)_{ij} U(\varphi)_{ll'}^{\dagger} U(\varphi)_{uu'}^{\dagger} \\ \times A^{[\tau]j}_{l'u'r'd'} U(\varphi)_{r'r}^{\dagger} U(\varphi)_{d'd}^{\dagger} \quad (30a)$$

for sites on the sublattices $\tau = 1, 2$, and

$$A^{[\tau]i}_{\text{lurd}} = \sum_j \sum_{l'u'r'd'} U(\varphi)_{ij}^{\dagger} U(\varphi)_{ll'} U(\varphi)_{uu'} \\ \times A^{[\tau]j}_{l'u'r'd'} U(\varphi)_{r'r} U(\varphi)_{d'd} \quad (30b)$$

for sites on the sublattices $\tau = 3, 4$. Here, $U(\varphi) = \begin{pmatrix} 1 & 0 \\ 0 & e^{i\psi} \end{pmatrix}$ is the U(1) rotation acting on a single fermion. Equation (30) implies that the virtual legs for the local tensors on the sublattices $\tau = 1, 2$ and $\tau = 3, 4$ transform as holes and particles. Hence, each pair of virtual legs associated with the same nearest-neighbor bond transforms oppositely under

the U(1) rotation, such that the PEPS is invariant under the physical part of Eq. (30). The charge associated with this symmetry is

$$\sum_{\mathbf{x}} \left[ \sum_{\tau=1,2} a_{\tau, \mathbf{x}}^{\dagger} a_{\tau, \mathbf{x}} - \sum_{\tau=3,4} a_{\tau, \mathbf{x}}^{\dagger} a_{\tau, \mathbf{x}} \right] \\ = \sum_{\mathbf{x}} \sum_{\tau=1}^{4} \hat{a}_{\tau, \mathbf{x}}^{\dagger} \hat{a}_{\tau, \mathbf{x}} - \frac{1}{2} N_{\mathbf{x}} N_{\mathbf{y}} \quad (31)$$

such that the U(1) symmetry ensures that the state lies exactly at half-filling of the lattice just as in the one-dimensional case.

## 2. Parent Hamiltonian

For all values of the parameters $\gamma, \alpha_x, \alpha_y, \beta_x$, and $\beta_y$, the PEPS from Eq. (29) can be expressed as a Gaussian TNS for free fermions. In Appendix H, we show that the PEPS with parameters $\alpha = \alpha_x = \alpha_y$ and $\beta = \beta_x = \beta_y$ is the unique ground state of an extended version of the quadrupole model $H_{\text{PEPS}}$ with $C_4$-symmetric hoppings, with a staggered chemical potential that breaks $C_4$ symmetry and with an additional next-to-nearest-neighbor hopping

$$\sum_{\tau} \sum_{\mathbf{x}} t_{\tau}^{(2)} [\hat{a}_{\tau, \mathbf{x}}^{\dagger} \hat{a}_{\tau, \mathbf{x}+\hat{s}} + \hat{a}_{\tau, \mathbf{x}}^{\dagger} \hat{a}_{\tau, \mathbf{x}+\hat{s}} + \text{H.c.}] \quad (32)$$

between sites on the same sublattice $\tau$ with amplitude $t_{\tau}^{(2)}$. The pattern of next-to-nearest-neighbor hoppings is shown in Fig. 6(c). The couplings of $H_{\text{PEPS}}$ depend on the parameters of the PEPS as

$$\mu = \frac{\gamma^4 - \alpha^4 - \beta^4}{\alpha^4 + \beta^4 + \gamma^4}, \quad (33a)$$

$$t_x^{(0)} = t_y^{(0)} = \frac{\sqrt{2}\beta^2\gamma^2}{\alpha^4 + \beta^4 + \gamma^4}, \quad (33b)$$

$$t_x^{(1)} = t_y^{(1)} = \frac{\sqrt{2}\alpha^2\gamma^2}{\alpha^4 + \beta^4 + \gamma^4}, \quad (33c)$$

$$t_1^{(2)} = t_2^{(2)} = -t_3^{(2)} = -t_4^{(2)} = -\frac{1}{2} \frac{\alpha^2 \beta^2}{\alpha^4 + \beta^4 + \gamma^4}. \quad (33d)$$

Similarly to the parent Hamiltonian $H_{\text{MPS}}$ of the SSH model MPS of Eq. (11), $H_{\text{PEPS}}$ describes different phases depending on the values of the parameters $\alpha, \beta$, and $\gamma$. When two out of the three parameters vanish, the system is in an atomic insulator state: If $\alpha = \beta = 0$, the particles are localized on the orbitals $\tau = 3, 4$, whereas they are localized on the orbitals $\tau = 1, 2$ if $\alpha = \gamma = 0$ or $\beta = \gamma = 0$. In contrast, when $\alpha = 0$ and $\beta, \gamma \neq 0$, all hoppings of the parent Hamiltonian vanish except for the nearest-neighbor hopping within unit cells. Similarly, when $\beta = 0$ and $\alpha, \gamma \neq 0$, the only nonzero hopping is the nearest-neighbor hopping between unit cells. In both cases, the system is in a dimerized phase with a staggered chemical potential. Setting $\beta = \gamma$ and $\alpha = \gamma$, respectively, we recover the trivial and OAI dimerized phases of the quadrupole model with vanishing chemical potential. Finally, when all three parameters are nonzero, the system has a nonvanishing nearest- and
next-nearest-neighbor hopping as well as a finite staggered chemical potential.

### 3. Corner states and entanglement spectrum

The characteristic ($d-1$)-dimensional gapless edge states of a conventional $d$-dimensional TI are reflected in the state’s bulk ES [22]. Similarly, the ES of the quadrupole model in its dimerized OAI phase hosts gapless corner states as long as the virtual cut is compatible with the protecting symmetries [24]. In order to further characterize the different phases described by the PEPS from Eq. (29), we therefore study the ES of the state defined on a torus with respect to a rectangular subsystem $A_{L_x \times L_y}$ of $L_x$ and $L_y$ unit cells in the horizontal and vertical direction.

For $\beta = 0$, the PEPS from Eq. (29) is in a dimerized phase where the system splits into four-site plaquettes shifted from the unit cell by one site in the horizontal and vertical directions. For $\alpha = \gamma$, we obtain the OAI dimerized phase of the quadrupole model, whereas for $\alpha \neq \gamma$ there is a nonzero staggered chemical potential that breaks the symmetries protecting the OAI phase. Due to the dimerization, with open boundaries the system has SSH chains with a staggered chemical potential at the edges, and corner sites which decouple from each other and the bulk. Correspondingly, the single-particle ES of the PEPS with respect to the subsystem $A_{L_x \times L_y}$ has three distinct contributions from the bulk, the edges, and the corners as sketched in Fig. 7. The bulk consists of $(L_x - 1)(L_y - 1)$ plaquettes decoupled from the rest of the system, which contribute $4(L_x - 1)(L_y - 1)$ nonentangled modes with levels $|\lambda_{\text{bulk}}| = 1$ (in the dimerized limit) to the single-particle ES.

On the other hand, the corner and edge sites belong to plaquettes crossed by the boundary of the subsystem $A_{L_x \times L_y}$. In Appendix H 3 we show that the four boundaries pro-

<figure><img src="image_9.png" /><figcaption>FIG. 7. Quadrupole model in the dimerized OAI phase with a rectangular subsystem $A$ marked by a blue rectangle. The four-site plaquettes coupled by the nearest-neighbor hopping $t^{(1)}$ (drawn in red) are shifted from the unit cells (denoted by green squares) by one site in both directions. The decoupled plaquettes in the bulk of $A$ (marked in blue) do not contribute to the ES. In contrast, nontrivial levels in the ES come from both the plaquettes at the corners with a single site in $A$ (marked in green), and from those at the edges with two sites in $A$ (marked in orange).</figcaption></figure>

vide $4(L_x - 1 + L_y - 1)$ entangled modes with levels $\pm \lambda_{\text{edge}}$, where

$$ \lambda_{\text{edge}} = \frac{\sqrt{\gamma^8 + \alpha^8}}{\gamma^4 + \alpha^4}. \qquad (34a) $$

Since the number of these entangled modes grows linearly with the size of the edge, the boundary leads to an area-law term in the entanglement entropy.

Furthermore, each corner site hosts exactly one mode corresponding to a level $\mp \lambda_{\text{corner}}$ with

$$ \lambda_{\text{corner}} = \frac{\gamma^4 - \alpha^4}{\gamma^4 + \alpha^4}. \qquad (34b) $$

Here, the negative sign holds for the top right and bottom left corner sites on the sublattices $\tau = 1$ and 2, respectively, whereas the positive signs applies to the top left and bottom right corner sites on the sublattices $\tau = 3$ and 4, respectively.

For $\alpha = \gamma$, the PEPS of Eq. (29) describes the OAI phase of the quadrupole model. Indeed, in this case the edges have entanglement levels $\lambda_{\text{edge}} = 1/\sqrt{2}$ and the four corners have degenerate levels $\lambda_{\text{corner}} = 0$ corresponding to maximally entangled corner modes. On the other hand, for $\alpha = 0$, both $\lambda_{\text{edge}} = 1$ and $\lambda_{\text{corner}} = 1$ such that the TNS describes an atomic state.

Finally, if the OAI dimerized phase is perturbed by a small nonzero value $\beta \neq 0$, one can check numerically that the corner modes acquire a finite splitting $\lambda_{\text{corner}} \neq 0$ and are no longer perfectly localized at the corner sites. This confirms that the TNS with $\beta \neq 0$ is not in the OAI phase of the quadrupole model as expected from the breaking of the mirror symmetries and $C_4$ symmetry by the next-nearest-neighbor hopping and the chemical potential.

### C. 3D chiral hinge PEPS from dipole pumping

The PEPS of Eq. (29) with parameters $\alpha = \alpha_x = \alpha_y$ and $\beta = \beta_x = \beta_y$ describes a dimerized dipole-pumping cycle of the quadrupole model if $\alpha, \beta$, and $\gamma$ follow the parametrization $\phi_{\text{pump}}$ which we found for the dimerized charge-pumping cycle of the SSH model MPS. The coupling constants of the parent Hamiltonian $H_{\text{PEPS}}$ derived from the parametrization $\phi_{\text{pump}}$ are shown in Fig. 8(a). They differ from those of the dipole-pumping Hamiltonian of Ref. [3] only by a factor of $1/\sqrt{2}$ for the hopping amplitudes, which does not affect the topology of the interpolation. The single-particle ES of the PEPS along $\phi_{\text{pump}}$ is shown in Fig. 8(b). In the first half of the cycle, the system is in a dimerized phase with each unit cell decoupled from the rest of the system. Hence, the single-particle ES contains only the bulk bands with $\lambda = \pm 1$. However, in the second half of the cycle for $0 < t < \pi$, the contributions from the edges and corners can be clearly distinguished. The four corner modes connect the bands with $\lambda = \pm 1$ and are degenerate in the OAI dimerized phase of the quadrupole model obtained for $t = \pi/2$. In the three-dimensional second-order TI of Ref. [2] obtained by the identification of the time $t$ along the dipole-pumping cycle with the momentum $k_z$ in the third direction, these corner modes generate the chiral modes localized at the hinges.
there is an additional reduction of the number of entangled modes since the left and down virtual legs decouple due to the vanishing parameter $\beta = 0$ along the trivial interpolation (see Appendix G).

Analogously to Sec. III E, we therefore conclude that the bond dimension $D_{z,\tau}$ in the third direction grows exponentially with $N_z$ for $\tau = 1, 2, 3, 4$. Due to the mirror symmetries of the quadrupole model, $D_{z,\tau}$ has the same value as the vertical bond dimension of the two-dimensional PEPSs obtained from cyclic interpolations of the SSH model. The increase in spatial dimensionality therefore does not cause an increase of the bond dimension in the $(d+1)$st direction.

## V. CONCLUSION

In this paper, we showed how to use charge pumping to define TNSs for $(d+1)$-dimensional conventional or higher-order TIs starting from TNSs of TIs in $d$-space dimensions. To that end, we constructed a Gaussian fermionic MPS for the SSH model with bond dimension $D_{SSH} = 2$ in $d = 1$ dimension, and a Gaussian fermionic PEPS for the topological quadrupole model with bond dimension $D_{Quad} = 2$ in $d = 2$ dimensions. We proved that these TNSs have local gapped parent Hamiltonians with up to next-nearest-neighbor hopping, and thereby showed that they describe the SSH model along a charge-pumping cycle and the quadrupole model along a dipole-pumping cycle, respectively. We employed these TNSs to construct a two-dimensional PEPS for a Chern insulator and a three-dimensional PEPS for a chiral hinge higher-order topological insulator (HOTI), respectively. The $(d+1)$-dimensional TNSs inherit the finite bond dimension $D_{SSH}$ and $D_{Quad}$ in the first $d$ dimensions, respectively. In a hybrid coordinate system where the $(d+1)$st dimension corresponds to momentum, the $(d+1)$-dimensional TNSs have a finite bond dimension in this direction. In contrast, we showed that the bond dimension in the $(d+1)$st direction grows exponentially in a real-space coordinate system.

Our results suggest several directions for future work. On one hand, it would be interesting to study if a real-space PEPS for the Chern insulator with a polynomially growing bond dimension can be found by truncating the Schmidt values of the real-space column in our construction. Such a result could potentially provide insight into the physical origin for the obstructions preventing the existence of chiral PEPSs with a finite bond dimension. On the other hand, we expect the TNSs constructed here to be useful for finite-size simulations despite their growing bond dimension since the bond dimension is finite in all but one direction. By their local nature, they could be employed as the building block of interacting $(d+1)$-dimensional TIs obtained by Gutzwiller projection or parton constructions [58–61].

## ACKNOWLEDGMENTS

We thank I. Cirac, C. Hubig, A. Molnar, S. Pinon, A. Sterdyniak, and E. Zohar for discussions. A.H. and N.S. acknowledge support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through the ERC Starting Grant WASCOSYS

(Grant No. 636201), and by the Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy (Grant No. EXC-2111–390814868). A.H. and N.R. were supported by Grants No. ANR-17-CE30-0013-01 and No. ANR-16-CE30-0025. B.A.B. was supported by the Department of Energy Grant No. DE-SC0016239, the National Science Foundation EAGER Grant No. DMR 1643312, Simons Investigator Grant No. 404513, ONR Grant No.N00014-14-1-0330, the Packard Foundation, the Schmidt Fund for Innovative Research, ARO MURI Grant No. W911NF-12-1-0461 and by the NSF through the Princeton Center for Complex Materials, a Materials Research Science and Engineering Center DMR-1420541.

## APPENDIX A: COVARIANCE MATRICES OF GAUSSIAN FERMIONIC STATES

Eigenstates and thermal states of free-fermion systems are given by Gaussian states which satisfy Wick’s theorem. They are hence fully characterized by their covariance matrix (CM) of two-point correlation functions [51]. In this Appendix, we review the definitions of the complex and real CM for pure and mixed states in Appendix A 1. We summarize the relation with the entanglement spectrum in Appendix A 2 and provide the concrete expression for the CM of a Gaussian state parametrized as the exponential of a fermion bilinear in Appendix A 4.

### 1. Definitions

We consider a generic system of $N$ fermionic DOFs with annihilation operators $a_j$ and creation operators $a_j^\dagger$ for $j = 1, \dots, N$. In a pure or mixed state of this system, its covariance matrix is defined as

$$G_{\mu\nu} = \frac{i}{2} \langle [\chi_{\mu}^{\dagger}, \chi_{\nu}] \rangle = \begin{pmatrix} R^* & Q^* \\ Q & R \end{pmatrix}, \quad (A1)$$

where $\mu, \nu = 1, \dots, 2N$ and $\chi = (a_1, \dots, a_N, a_1^\dagger, \dots, a_N^\dagger)$ is the mode vector. The blocks $R$ and $Q$ of dimension $N \times N$ in Eq. (A1) are anti-Hermitian and antisymmetric, respectively, i.e., $R^\dagger = -R$ and $Q^T = -Q$, such that $G$ is anti-Hermitian. For a generic mixed state, $G$ satisfies the inequality $GG^\dagger \le \frac{1}{4}1$ and its eigenvalues come in complex-conjugate pairs $\pm \frac{i}{2}|\lambda_j|$ with $0 \le |\lambda_j| \le 1$ for $1 \le j \le N$. For pure states, $GG^\dagger = \frac{1}{4}1$ such that its eigenvalues are given by $|\lambda_j| = 1$.

For a state with charge conservation, the levels $-1 \le \lambda_j \le 1$ can be computed directly including their sign [62]: Since $\langle a_j^\dagger a_k^\dagger \rangle = 0$ for all $1 \le j, k \le N$, the off-diagonal block $Q^*$ vanishes in the complex CM $G$ of Eq. (A1), and the eigenvalues of the diagonal block $R^*$ are given by $\frac{i}{2}\lambda_j$.

For a generic system with pair creation and annihilation, it is convenient to define Majorana fermionic modes $c_{2j} = a_j^\dagger + a_j$, $c_{2j-1} = (-i)(a_j^\dagger - a_j)$ for $j = 1, \dots, N$ with $\{c_\mu, c_\nu\} = 2\delta_{\mu\nu}$ [63]. We denote the matrix corresponding to this basis change by

$$S_{\mu\nu} = \begin{cases} \delta_{\nu,\mu/2} + \delta_{\nu,\mu/2+N}, & \mu \text{ even} \\ i(\delta_{\nu,(\mu+1)/2} - \delta_{\nu,(\mu+1)/2+N}), & \mu \text{ odd} \end{cases} \quad (A2)$$
with $c_μ = \sum_ν S_{μν}χ_ν$ and $S^† = 2S^{-1}$. In terms of the Majorana modes, the CM

$$ \Gamma_{\mu\nu} = \frac{i}{2} \langle [c_{\mu}, c_{\nu}] \rangle = [S^* G S^T]_{\mu\nu} \quad (A3) $$

of size $2N \times 2N$ is real and antisymmetric. Hence, each of its singular values $|λ_j|$ with $1 \le j \le N$ is doubly degenerate. Moreover, $Γ$ satisfies $ΓΓ^† \le 1$ with equality for a pure state.

## 2. Relation to entanglement spectrum

We frequently consider the restriction of a pure quantum state $|ψ⟩$ to a subsystem $A$ of the entire system, which is described by the reduced density matrix $ρ_A = tr_{\bar{A}}[|ψ⟩⟨ψ|]$ obtained by tracing over the DOFs in the complement $\bar{A}$ of $A$. The many-body entanglement Hamiltonian $H_{Ent}$ with respect to this partition is given by the logarithm of the reduced density matrix as [23]

$$ ρ_A = \frac{1}{Z} e^{-H_{Ent}}, \quad (A4) $$

where $Z = tr(e^{-H_{Ent}})$.

If $|ψ⟩$ is a Gaussian state with CM $G$, then $ρ_A$ is Gaussian with its CM $G_A = (G_{μν})_{μ,ν∈A}$ given by the restriction of $G$ to the modes of $A$. In this case, the many-body entanglement Hamiltonian is a bilinear function of the fermionic mode operators defined by a square matrix referred to as the single-particle entanglement Hamiltonian. The eigenvalues $β_j$ of the single-particle entanglement Hamiltonian are related to the eigenvalues $±i/2|λ_j|$ of $G_A$ as $|λ_j| = tanh β_j/2$, where $1 \le j \le L$ and $L$ is the number of modes in the subsystem $A$ [52]. We therefore refer to the collection

$$ \{|λ_j|\}_{1≤j≤L} \quad (A5) $$

as the single-particle ES of $ρ_A$.

The single-particle ES can also be computed from the Majorana CM $Γ$ of Eq. (A3) for the Gaussian state $|ψ⟩$. One obtains the levels $|λ_j|$ for $1 \le j \le L$ directly as the singular values of the block $(Γ_{μν})_{μ,ν∈A}$ associated with the DOFs in $A$. On the other hand, the levels can be computed indirectly from the off-diagonal block $(Γ_{μν})_{μ∈A,ν∈\bar{A}}$ describing the correlations between $A$ and its complement $\bar{A}$. Indeed, the two blocks are coupled by the constraint $ΓΓ^T = 1$ originating in the purity of the state $|ψ⟩$. Concretely, the singular values of the off-diagonal block are given by

$$ μ_j = \sqrt{1 - λ_j^2} \quad (A6) $$

for $1 \le j \le L$ [64]. This allows to infer the single-particle entanglement energies for weakly entangled modes $|λ_j| \approx 1$ with an improved numerical accuracy.

## 3. Gaussian projections and Schur complements

In the construction of GfTNSs to be discussed in Appendix B, we will frequently encounter the following situation: Let us take a total system of $N + M$ fermionic modes with mode operators $a_j, a_j^†$ for $j = 1, ..., N + M$, and consider the subsystem defined by the $M$ last modes $j = N + 1, ..., N + M$. Let $|Q⟩$ be a Gaussian pure state of

the total system and $|ω⟩$ a Gaussian pure state of the subsystem of the last $M$ modes. Then, the projection

$$ |ψ⟩ = \langle ω|Q⟩, \quad (A7) $$

which is a state of only the first $N$ modes, is again a Gaussian state. The Majorana CMs $Γ_ψ, Γ_Q$, and $Γ_ω$ of the three Gaussian states are given by the expression in Eq. (A3), where the expectation values are taken in the states $|ψ⟩, |Q⟩$, and $|ω⟩$, respectively.

$Γ_ψ$ can be computed from $Γ_Q$ and $Γ_ω$ using a Schur complement coming from the Gaussian integration over the last $M$ modes [51,65–67]. Indeed, we can write the antisymmetric Majorana CM $Γ_Q$ of the unprojected state of the total system as

$$ Γ_Q = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ -\mathbf{B}^T & \mathbf{D} \end{pmatrix}, \quad (A8) $$

where the block $A$ of size $2N \times 2N$ refers to the Majorana modes corresponding to $a_j, a_j^†$ for $j = 1, ..., N$, and the block $D$ of size $2M \times 2M$ refers to the Majorana modes corresponding to $a_j, a_j^†$ for $j = N + 1, ..., N + M$. The block $B$ of size $2N \times 2M$ encodes the correlations between the first $N$ and last $M$ modes. In terms of these blocks and the CM $Γ_ω$ of size $2M \times 2M$, the CM of the projected state is

$$ Γ_ψ = \mathbf{A} + \mathbf{B}[\mathbf{D} + Γ_ω]^{-1}\mathbf{B}^T. \quad (A9) $$

Mathematically, this is the Schur complement of the block corresponding to the last $M$ modes.

## 4. Parametrization of Gaussian states

We consider a normalized Gaussian state with even fermionic parity which, anticipating Appendix B, we will denote by $|ψ⟩ = |Q⟩$. $|Q⟩$ is parametrized in terms of the antisymmetric complex matrix $M$ of dimension $N \times N$ as [51,68]

$$ |Q⟩ = \mathcal{N} \times e^{\sum_{ij} M_{ij} a_i^\dagger a_j^\dagger} |\Omega⟩, \quad (A10) $$

where $|Ω⟩$ is the fermionic vacuum and $\mathcal{N}$ is a normalization factor.

A relevant case is given by states whose parametrization matrices have nonzero entries only in the first row and column (or more generally, in the $n$th row and column). In particular, we will see in Appendices E and H that the local tensors of the MPS of Eq. (11) and the PEPS of Eq. (29) define Gaussian states of this form. In this case, $M_{ij}M_{kl}a_i^†a_j^†a_k^†a_l^† = 0$ for all $1 \le i, j, k, l \le N$ since all nonvanishing terms contain the factor $(a_1^†)^2 = 0$ [$(a_n^†)^2 = 0$ in the general case]. Therefore, the series expression of the exponential in Eq. (A10) terminates after first order such that

$$ |Q⟩ = \mathcal{N} \left[ 1 + \sum_{i,j=1}^{N} M_{ij} a_i^\dagger a_j^\dagger \right] |\Omega⟩, \quad (A11) $$

leading to the CM

$$ G = \frac{i}{2} \begin{pmatrix} -1_N + 8\mathcal{N}^2 \times M^\dagger M & -4\mathcal{N}^2 \times M^\dagger \\ -4\mathcal{N}^2 \times M & 1_N - 8\mathcal{N}^2 \times M M^\dagger \end{pmatrix}, \quad (A12) $$

where the normalization is $\mathcal{N}^{-2} = 1 + 2 \text{tr}(MM^\dagger)$.
APPENDIX B: CONSTRUCTION OF GfTNS

Gaussian fermionic tensor network states (GfTNS) describe free-fermion systems. They have the advantage that the contraction of the network can be performed in terms of covariance matrices, allowing efficient computations even for large systems. In this Appendix, we review the construction of Gaussian TNSs in one and two spatial dimensions. We begin in Appendix B1 by illustrating the construction of a TNS via fiducial states using the simple example of a bosonic MPS, where we do not need to take care of fermionic signs in tensor products. In Appendix B2, we move to the case of fermionic physical and virtual particles. In Appendix E, the concepts introduced here are illustrated using the pedagogical example of the SSH model MPS from Sec. II C.

1. Fiducial state approach

We recall that the construction of a TNS can be performed in several different but equivalent ways. Within the most well-known approach, the local information on the state is encoded in the local tensor with physical and virtual legs. For example, for a one-dimensional bosonic MPS $|\psi_{bMPS}\rangle$ the local tensor at position x takes the form $A[x]_{l_x, r_x}^{l_x}$, where $l_x$ is the physical index, $l_x$ is the left virtual index, and $r_x$ is the right virtual index. The global state is obtained by contracting, for each nearest-neighbor bond, the two virtual legs associated with this bond which belong to neighboring tensors. This is done by identifying and summing over the corresponding virtual indices. For example, the contraction of two neighboring MPS tensors gives $\sum_{r_x, l_x+1} A[x]_{l_x, r_x}^{l_x} \delta_{r_x, l_x+1} A[x+1]_{l_x+1, r_x+1}^{l_x+1}$.

More formally, this can be expressed as follows: First, we create a total maximally entangled state in the virtual layer whose role it is to implement the contraction of bonds. The total maximally entangled state is the tensor product over all nearest-neighbor bonds of a maximally entangled state of the two virtual particles for this bond. In our MPS example, the maximally entangled state for the bond between sites x and x + 1 is

$$|\omega_{x,x+1}\rangle = \sum_{r_x, l_x+1} \delta_{r_x, l_x+1} |r_x, l_x+1\rangle. \qquad (B1)$$

Second, at each site we translate the local tensor into a local projection map, which maps the virtual particles onto the physical particle at this site. The representation matrix of the local projection map is given by the local tensor. For instance, the MPS local projection map at site x is

$$\hat{A}[x] = \sum_{i_x, l_x, r_x} A[x]_{i_x, l_x}^{i_x} |i_x\rangle \langle l_x, r_x|. \qquad (B2)$$

The bond between sites x and x + 1 is contracted by applying the projection maps $\hat{A}[x]$ and $\hat{A}[x+1]$ to $|\omega_{x,x+1}\rangle$, giving

$$\begin{aligned} \hat{A}[x]\hat{A}[x+1] |\omega_{x,x+1}\rangle &= \sum_{\substack{i_x, i_{x+1} \\ l_x, r_x+1}} |i_x, i_{x+1}\rangle \langle l_x, r_{x+1}| \\ & \times \sum_{r_x, l_x+1} A[x]_{i_x, l_x}^{i_x} \delta_{r_x, l_x+1} A[x+1]_{l_x+1, r_x+1}^{l_x+1}. \end{aligned} \qquad (B3)$$

The global state is then obtained by contracting all bonds, i.e., by applying the product of all local projection maps to the total virtual maximally entangled state. For instance, the MPS on a chain with $N_x$ sites and periodic boundaries takes the form

$$|\psi_{bMPS}\rangle = \prod_{x=0}^{N_x-1} \hat{A}[x] \otimes \prod_{x=0}^{N_x-1} |\omega_{x,x+1}\rangle, \qquad (B4)$$

where the lattice site indices are $x = 0, \dots, N_x - 1$. This is a state of only the physical particles.

For GfTNSs, we will deal with fermionic particles and therefore follow a slightly different, but equivalent, approach to construct the TNSs using so-called fiducial states (see Refs. [47,68] for pedagogical introductions). In this approach, we consider fiducial states instead of local projection maps, which contain only creation operators. The fiducial state on each lattice site lies in the joint Hilbert space of the physical and virtual particles on this site. Its basis coefficients are given by the entries of the local tensor. For example, the local fiducial states of the MPS are

$$|Q_x\rangle = \sum_{i_x, l_x, r_x} A[x]_{i_x, l_x}^{i_x} |i_x\rangle |l_x, r_x\rangle. \qquad (B5)$$

Hence, the fiducial states are equivalent to the local projection maps or local tensors; in particular, they contain all the local information about the TNSs. The total state is obtained by projecting the tensor product of all local fiducial states on the total virtual maximally entangled state. For the MPS, we have

$$|\psi_{bMPS}\rangle = \bigotimes_{x=0}^{N_x-1} \langle\omega_{x,x+1}| \bigotimes_{x=0}^{N_x-1} |Q_x\rangle \qquad (B6)$$

which is equivalent to the expression in Eq. (B4). We are now ready to generalize the fiducial state formalism to fermionic physical and virtual particles.

2. Fermionic particles

We consider a one- or two-dimensional lattice system of free fermions with f fermionic modes per lattice site x, which are associated to the physical creation operators $a_{\tau, x}^{\dagger}$ with $\tau = 1, \dots, f$. A Gaussian fermionic TNS for this system with physical dimension $2^f$ and bond dimension $2^\xi$ is obtained by associating $\xi$ complex virtual fermionic modes with each physical lattice site x and nearest-neighbor direction $\alpha$, where $\alpha = L, R$ for one-dimensional MPSs with left and right nearest-neighbor bonds and $\alpha = L, U, R, D$ for two-dimensional PEPSs on the square lattice with left, up, right, and down nearest-neighbor bonds [33,68]. We denote the creation operators for these virtual modes by $b_{\alpha, j, x}^{\dagger}$ where $j = 1, \dots, \xi$ labels the different modes per bond and lattice site. For each lattice site, there are hence $n_{\text{modes}} = f + 2\xi$ modes for a fermionic MPS and $n_{\text{modes}} = f + 4\xi$ modes for a two-dimensional fermionic PEPS. The mode operators needed for the construction of the SSH MPS are discussed in the first paragraph of Appendix E.

We collect the mode operators associated with one lattice site into a mode vector

$$\chi_x = (a_{1,x}, a_{2,x}, \dots, a_{f,x}, b_{L,1,x}, b_{L,2,x}, \dots, b_{R,1,x}, \dots, \\ a_{1,x}^{\dagger}, \dots, b_{L,1,x}^{\dagger}, \dots, b_{R,1,x}^{\dagger}, \dots)^T \qquad (B7)$$
of length $2n_{\text{modes}}$. The physical and virtual mode operators obey canonical anticommutation relations

$$ \{ (\chi_x)_\mu, (\chi_x')_\nu \} = \delta_{x,x'} \delta_{|\mu-\nu|, n_{\text{modes}}} \quad (\text{B8}) $$

for $1 \le \mu, \nu \le 2n_{\text{modes}}$. Their joint vacuum state $|\Omega\rangle$ satisfies

$$ a_{\tau, x} |\Omega\rangle = b_{\alpha, j, x} |\Omega\rangle = 0 \quad (\text{B9}) $$

for all $\mathbf{x}, \tau = 1, \dots, f, \alpha$, and $j = 1, \dots, \xi$.

The local information about the TNS is contained in the local fiducial states $|Q_x\rangle$, introduced in Eq. (B5), which are equivalent to the local tensors $A[\mathbf{x}]$ used in the main text. We can easily translate between the two approaches since the basis elements of the local fiducial state are by definition equal to the local tensors [cf. Eq. (B5)]. For fermions, we write $|Q_x\rangle = Q_x|\Omega\rangle$ where $Q_x$ is a polynomial of creation operators which acts on the vacuum to create the fiducial state. For fermionic MPSs and PEPSs in one and two dimensions we have [49,57]

$$ A[\mathbf{x}]_{lr}^i = [\langle r | \otimes \langle i | \otimes \langle l | ] Q_x | \Omega \rangle, \quad (\text{B10}) $$

$$ A[\mathbf{x}]_{\text{lurd}}^{i} = [\langle d | \otimes \langle r | \otimes \langle u | \otimes \langle l | \otimes \langle i | ] Q_x | \Omega \rangle, \quad (\text{B11}) $$

where $|i\rangle$ with $i = 0, \dots, 2^f - 1$ is a basis for the Fock space associated with the physical mode operators $\{a_{\tau,x}^{\dagger}\}_{1\le\tau\le f}$, and $|l\rangle$ with $l = 0, \dots, 2^\xi - 1$ is a basis of the Fock space associated with the left virtual mode operators $\{b_{L,j,x}^{\dagger}\}_{1\le j\le \xi}$ on site $\mathbf{x}$ (similarly $|u\rangle$ for the up mode operators $\{b_{U,j,x}^{\dagger}\}_{1\le j\le \xi}$, $|r\rangle$ for the right mode operators $\{b_{R,j,x}^{\dagger}\}_{1\le j\le \xi}$, and $|d\rangle$ for the down mode operators $\{b_{D,j,x}^{\dagger}\}_{1\le j\le \xi}$). For the SSH model MPS, the fiducial states obtained thus are given in Eq. (E1).

A Gaussian fermionic TNS has the property that all local fiducial states $|Q_x\rangle$ satisfy Wick's theorem. In this case, the global physical state $|\psi\rangle$ is also Gaussian [33,51]. We denote by $\Gamma_Q$ the Majorana CM of the product of all fiducial states, $\prod_x Q_x|\Omega\rangle$, also referred to as total fiducial state. From now onward, we consider GfTNSs with parity-even local tensors as discussed in Sec. II B, whose local fiducial states therefore have an even number of physical and virtual fermions. The maps $Q_x$ can thus be expressed as in Eq. (A10) as the exponential of a quadratic form of the physical and virtual creation operators on the site $\mathbf{x}$, which are contained in the last $n_{\text{modes}}$ entries $[(\chi_x)_{n_{\text{modes}}+m}]_{1\le m\le n_{\text{modes}}}$ of the mode vector of Eq. (B7) [51,68]. Concretely, the local fiducial state is parametrized as

$$ Q_x = \exp \left[ \sum_{m, m'=1}^{n_{\text{modes}}} (M_x)_{mm'} (\chi_x)_{n_{\text{modes}}+m} (\chi_x)_{n_{\text{modes}}+m'} \right], \quad (\text{B12}) $$

with an antisymmetric square matrix $M_x$ of dimension $n_{\text{modes}}$. For the SSH model MPS, this matrix is given in Eq. (E2).

In order to illustrate the concepts introduced above, let us construct a Gaussian maximally entangled state of the virtual fermions for each nearest-neighbor bond $\langle\mathbf{x}'\mathbf{x}''\rangle$. Let us denote by $\alpha'$ and $\alpha''$ the type of virtual fermion involved in the bond $\langle\mathbf{x}'\mathbf{x}''\rangle$ on the site $\mathbf{x}'$ and $\mathbf{x}''$, respectively. For instance, if $\mathbf{x}'' = \mathbf{x}' + \hat{\mathbf{x}}$, then $\alpha' = R$ and $\alpha'' = L$. A fermionic maximally entangled state for this bond is then given by [33]

$$ |\omega_{\mathbf{x}'\mathbf{x}''}\rangle = \prod_{j=1}^{\xi} \frac{1}{\sqrt{2}} (1 + b_{\alpha', j, \mathbf{x}'}^{\dagger} b_{\alpha'', j, \mathbf{x}''}^{\dagger}) |\Omega\rangle. \quad (\text{B13}) $$

Due to the fermionic anticommutation relations, this expression is not symmetric under exchange of $(\alpha', \mathbf{x}')$ and $(\alpha'', \mathbf{x}'')$; we say that the bond points from the initial site $\mathbf{x}'$ to the final site $\mathbf{x}''$.

The state $|\omega_{\mathbf{x}'\mathbf{x}''}\rangle$ is a Gaussian state satisfying Wick's theorem. To see this, we first consider the simple case of a one-dimensional MPS with $\xi = 1$ virtual fermion per nearest-neighbor bond and lattice site. In this case, the virtual maximally entangled state from Eq. (B13) for the bond $\langle\mathbf{x}, \mathbf{x} + \hat{\mathbf{x}}\rangle$ becomes

$$ \begin{aligned} |\omega_{\mathbf{x}, \mathbf{x}+\hat{\mathbf{x}}}\rangle &= \frac{1}{\sqrt{2}}[1 + b_{R, \mathbf{x}}^{\dagger} b_{L, \mathbf{x}+\hat{\mathbf{x}}}^{\dagger}]|\Omega\rangle \\ &= \frac{1}{\sqrt{2}}\left[1 + (b_{R, \mathbf{x}}^{\dagger} b_{L, \mathbf{x}+\hat{\mathbf{x}}}) \frac{1}{2} \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} b_{R, \mathbf{x}}^{\dagger} \\ b_{L, \mathbf{x}+\hat{\mathbf{x}}}^{\dagger} \end{pmatrix}\right]|\Omega\rangle, \end{aligned} \quad (\text{B14}) $$

where the bond points from site $\mathbf{x}$ to site $\mathbf{x} + \hat{\mathbf{x}}$. This is a Gaussian state of the form of Eq. (A11) parametrized by the antisymmetric matrix $M = (i/2)\sigma_2$. Hence, its complex CM $G$ can be computed from Eq. (A12) and is given by

$$ G = -\frac{1}{2} \begin{pmatrix} 0 & \sigma_2 \\ -\sigma_2 & 0 \end{pmatrix}, \quad (\text{B15}) $$

where $\sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}$ denotes the second Pauli matrix. According to Eq. (A3), the corresponding Majorana CM $\Gamma$ is obtained by conjugation with the matrix $S$, and is found to be

$$ \Gamma = \begin{pmatrix} 0 & \sigma_1 \\ -\sigma_1 & 0 \end{pmatrix}, \quad (\text{B16}) $$

where $\sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$ denotes the first Pauli matrix. For a two-dimensional TNS with $\xi = 1$, the real CM of the virtual state for the vertical nearest-neighbor bond $\langle\mathbf{x}, \mathbf{x} - \hat{\mathbf{y}}\rangle$ is also given by Eq. (B16). In this case, we choose the bond to be oriented downward from $\mathbf{x}$ to $\mathbf{x} - \hat{\mathbf{y}}$ to comply with our ordering of the virtual legs as $L, U, R, D$. For TNS with more virtual fermions, $\xi > 1$, the virtual maximally entangled state from Eq. (B13) is a tensor product of multiple states of the form of Eq. (B14). Hence, the corresponding real CM is a direct sum of multiple copies of the $\Gamma$ from Eq. (B16). We denote by $\Gamma_\omega$ the Majorana CM of the total virtual maximally entangled state $\otimes_{\langle\mathbf{x}'\mathbf{x}''\rangle} |\omega_{\mathbf{x}'\mathbf{x}''}\rangle$. Since the total virtual maximally entangled state is a tensor product over all bonds, $\Gamma_\omega$ is a direct sum of multiple copies of the $\Gamma$ from Eq. (B16).

As explained in Appendix B1, the global physical state $|\psi\rangle$ is obtained from the constituents introduced above by projecting the total fiducial state onto the virtual maximally entangled state

$$ |\psi\rangle = \left[ \bigotimes_{\langle\mathbf{x}'\mathbf{x}''\rangle} \langle\omega_{\mathbf{x}'\mathbf{x}''}| \right] \prod_{\mathbf{x}} Q_{\mathbf{x}} |\Omega\rangle. \quad (\text{B17}) $$

We recall from Appendix A3 that for Gaussian states, this projection can be formulated in terms of CMs and gives rise to a Schur complement. This is the approach we take in the following.

Let us see how we can apply the general Schur complement expression of Eq. (A9) in order to compute the Majorana
CM $\Gamma_{|\psi\rangle}$ for the physical state from the CMs $\Gamma_Q$ and $\Gamma_\omega$ defined above for the fiducial and virtual maximally entangled states. We introduce the symbols $p$ and $v$ to collectively refer to all Majorana mode operators for the physical and virtual fermions. Therefore, $(\Gamma_Q)_{pp}$ and $(\Gamma_Q)_{vv}$ denote the blocks of the CM of the total fiducial state that describe the reduced state of only the physical and virtual degrees of freedom (DOFs), respectively. On the other hand, the blocks $(\Gamma_Q)_{pv} = -(\Gamma_Q)_{vp}^T$ encode the correlations between physical and virtual fermions. In Eq. (B17), we are projecting the fiducial state of physical and virtual modes onto a maximally entangled state of the virtual modes in order to obtain a state of only the physical modes. According to Eq. (A9), this is represented by the Schur complement of the virtual block $(vv)$ given by

$$ \Gamma_{|\psi\rangle} = (\Gamma_Q)_{pp} + (\Gamma_Q)_{pv}[(\Gamma_Q)_{vv} + \Gamma_\omega]^{-1}(\Gamma_Q)_{pv}^T. \quad (B18) $$

Note that the CM on the left-hand side of this equation tracks every physical mode independently. Hence, its dimension is proportional to the system size and becomes large for our cases of interest. We will now specialize Eq. (B18) to translation-invariant states, where we can achieve a massive reduction of the size of the physical CM down to the number of Bloch bands independent of the system size.

## APPENDIX C: TRANSLATION-INVARIANT GfTNSs

In this Appendix, we specialize the formalism reviewed in Appendix B to translation-invariant GfTNSs, where the contraction of the network can often be performed analytically. In Appendix C1, we compute the CM of a translation-invariant GfTNS, which we then use in Appendix C2 to construct a parent Hamiltonian for the state.

### 1. Covariance matrix

For translation-invariant Gaussian TNSs, we introduce the Fourier transform (FT) of the physical and virtual mode operators as

$$ a_{\tau, \mathbf{k}} = \sum_{\mathbf{x}} \mathcal{F}_{\mathbf{k}, \mathbf{x}} a_{\tau, \mathbf{x}}, \quad (C1a) $$

$$ b_{\alpha, j, \mathbf{k}} = \sum_{\mathbf{x}} \mathcal{F}_{\mathbf{k}, \mathbf{x}} b_{\alpha, j, \mathbf{x}} \quad (C1b) $$

for all $\tau = 1, \dots, f$, $j = 1, \dots, \xi$, $\alpha = L, R$ for MPSs, and $\alpha = L, U, R, D$ for two-dimensional PEPs. Here, the FT in two spatial dimensions with $N_x$ and $N_y$ sites in the horizontal and vertical direction, respectively, position vector $\mathbf{x} = (x, y)$ and momentum vector $\mathbf{k} = (k_x, k_y)$, is given by

$$ \mathcal{F}_{\mathbf{k}, \mathbf{x}} = \frac{1}{\sqrt{N_x N_y}} e^{-i \mathbf{k} \cdot \mathbf{x}}. \quad (C1c) $$

The momenta in the horizontal ($k_x$) and vertical ($k_y$) direction take values $k_x = \frac{2\pi j}{N_x}$ with $0 \le j \le N_x - 1$ and $k_y = \frac{2\pi j}{N_y}$ with $0 \le j \le N_y - 1$. The FT for a single spatial direction is analogous.

We define the FT of the mode vector $\chi_x$ from Eq. (B7) to be

$$ \chi_{\mathbf{k}} = \sum_{\mathbf{x}} \mathcal{F}_{\mathbf{k}, \mathbf{x}} \chi_{\mathbf{x}}. \quad (C2) $$

Therefore,

$$ \chi_{\mathbf{k}} = (a_{1,\mathbf{k}}, a_{2,\mathbf{k}}, \dots, a_{f,\mathbf{k}}, b_{L,1,\mathbf{k}}, b_{L,2,\mathbf{k}}, \dots, b_{R,1,\mathbf{k}}, \dots, a_{1,-\mathbf{k}}^{\dagger}, \dots, b_{L,1,-\mathbf{k}}^{\dagger}, \dots, b_{R,1,-\mathbf{k}}^{\dagger}, \dots)^T \quad (C3) $$

mixes mode operators at momenta $\mathbf{k}$ and $-\mathbf{k}$ similarly to a Nambu spinor.

Due to the translation invariance of the GfTNS, the FT brings the CMs of the physical state, the total fiducial state, and the total virtual maximally entangled state into a block-diagonal form. We denote the Majorana CM of the total fiducial state with respect to the Fourier transform of the mode operators by

$$ (\Gamma_Q)_{\mu, \nu; \mathbf{k}\mathbf{q}} = \sum_{\mathbf{x}, \mathbf{x}'} \mathcal{F}_{\mathbf{k}, \mathbf{x}}^*(\Gamma_Q)_{\mu, \nu; \mathbf{x}\mathbf{x}'} \mathcal{F}_{\mathbf{q}, \mathbf{x}'} = \delta_{\mathbf{k}\mathbf{q}}(\tilde{\Gamma}_Q(\mathbf{k}))_{\mu, \nu}, \quad (C4) $$

where the last equality defines the Majorana CM $\tilde{\Gamma}_Q(\mathbf{k})$ restricted to the block of momentum $\mathbf{k}$. Note that the size of $\tilde{\Gamma}_Q(\mathbf{k})$ is given by twice the number of Bloch bands, $2 \times n_{\text{modes}}$, and therefore no longer grows with the system size. Analogous statements hold for the physical CM $\Gamma_{|\psi\rangle}$ with the Fourier block matrix $\tilde{\Gamma}_{|\psi\rangle}(\mathbf{k})$ of size $2f$, and for the CM of the total virtual maximally entangled state $\Gamma_\omega$ with Fourier block $\tilde{\Gamma}_\omega(\mathbf{k})$ of size $2(n_{\text{modes}} - f)$.

Since the TNS is translation invariant, the local fiducial state $|Q_\chi\rangle$ is the same on every unit cell. Hence, the Fourier CM of the total fiducial state is localized at momentum $\mathbf{k} = 0$,

$$ \tilde{\Gamma}_Q(\mathbf{k}) = \Gamma_{Q_x} \times \delta_{\mathbf{k},0}. \quad (C5) $$

Here, $\Gamma_{Q_x}$ is the CM of the local fiducial state $|Q_\chi\rangle$ on a single site of dimension $2n_{\text{modes}}$. It has the block structure

$$ \Gamma_{Q_x} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ -\mathbf{B}^T & \mathbf{D} \end{pmatrix}, \quad (C6) $$

where the real antisymmetric blocks $\mathbf{A}$ and $\mathbf{D}$ of dimension $2f$ and $2(n_{\text{modes}} - f)$ describe the physical and virtual subspaces, respectively, whereas the block $\mathbf{B}$ encodes the coupling between physical and virtual modes.

On the other hand, the CM of the total virtual maximally entangled state $\tilde{\Gamma}_\omega(\mathbf{k})$ has a nontrivial momentum dependence since the maximally entangled states from Eq. (B13) connect different unit cells. It is a direct sum of the contributions from the different spatial directions. For a single virtual fermion with $\xi = 1$, the Majorana Fourier CM of the horizontal bonds oriented from left to right is the Fourier transform of the matrix $\Gamma$ from Eq. (B16). It reads as [40]

$$ \tilde{\Gamma}_{\omega}(k_{\chi}) = \begin{pmatrix} 0 & -\sigma_{1}e^{-ik_{\chi}} \\ \sigma_{1}e^{ik_{\chi}} & 0 \end{pmatrix}. \quad (C7) $$

This matrix is written in the basis of the FT of the Majorana operators constructed from the complex modes $(b_{L,x}, b_{R,x}, b_{L,x}^{\dagger}, b_{R,x}^{\dagger})$ where we have omitted the index $j$ since $\xi = 1$. The contribution from the vertical bonds is given by Eq. (C7) with $k_x \mapsto -k_y$ since the vertical bonds are oriented downward and hence in the direction of negative $k_y$. Therefore, the Fourier CM of the virtual bonds is anti-Hermitian and satisfies the identities $\tilde{\Gamma}_\omega(\mathbf{k})^* = \tilde{\Gamma}_\omega(-\mathbf{k})$ and $\tilde{\Gamma}_\omega(\mathbf{k})^T = -\tilde{\Gamma}_\omega(-\mathbf{k})$.
We are now in a position to compute the physical Majorana CM in momentum space by taking the FT of Eq. (B18). In momentum space, the right-hand side simplifies significantly due to the expression for the fiducial state CM given in Eqs. (C5) and (C6). We thus find [33]

$$ \tilde{\Gamma}_{|\psi\rangle}(\mathbf{k}) = \mathbf{A} + \mathbf{B}[\mathbf{D} + \tilde{\Gamma}_{\omega}(\mathbf{k})]^{-1}\mathbf{B}^T. \quad (C8) $$

The physical state is well defined if the matrix inversion can be carried out, i.e., if the determinant

$$ q(\mathbf{k}) = \det[\mathbf{D} + \tilde{\Gamma}_{\omega}(\mathbf{k})] \quad (C9) $$

is not equal to zero. We emphasize that the matrices in Eq. (C8) have a constant size given by the number of Bloch bands, such that Eq. (C8) can typically be evaluated analytically. In Appendix E, we use Eq. (C8) to evaluate the Bloch CM of the SSH model MPS from Sec. II A, and in Appendix H, we use Eq. (C8) to evaluate the Bloch CM of the quadrupole model PEPS of Sec. IV B.

## 2. Parent Hamiltonian

A translation-invariant Gaussian fermionic TNS $|\psi\rangle$ has an infinite number of parent Hamiltonians for which it is an exact ground state. For any non-negative scalar function $\epsilon(\mathbf{k}) \ge 0$ on the Brillouin zone,

$$ H_{\epsilon} = \frac{i}{4} \sum_{\mathbf{k}} \epsilon(\mathbf{k}) (\chi_{\mathbf{k}}^{(p)})^{\dagger} [S^T \tilde{\Gamma}_{|\psi\rangle}(\mathbf{k}) S^{*}] \chi_{\mathbf{k}}^{(p)} \quad (C10) $$

is a parent Hamiltonian for the Gaussian fermionic TNS $|\psi\rangle$ [33,40]. Here,

$$ \chi_{\mathbf{k}}^{(p)} = (a_{1,\mathbf{k}}, \dots, a_{f,\mathbf{k}}, a_{1,-\mathbf{k}}^{\dagger}, \dots, a_{f,-\mathbf{k}}^{\dagger})^{\mathrm{T}} \quad (C11) $$

is the physical part of the mode vector of Eq. (C3), and S is the transformation matrix defined in Eq. (A2) from the Majorana fermions to the complex fermions.

The properties of the parent Hamiltonian $H_{\epsilon}$ depend on the dispersion function $\epsilon(\mathbf{k})$. If $\epsilon(\mathbf{k}) > 0$ is strictly positive throughout the Brillouin zone, $H_{\epsilon}$ is gapped. Moreover, if all matrix entries of the product $\epsilon(\mathbf{k})\tilde{\Gamma}_{|\psi\rangle}(\mathbf{k})$ are polynomials in $e^{\pm ik_{\tau}}$ and $e^{\pm ik_{\tau}}$, the parent Hamiltonian is strictly local.

A natural, but not unique, choice for the dispersion function is given by $\epsilon(\mathbf{k}) = q(\mathbf{k})$ from Eq. (C9). Indeed, since D and $\tilde{\Gamma}_{\omega}(\mathbf{k})$ are anti-Hermitian and of even dimension [see Eqs. (C6) and (C7)], it follows that $q(\mathbf{k})$ is real. If moreover $q(\mathbf{k})$ is strictly positive throughout the Brillouin zone, implying that the PEPS has exponentially decaying real-space correlations, the parent Hamiltonian $H_q$ is gapped and strictly local with all terms acting on at most $2\xi$ successive unit cells [40,47].

## APPENDIX D: GfTNSs WITH CONSERVED PARTICLE NUMBER

In the main text, we consider Gaussian fermionic TNSs with a conserved particle number: the ground states of both the SSH model and the quadrupole model lie at half-filling. The TNSs are written in a basis related to the physical basis by a staggered particle-hole conjugation [cf. Eq. (7) for the SSH MPS and Eq. (28) for the quadrupole PEPS]. Hence, the U(1) symmetry of the local tensors, which imposes the

conservation of the physical particle number, also takes a staggered form [cf. Eq. (12) for the SSH MPS and Eq. (30) for the quadrupole PEPS]. In this Appendix, we rephrase this U(1) symmetry in the language of fiducial states, and show that it enforces many vanishing elements for the CM of the local fiducial state. These are equivalent to the vanishing of the off-diagonal block $Q = 0$ of the complex CM of a state with conserved particle number (see Appendix A 1), but expressed in the basis after the staggered particle-hole transformation. We focus on the one-dimensional case since the computation for two-dimensional GfTNSs is analogous.

For one-dimensional Gaussian fermionic MPS, we consider a U(1) symmetry of the local tensor of the general form

$$ A[\mathbf{x}]_{lr}^{i} = \sum_{i' l' r'} \left( \bigotimes_{\tau=0}^{f} U(\eta_{\tau, x} \varphi) \right) \left( \bigotimes_{j=0}^{\xi} U(\eta_{L, j, x} \varphi) \right)_{ll'} \\ \times \left( \bigotimes_{j=0}^{\xi} U(\eta_{R, j, x} \varphi) \right)_{rr'}, \quad (D1) $$

where $\eta_{\tau, x}$, $\eta_{L, j, x}$, $\eta_{R, j, x} \in \{\pm 1\}$ for $\tau = 1, \dots, f$ and $j = 1, \dots, \xi$. Here,

$$ U(\varphi) = \begin{pmatrix} 1 & 0 \\ 0 & e^{i\varphi} \end{pmatrix} \quad (D2) $$

is the U(1) rotation acting on a single spinless fermion. Positive and negative values for $\eta$ indicate that the corresponding physical or virtual modes transform as particles and holes, respectively. The symmetries of Eq. (12) for the local tensors on the A and B sublattices of the SSH charge-pumping MPS are examples with $f = \xi = 1$.

Since the elements of the local tensor $A[\mathbf{x}]$ are the basis coefficients of the local fiducial state $|Q_{\mathbf{x}}\rangle$ [cf. Eq. (B10)], Eq. (D1) is equivalent to the invariance of the local fiducial state under the U(1) symmetry

$$ \hat{U}_{\mathbf{x}}(\varphi) = \prod_{\tau=0}^{f} \hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}} \varphi) \prod_{j=0}^{\xi} \hat{U}_{L, j, \mathbf{x}}(\eta_{L, j, \mathbf{x}} \varphi) \hat{U}_{R, j, \mathbf{x}}(\eta_{R, j, \mathbf{x}} \varphi) \quad (D3) $$

whose many-body basis representation is given in Eq. (D1). Here, each individual operator $\hat{U}$ in the product acts on exactly one fermion. For instance, the operator acting on the physical fermion $\tau$ is given by

$$ \hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}} \varphi) = e^{i \eta_{\tau, \mathbf{x}} \varphi} a_{\tau, \mathbf{x}}^{\dagger} a_{\tau, \mathbf{x}} + a_{\tau, \mathbf{x}} a_{\tau, \mathbf{x}}^{\dagger} \quad (D4) $$

and similarly for the virtual fermions. This agrees with the matrix representation of the U(1) rotation acting on a single mode given in Eq. (D2).

We observe that the U(1) operator of the physical fermion $\tau$ from Eq. (D4) satisfies $a_{\tau, \mathbf{x}}\hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}}\varphi) = e^{i\eta_{\tau, \mathbf{x}}\varphi}\hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}}\varphi)a_{\tau, \mathbf{x}}$ and $a_{\tau, \mathbf{x}}^{\dagger}\hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}}\varphi) = e^{-i\eta_{\tau, \mathbf{x}}\varphi}\hat{U}_{\tau, \mathbf{x}}(\eta_{\tau, \mathbf{x}}\varphi)a_{\tau, \mathbf{x}}^{\dagger}$. Extending this to all of the modes in the mode operator of Eq. (B7), we find that

$$ (\chi_{\mathbf{x}})_{\mu} \hat{U}_{\mathbf{x}}(\varphi) = e^{i\varphi(\eta_{\mathbf{x}})_{\mu}} \hat{U}_{\mathbf{x}}(\varphi)(\chi_{\mathbf{x}})_{\mu}, \quad (D5) $$

where $1 \le \mu, \nu \le 2n_{\text{modes}}$. Here, we collected all the parameters $\eta$ into a vector of length $2n_{\text{modes}}$ using the same ordering
as in the mode vector of Eq. (B7),

$$ \eta_x = (\eta_{1,x}, \eta_{2,x}, \dots, \eta_{f,x}, \eta_{L,1,x}, \eta_{L,2,x}, \dots, \eta_{R,1,x}, \dots, \\ - \eta_{1,x}, \dots, -\eta_{L,1,x}, \dots, -\eta_{R,1,x}, \dots), \quad (D6) $$

where the last $n_{\text{modes}}$ entries describe the creation operators which transform with the opposite sign under the U(1) symmetry.

We are now ready to infer the consequences of the local fiducial state's invariance under the U(1) symmetry from Eq. (D3) for its complex CM. Indeed, due to the relation of Eq. (D5), correlation functions transform under the symmetry as

$$ \begin{aligned} \langle Q_x | (\chi_x^\dagger)_\nu (\chi_x)_\mu | Q_x \rangle &= \langle Q_x | \hat{U}_x^\dagger(\varphi) (\chi_x^\dagger)_\nu (\chi_x)_\mu \hat{U}_x(\varphi) | Q_x \rangle \\ &= e^{i\varphi[(\eta_x)_\mu - (\eta_x)_\nu]} \langle Q_x | (\chi_x^\dagger)_\nu (\chi_x)_\mu | Q_x \rangle. \end{aligned} \quad (D7) $$

Therefore, they vanish unless $(\eta_x)_\mu = (\eta_x)_\nu$. This shows that the symmetry of Eq. (D1) forces the vanishing of half the elements of the complex CM of the local fiducial state

$$ (G_x)_{\mu\nu} = 0 \quad \text{if } (\eta_x)_\mu \neq (\eta_x)_\nu, \quad (D8) $$

where $1 \le \mu, \nu \le 2n_{\text{modes}}$.

## APPENDIX E: COVARIANCE MATRIX FOR SSH MPS

In this Appendix, we illustrate the formalism of GfTNS introduced in Appendices B and C using the example of the MPS from Eq. (11) describing charge pumping in the SSH model. After expressing the state as a GfTNS in Appendix E 1, we demonstrate the computation of its Bloch CM and parent Hamiltonian in Appendix E 2.

### 1. Expression as GfTNS

Our goal is to express the SSH charge-pumping MPS, which is already fully defined by Eq. (11) in the language of local tensors, as a GfTNS using the formalism from Appendix B. Since the MPS has physical dimension 2 and bond dimension 2, it is described by $f=1$ physical fermion and $\xi=1$ virtual fermion per nearest-neighbor bond and lattice site (cf. the first paragraph of Appendix B 2). On the A sublattice, the annihilation operators for the physical, left virtual, and right virtual modes are $a_{A,x}$, $b_{L,A,x}$, and $b_{R,A,x}$, respectively, and similarly for the B sublattice. Following the recipe given above, we need to find the local fiducial states $|Q_{A,x}\rangle = Q_{A,x}|\Omega\rangle$ and $|Q_{B,x}\rangle = Q_{B,x}|\Omega\rangle$ that match the local tensors A and B from Eq. (11) at each unit cell x. The link between the local tensors and fiducial states is then given by Eq. (B10) stating that the local tensor entries are the basis coefficients of the fiducial state. For example, the tensor element $A_{01}^1 = \beta$ tells us that the local fiducial map $Q_{A,x}$ contains a term $\beta a_{A,x}^\dagger b_{R,A,x}^\dagger$. Performing this matching for every nonzero tensor entry, we find that the fiducial maps are given by

$$ Q_{A,x} = \gamma + \alpha a_{A,x}^{\dagger} b_{L,A,x}^{\dagger} + \beta a_{A,x}^{\dagger} b_{R,A,x}^{\dagger}, \quad (E1a) $$

$$ Q_{B,x} = \gamma + \beta a_{B,x}^{\dagger} b_{L,B,x}^{\dagger} - \alpha a_{B,x}^{\dagger} b_{R,B,x}^{\dagger}. \quad (E1b) $$

In a second step, we want to write the local fiducial states as Gaussian states satisfying Wick's theorem in order to express the MPS as a GfTNS. Since the local tensors are parity even, the local fiducial states $|Q_{A,x}\rangle$ and $|Q_{B,x}\rangle$ can be parametrized as in Eq. (B12) using the exponential of antisymmetric coefficient matrices $M_{A,x}$ and $M_{B,x}$. In the case of the SSH pumping MPS, this is very simple since the fiducial states we derived in Eq. (E1) have the form of Eq. (A11), allowing us to directly read off $M_{A,x}$ and $M_{B,x}$ (imposing antisymmetry). We find that the coefficient matrices are given by

$$ M_{A,x} = \frac{1}{2} \begin{pmatrix} 0 & a & b \\ -a & 0 & 0 \\ -b & 0 & 0 \end{pmatrix}, \quad (E2a) $$

$$ M_{B,x} = \frac{1}{2} \begin{pmatrix} 0 & b & -a \\ -b & 0 & 0 \\ a & 0 & 0 \end{pmatrix}. \quad (E2b) $$

Here, we defined the quotients $a = \alpha/\gamma$ and $b = \beta/\gamma$ of the parameters used in Eq. (11), and we absorbed the remaining factor $\gamma$ into the normalization constant $N$ in Eq. (A11). The complex CMs of these local fiducial states are computed from $M_{A,x}$ and $M_{B,x}$ using Eq. (A12), and are then transformed to the Majorana representation using Eq. (A3). One finds that the Majorana CM $\Gamma_{Q_{A,x}}$ is

$$ \Gamma_{Q_{A,x}} = \frac{1}{c} \begin{pmatrix} (c-2)i\sigma_2 & 4a\sigma_1 & 4b\sigma_1 \\ -4a\sigma_1 & (8a^2-c)i\sigma_2 & 8abi\sigma_2 \\ -4b\sigma_1 & 8abi\sigma_2 & (8b^2-c)i\sigma_2 \end{pmatrix}, \quad (E3) $$

where we introduced the shorthand notation $c \equiv 1 + 4a^2 + 4b^2$, and the basis is the Majorana basis derived from $(a_{A,x}, b_{L,A,x}, b_{R,A,x}, a_{A,x}^\dagger, b_{L,A,x}^\dagger, b_{R,A,x}^\dagger)$. The Majorana CM $\Gamma_{Q_{B,x}}$ in the Majorana basis derived from $(a_{B,x}, b_{L,B,x}, b_{R,B,x}, a_{B,x}^\dagger, b_{L,B,x}^\dagger, b_{R,B,x}^\dagger)$ is given by Eq. (E3) with the replacements $a \mapsto b$ and $b \mapsto -a$. We have thus successfully expressed the MPS from Eq. (11) as a GfTNS.

The virtual maximally entangled states for the SSH MPS are of the form discussed in Appendix B 2. In particular, the Majorana CM $\Gamma_{|\omega_{A,x,B,x}\rangle}$ of the state $|\omega_{A,x,B,x}\rangle$ within a unit cell is given by Eq. (B16) in the Majorana basis obtained from $(b_{R,A,x}, b_{L,B,x}, b_{R,A,x}^\dagger, b_{L,B,x}^\dagger)$. Similarly, the Majorana CM of the state $|\omega_{B,x,A,x+1}\rangle$ between unit cells is given by Eq. (B16) in the Majorana basis obtained from $(b_{R,B,x}, b_{L,A,x+1}, b_{R,B,x}^\dagger, b_{L,A,x+1}^\dagger)$.

### 2. Bloch CM and parent Hamiltonian

Now that we have expressed the MPS from Eq. (11) as a GfTNS, we want to use the power of the formalism introduced in Appendix C to derive its CM and parent Hamiltonian on a chain with periodic boundary conditions. We will compute the Bloch CM by evaluating Eq. (C8), and from there obtain the parent Hamiltonian via Eq. (C10).

#### a. CM for a unit cell

The expression for the Bloch CM from Eq. (C8) is valid for a translation-invariant GfTNS. In particular, the CM $\Gamma_{Q_x}$ with its blocks A, B, and D from Eq. (C5) is the Majorana CM of the fiducial state of a unit cell, not a single site. In
order to proceed, we therefore need to derive the fiducial state of an entire unit cell by contracting the virtual bond within a unit cell. In the language of GfTNSs, this is described by projecting the fiducial states $Q_{A,x}Q_{B,x}|\Omega\rangle$ of the two sites in one unit cell onto the virtual maximally entangled state connecting them:

$$ \langle\omega_{AB}|[Q_{A,x}Q_{B,x}|\Omega]\rangle. \qquad (E4) $$

This projection is a special case of Eq. (A7). Hence, the CM $\Gamma_{Q_x}$ of the resulting state is given, as in Eq. (A9), by the Schur complement of the block of size $4 \times 4$ corresponding to the Majorana mode operators constructed from $(b_{R,A}, b_{L,B}, b_{R,A}^{\dagger}, b_{L,B}^{\dagger})$:

$$ \Gamma_{Q_x} = \mathbf{A}^{\prime} + \mathbf{B}^{\prime} [\mathbf{D}^{\prime} + \Gamma_{[\omega_{A,x,B,x}]}]^{-1} \mathbf{B}^{\prime T}. \qquad (E5) $$

Here, $\mathbf{A}'$ and $\mathbf{D}'$ are the blocks on the diagonal of the direct sum $\Gamma_{Q_{A,x}} \oplus \Gamma_{Q_{B,x}}$ corresponding to the Majorana modes derived from $(a_{A,x}, a_{B,x}, b_{L,A,x}, b_{R,B,x}, a_{A,x}^{\dagger}, a_{B,x}^{\dagger}, b_{L,A,x}^{\dagger}, b_{R,B,x}^{\dagger})$ and $(b_{R,A,x}, b_{L,B,x}, b_{R,A,x}^{\dagger}, b_{L,B,x}^{\dagger})$, respectively. Correspondingly, $\mathbf{B}'$ is the off-diagonal block.

$\Gamma_{Q_x}$ has the block form

$$ \Gamma_{Q_x} = \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ -\mathbf{B}^T & \mathbf{D} \end{pmatrix}, \qquad (E6) $$

where $\mathbf{A}$ and $\mathbf{D}$ are real antisymmetric blocks of size $4 \times 4$, and $\mathbf{B}$ is a real block of size $4 \times 4$. We find that the diagonal blocks $\mathbf{A}$ and $\mathbf{D}$ are of the form

$$ Z^{(1)}(r, s) = \begin{pmatrix} 0 & -r & 0 & -s \\ r & 0 & -s & 0 \\ 0 & s^* & 0 & -r \\ s^* & 0 & r & 0 \end{pmatrix}, \qquad (E7) $$

where $r$ and $s$ are parameters. Specifically, the physical and virtual blocks are

$$ \mathbf{A} = Z^{(1)}(r_p, s_p), \qquad (E8a) $$

$$ \mathbf{D} = Z^{(1)}(r_v, s_v) \qquad (E8b) $$

with

$$ r_p = \frac{1 - a^4 - b^4}{1 + 2a^2 + a^4 + b^4}, \qquad (E8c) $$

$$ s_p = \frac{2b^2}{1 + 2a^2 + a^4 + b^4}, \qquad (E8d) $$

$$ r_v = \frac{1 - a^4 + b^4}{1 + 2a^2 + a^4 + b^4}, \qquad (E8e) $$

$$ s_v = \frac{2a^2b^2}{1 + 2a^2 + a^4 + b^4}. \qquad (E8f) $$

Here, the denominator is a consequence of the matrix inverse in Eq. (E5). In addition, the block containing the coupling between physical and virtual fermions is

$$ \mathbf{B} = a \begin{pmatrix} 0 & r_p + 1 & 0 & -s_p \\ r_p + 1 & 0 & s_p & 0 \\ 0 & -s_p & 0 & -r_p - 1 \\ s_p & 0 & -r_p - 1 & 0 \end{pmatrix}. \qquad (E9) $$

b. Bloch CM

We can now directly compute the Fourier Majorana CM $\tilde{\Gamma}_{|\psi\rangle}(k_x)$ of the physical state defined by the MPS on a chain with $N_x$ sites and periodic boundary conditions. $\tilde{\Gamma}_{|\psi\rangle}(k_x)$ is given by the Schur complement in Eq. (C8), where the Fourier Majorana CM $\tilde{\Gamma}_{\omega}(k_x)$ for the virtual bonds is given by Eq. (C7), and $\mathbf{A}, \mathbf{B}$, and $\mathbf{D}$ are given in the previous section.

For the SSH pumping MPS, Eq. (C8) is a matrix equation of size $4 \times 4$ since there are two physical Bloch bands. The matrix inverse can be evaluated analytically using the special parametrization $Z^{(1)}(r, s)$ from Eq. (E7). Indeed, the Fourier Majorana CM for the virtual bonds from Eq. (C7) can be written as $\tilde{\Gamma}_{\omega}(k_x) = Z^{(1)}(0, e^{ik_x})$. One easily checks that

$$ Z^{(1)}(r, s) + Z^{(1)}(r', s') = Z^{(1)}(r + r', s + s'), \qquad (E10a) $$

$$ \det[Z^{(1)}(r, s)] = (r^2 + ss^*)^2, \qquad (E10b) $$

$$ (Z^{(1)}(r, s))^{-1} = - \frac{Z^{(1)}(r, s)}{\sqrt{\det[Z^{(1)}(r, s)]}}. \qquad (E10c) $$

Using these identities, the matrix inverse in Eq. (C8) can be performed by hand,

$$ \begin{aligned} (\mathbf{D} + \tilde{\Gamma}_{\omega}(k_x))^{-1} &= (Z^{(1)}(r_v, s_v + e^{ik_x}))^{-1} \\ &= - \frac{Z^{(1)}(r_v, s_v + e^{ik_x})}{\sqrt{r_v^2 + (s_v + e^{ik_x})^2}}, \end{aligned} \qquad (E11) $$

and we evaluate the determinant $q(k_x)$ from Eq. (C9) as

$$ \begin{aligned} q(k_x) &= [r_v^2 + (s_v + e^{ik_x})^2]^2 \\ &= \frac{4(1 + a^4 + b^4 + 2a^2b^2 \cos k_x)^2}{((a^2 + 1)^2 + b^4)^2}. \end{aligned} \qquad (E12) $$

Unless $\gamma = 0$ and $|\alpha| = |\beta|$, $q(k_x)$ is strictly positive such that the MPS is well defined everywhere except for these parameter values.

After the matrix inversion in Eq. (C8), the remaining matrix multiplications in Eq. (C8) can be performed using a computer algebra system. We thus find that the Fourier Majorana CM $\tilde{\Gamma}_{|\psi\rangle}(k_x)$ of the physical state is again of the form of Eq. (E7):

$$ \tilde{\Gamma}_{|\psi\rangle}(k_x) = Z^{(1)}(r(k_x), s(k_x)) \qquad (E13a) $$

with parameters

$$ r(k_x) = \frac{1 - a^4 - b^4 - 2a^2b^2 \cos k_x}{1 + a^4 + b^4 + 2a^2b^2 \cos k_x}, \qquad (E13b) $$

$$ s(k_x) = \frac{2(b^2 + a^2e^{-ik_x})}{1 + a^4 + b^4 + 2a^2b^2 \cos k_x}. \qquad (E13c) $$

c. Parent Hamiltonian

Using Eq. (C10), we can now find a Bloch parent Hamiltonian $H_e$ for the MPS from its Majorana Fourier CM which we computed in Eq. (E13). In order to gain a physical understanding of the parent Hamiltonian, we express it in terms of the original complex physical modes before the particle-hole transformation of Eq. (7) given by

$$ \hat{a}_{A,k} = a_{A,k}, \qquad (E14a) $$

$$ \hat{a}_{B,k} = a_{B,-k}^{\dagger} \qquad (E14b) $$
after the FT of Eq. (C1). We find

$$H_{\epsilon} = \sum_{k_x} \epsilon(k_x) \begin{pmatrix} \hat{a}_{A,k_x} \\ \hat{a}_{B,k_x} \end{pmatrix}^{\dagger} \begin{pmatrix} r(k_x) & s(k_x) \\ s(k_x)^* & -r(k_x) \end{pmatrix} \begin{pmatrix} \hat{a}_{A,k_x} \\ \hat{a}_{B,k_x} \end{pmatrix}, \quad (E15)$$

where the functions $s(k_x)$ and $r(k_x)$ are defined in Eq. (E13).

In order for $H_\epsilon$ to be gapped and strictly local, we need to find a dispersion function $\epsilon(k_x)$ which is strictly positive such that $\epsilon(k_x)s(k_x)$ and $\epsilon(k_x)r(k_x)$ are polynomials in $e^{\pm ik_x}$. As explained in Appendix C2, a natural choice is $\epsilon(k_x) = q(k_x)$. Indeed, $q(k_x)$ computed in Eq. (E12) cancels the denominator of Eq. (E13). Since $q(k_x)s(k_x)$ and $q(k_x)r(k_x)$ contain the factor $e^{\pm ik_x}$ up to second order, the parent Hamiltonian $H_q$ contains hopping terms between up to second-nearest-neighbor unit cells.

Due to the special structure of Eq. (E8), we can in fact obtain a more short-ranged parent Hamiltonian from the choice

$$\epsilon(k_x) = \frac{1 + a^4 + b^4 + 2a^2b^2 \cos k_x}{a^4 + b^4 + 1} \quad (E16)$$

proportional to $\sqrt{q(k_x)}$, which is in turn proportional to the denominator of $r$ and $s$ from Eq. (E13b). $\epsilon(k_x)$ is strictly positive for all parameter values that lead to a well-defined state. The factor $1/(a^4 + b^4 + 1)$ is a normalization ensuring that the parent Hamiltonian matches Eq. (3) if the MPS is given by the parametrization $\phi_{\text{pump}}$ of Eq. (15). We see that the parent Hamiltonian $H_\epsilon$ with Bloch representation

$$H_{\epsilon}(k_x) = \frac{1}{a^4 + b^4 + 1} \left[ (1 - a^4 - b^4 - 2a^2b^2 \cos k_x) \sigma_3 + 2(b^2 + a^2 \cos k_x) \sigma_1 + 2a^2 \sin k_x \sigma_2 \right] \quad (E17)$$

with respect to the original complex physical modes has hopping only up to nearest-neighbor unit cells.

## APPENDIX F: COLUMN COVARIANCE MATRIX OF REAL-SPACE (d + 1)-DIMENSIONAL TNS FROM CHARGE PUMPING OF d-DIMENSIONAL TNS WITH CONSERVED PARTICLE NUMBER

In Sec. III C, we introduced the tensors $A^{\text{col}}$ and $B^{\text{col}}$ describing the real-space Chern PEPS restricted to a column of sites on the A and B sublattices at positions $\{(x, y)\}_{0\le y \le N_y-1}$. $A^{\text{col}}$ and $B^{\text{col}}$ are defined by the application of the inverse FT of Eq. (22) to the physical and horizontal virtual legs of a column of the hybrid Chern PEPS at positions $\{(x, k_y^{(j)})\}_{0\le j \le N_y-1}$. In this Appendix, using the representation of the SSH pumping MPS as a GfTNS from Appendix E, we compute $A^{\text{col}}$ and $B^{\text{col}}$ explicitly in terms of their CMs. In order to demonstrate the generality of the result, we will consider a general (d + 1)-dimensional TNS constructed from charge pumping of a d-dimensional TNS, which is assumed to have a conserved number of physical particles, and hence possess a U(1) symmetry of the form discussed in Appendix D.

Thus, let $|\psi_d(t)\rangle$ be the Gaussian fermionic TNS in d spatial dimensions along a cyclic interpolation parametrized

by the time $t \in (-\pi, \pi]$. With the same basis as in Eq. (B7), we collect the physical and virtual mode operators for one unit cell $\mathbf{x} \in \mathbb{Z}^d$ of $|\psi_d(t)\rangle$ into the mode vector

$$\chi_{\mathbf{x}}(t) = (a_{1,\mathbf{x}}(t), a_{2,\mathbf{x}}(t), \dots, a_{f,\mathbf{x}}(t), \\ b_{L,1,\mathbf{x}}(t), b_{L,2,\mathbf{x}}(t), \dots, b_{R,1,\mathbf{x}}(t), \dots, \\ a_{1,\mathbf{x}}^{\dagger}(t), \dots, b_{L,1,\mathbf{x}}^{\dagger}(t), \dots, b_{R,1,\mathbf{x}}^{\dagger}(t), \dots)^T \quad (\text{F1})$$

of length $2n_{\text{modes}}$. The physical and virtual mode operators now depend on the time $t$ along the interpolation. As explained in Appendix B, the d-dimensional TNS is defined by its Gaussian local fiducial state $Q_x(t)|\Omega\rangle$, which is characterized by its complex CM $G_x(t)$ of dimension $2n_{\text{modes}}$.

We assume that the TNS has a conserved number of particles, such that the local fiducial state $Q_x(t)|\Omega\rangle$ of $|\psi_d(t)\rangle$ has a U(1) symmetry of the form discussed in Appendix D. This symmetry determines which physical and virtual modes correspond to holes and particles: for each $1 \le \mu \le n_{\text{modes}}$, an entry $(\eta_{\mathbf{x}})_{\mu} = 1$ or $(\eta_{\mathbf{x}})_{\mu} = -1$ in the vector $\eta_{\mathbf{x}}$ from Eq. (D6) indicates that the mode $\mu$ has a particlelike or holelike character, respectively. Note that $\eta_{\mathbf{x}}$ does not depend on the time $t$, such that the holelike or particlelike character of the modes remains unchanged along the interpolation.

We can now move to the hybrid (d + 1)-dimensional TNS, which is defined by Eq. (20) of the main text. From Sec. III B we recall that the local fiducial state (equivalent to the local tensor) of the hybrid state at the position $(\mathbf{x}, k_{d+1}^{(j)})$ is given by $Q_x(t^{(j)})|\Omega\rangle$ containing the modes $\chi_{\mathbf{x}}(t^{(j)})$ from Eq. (F1). In particular, due to the tensor product in the (d + 1)st direction, virtual fermions in this direction are not needed.

We can now easily write the CM of one column of the hybrid (d + 1)-dimensional TNS, given by the sites at positions $\{(\mathbf{x}, k_{d+1}^{(j)})\}_{0\le j \le N_{d+1}-1}$. Indeed, due to the absence of virtual legs in the direction d + 1, the contraction of the bonds in this direction of the hybrid column amounts to a trivial tensor product in the language of local tensors. In terms of fiducial states, this corresponds to a direct sum of CMs. Hence, the complex CM $G_{\mathbf{x}}^{\text{hybrid}}$ of the column $\{(\mathbf{x}, k_{d+1}^{(j)})\}_{0\le j \le N_{d+1}-1}$ of the hybrid state is block diagonal:

$$(G_{\mathbf{x}}^{\text{hybrid}})_{\mu, \mu'; t^{(j)}, t^{(j')}} = \delta_{j, j'} (G_{\mathbf{x}}(t^{(j)}))_{\mu, \mu'}. \quad (\text{F2})$$

Here, $G_{\mathbf{x}}^{\text{hybrid}}$ is written in the basis $\{(\chi_{\mathbf{x}})_{\mu}(t^{(j)})\}_{1\le \mu \le 2n_{\text{modes}}, 0\le j \le N_{d+1}-1}$ of all operators for the physical and virtual modes in the first d directions in the column.

We now consider the (d + 1)-dimensional real-space state restricted to a column $\{(\mathbf{x}, x_{d+1})\}_{0\le x_{d+1} \le N_{d+1}-1}$, which is obtained by applying the inverse FT $\tilde{\mathcal{F}}$ in direction d + 1 to the physical and virtual legs of the hybrid column. The complex CM $G_{\mathbf{x}}^{\text{col}}$ of this at position $\mathbf{x}$ is therefore given by

$$G_{\mathbf{x}}^{\text{col}} = \tilde{\mathcal{F}}^{*} G_{\mathbf{x}}^{\text{hybrid}} \tilde{\mathcal{F}}^{T}. \quad (\text{F3})$$

The inverse FT is the d-dimensional generalization of Eq. (22),

$$\tilde{\mathcal{F}}_{\mu, \mu'; x_{d+1}, t^{(j)}} = \frac{\delta_{\mu, \mu'}}{\sqrt{N_{d+1}}} e^{i(\eta_{\mathbf{x}})_{\mu} x_{d+1} t^{(j)}} \quad (\text{F4})$$

with $0 \le j, x_{d+1} \le N_{d+1}-1$.
From Eqs. (F2), (F3), and (F4), $G_x^{\text{col}}$ becomes

$$ (G_x^{\text{col}})_{\mu, \mu', x_{d+1}, x'_{d+1}} = \frac{1}{N_{d+1}} \sum_{j=0}^{N_{d+1}-1} e^{-i[(\eta_x)_\mu x_{d+1} - (\eta_x)_{\mu'} x'_{d+1}]t^{(j)}} (G_x(t^{(j)}))_{\mu, \mu'}. \quad (\text{F5}) $$

This expression can be further simplified due to the constraint of Eq. (D8) imposed on $G_x(t^{(j)})$ by its U(1) symmetry, implying $G_x(t^{(j)})_{\mu, \mu'} = 0$ unless $(\eta_x)_\mu = (\eta_x)_{\mu'}$. Indeed, we may thus define the matrix

$$ \tilde{G}_x(t^{(j)})_{\mu, \mu'} \equiv G_x((\eta_x)_\mu t^{(j)})_{\mu, \mu'} \quad (\text{F6}) $$

that mixes elements of the complex CM of the d-dimensional state at times $t^{(j)}$ and $-t^{(j)}$ according to whether the modes $\mu$ and $\mu'$ transform as particles or holes, respectively. The $2n_{\text{modes}}$-dimensional blocks of the column CM $G_x^{\text{col}}$ are then given by the FT of this matrix:

$$ (G_x^{\text{col}})_{x_{d+1}, x'_{d+1}} = \sum_{j=0}^{N_{d+1}-1} \frac{e^{-i[x_{d+1} - x'_{d+1}]t^{(j)}}}{N_{d+1}} \tilde{G}_x(t^{(j)}). \quad (\text{F7}) $$

This expression is explicitly invariant under real-space translations $x_{d+1} \mapsto x_{d+1} + 1$ acting on both the physical modes and the virtual modes in the first d directions. Hence, the inverse FT of Eq. (F4) guarantees the translation invariance of the column CM $G_x^{\text{col}}$ in the direction $d+1$.

Finally, we want to investigate which form the global U(1) symmetry related to particle-number conservation takes for $G_x^{\text{col}}$. This result will be used in Appendix G. The U(1) is inherited from the invariance of each fiducial state $Q_x(t)|\Omega\rangle$ under the U(1) symmetry of Eq. (D3), which is independent of $t$. From the expression of the latter in second quantization [cf. Eq. (D4)], we see that the generator of the global U(1) symmetry of the column state is

$$ \left\langle \sum_{x_{d+1}=0}^{N_{d+1}-1} \left[ \sum_{\tau=0}^{f} \eta_{\tau, \mathbf{x}} a_{\tau, (\mathbf{x}, x_{d+1})}^{\dagger} a_{\tau, (\mathbf{x}, x_{d+1})} + \sum_{\alpha} \sum_{j=0}^{\xi} \eta_{\alpha, j, \mathbf{x}} b_{\alpha, j, (\mathbf{x}, x_{d+1})}^{\dagger} b_{\alpha, j, (\mathbf{x}, x_{d+1})} \right] \right\rangle = 0, \quad (\text{F8}) $$

whose expectation value vanishes.

## APPENDIX G: DISENTANGLED MODES IN SINGLE-PARTICLE ES AT FIXED PARTICLE NUMBER

In this Appendix, we show that as stated in Eq. (25), the maximal number of entangled modes compatible with the U(1) symmetry of the SSH model MPS in the ES of the Chern PEPS column state is given by $\max\{3L, N_y\}$, where L is the number of sites in the subsystem. We proceed by deriving a lower bound on the number of disentangled modes in the ES of a state with a conserved particle number (and therefore an upper bound on the number of entangled modes). We first consider a generic state in Appendix G1, before specializing to the fiducial state of a real-space column of the (d+1)-dimensional TNS in Appendix G2.

### 1. Insulator with filling fraction q

We consider a noninteracting system of N fermionic DOFs with creation and annihilation operators $a_j^\dagger, a_j$ for $j = 1, \dots, N$. We define a bipartition of the system into the subsystem $\mathcal{A}$ and its complement $\bar{\mathcal{A}}$, where the DOFs of $\mathcal{A}$ are described by the first $N_A$ modes $j = 1, \dots, N_A$. Let $|\psi\rangle$ be a pure state of this system with a conserved particle number and filling fraction $q$, such that the total number of occupied modes in $|\psi\rangle$ is $qN$. We denote by $n_A^{\lambda=1}$ the number of entanglement levels with the value $\lambda = 1$ in the single-particle ES of $|\psi\rangle$ restricted to $\mathcal{A}$. We now want to show that this number is bounded below by the filling fraction as

$$ n_A^{\lambda=1} \ge \max\{qN - (N - N_A), 0\}. \quad (\text{G1}) $$

*Proof.* Let $H = \sum_{i,j=1}^N h_{ij} a_i^\dagger a_j$ be a noninteracting flat-band Hamiltonian whose ground state is $|\psi\rangle$, where $h$ is a Hermitian matrix of dimension $N \times N$. The occupied modes in $|\psi\rangle$ are given by the $qN$ orthogonal eigenstates $u^{(k)}$ of $h$ with energy $-1$, i.e.,

$$ \sum_{j=1}^{N} h_{ij} u_{j}^{(k)} = -u_{i}^{(k)}, \quad (\text{G2}) $$

where $k = 1, \dots, qN$.

With respect to the bipartition into $\mathcal{A}$, $\bar{\mathcal{A}}$, each basis state $u^{(k)}$ falls into exactly one of the following three categories (assuming that the eigenstates are ordered accordingly):

(1) For $k = 1, \dots, m_1$, the states satisfy $u_i^{(k)} = 0$ for $i = N_A + 1, \dots, N$ such that the corresponding occupied modes are composed of DOFs of the subsystem $\mathcal{A}$. According to Theorem 1 of Ref. [69], each such state leads to an entanglement level $\lambda = 1$ in the single-particle ES of $|\psi\rangle$ restricted to $\mathcal{A}$. Hence, $m_1 \le n_A^{\lambda=1}$.

(2) For $k = m_1 + 1, \dots, m_2$, $u_i^{(k)} = 0$ for $i = 1, \dots, N_A$ such that the corresponding occupied modes are localized in $\bar{\mathcal{A}}$.

(3) For $k = m_2 + 1, \dots, qN$, $u_i^{(k)} \ne 0$ both for some $i \in \{1, \dots, N_A\}$ and some $i \in \{N_A + 1, \dots, N\}$ such that the corresponding occupied modes are localized neither in $\mathcal{A}$ nor in $\bar{\mathcal{A}}$.

The numbers $m_1, m_2$ are assumed to be maximal in the sense that no linear combination of eigenstates of the category (3) lies either purely in $\mathcal{A}$ or purely in $\bar{\mathcal{A}}$.

In the next paragraph, we will show that the number $qN - m_1$ of states from the categories (2) and (3) is no larger than the number $N - N_A$ of DOFs in $\bar{\mathcal{A}}$. This proves the claim since the number of states from the different categories can then be estimated as

$$ qN = m_1 + (qN - m_1) \le n_A^{\lambda=1} + (N - N_A), \quad (\text{G3}) $$

leading to Eq. (G1).

As a final step, we need to show that $qN - m_1 \le N - N_A$. This follows from the linear independence of the $qN - m_1$ vectors $\{\tilde{u}^{(k)}\}_{k \ge m_1}$ of dimension $N - N_A$, where $\tilde{u}_i^{(k)} = u_{N_A+i}^{(k)}$ for $i = N_A + 1, \dots, N$ is the restriction of the eigenvector to the DOFs of $\bar{\mathcal{A}}$. Indeed, let us assume that we have a vanishing
linear combination

$$0 = \sum_{k=m_1+1}^{qN} \mu_k \tilde{u}^{(k)} \qquad (G4)$$

with coefficients $\mu_k$. The scalar product with $\tilde{u}^{(l)}$ shows that $0 = \mu_l$ for $l = m_1 + 1, \dots, m_2$: indeed, the orthogonality of the $u^{(k)}$ implies $\sum_{i>N_A} \tilde{u}_i^{(l)} \tilde{u}_i^{(k)} = \sum_{i \ge 1} u_i^{(l)} u_i^{(k)} = \delta_{kl}$. Equation (G4) therefore implies $0 = \sum_{k>m_2}^{qN} \mu_k \tilde{u}^{(k)}$ where the sum runs only over category (3). Due to the maximality of $m_1, m_2$ as defined above just below point (3), we then have $0 = \sum_{k>m_2}^{qN} \mu_k u^{(k)}$. Otherwise, this would be a linear combination purely in $\mathcal{A}$ since the part in $\bar{\mathcal{A}}$ vanishes by assumption. We thus also get $0 = \mu_l$ for $l = m_2 + 1, \dots, qN$, proving that $qN - m_1 \le N - N_A$.

## 2. Real-space column of (d + 1)-dimensional TNS

We now apply Eq. (G1) to the case where $|\psi\rangle$ is the fiducial state of one real-space column of the (d + 1)-dimensional pumping TNS, whose CM is computed in Eq. (F5). The column state has $N = n_{\text{modes}}N_{d+1}$ degrees of freedom, where $N_{d+1}$ is the number of sites in the direction $d+1$ and $n_{\text{modes}}$ the number of physical and virtual particles per lattice site of the d-dimensional TNS. Below we will show that in a suitable basis, the column state has a conserved particle number $qN = N_{d+1}n_{\eta_{-}}$, where $n_{\eta_{-}}$ is the number of values $-1$ in the first $n_{\text{modes}}$ entries of the vector $\eta_{\chi}$ from Eq. (D6). We consider the single-particle ES of the column state with respect to the subsystem $\mathcal{A}_L$ of the first $L$ sites $0 \le x_{d+1} \le L-1$ of the column, which has $N_{\mathcal{A}_L} = Ln_{\text{modes}}$ DOFs. By Eq. (G1), the number $n_{\mathcal{A}_L}^{\lambda=1}$ of entanglement levels with the value $\lambda = 1$ in this spectrum is lower bounded as

$$n_{\mathcal{A}_L}^{\lambda=1} \ge \max\{N_{d+1}n_{\eta_{-}} - (N_{d+1} - L)n_{\text{modes}}, 0\}. \qquad (G5)$$

We now show that with a suitable particle-hole transformation we can find a single-particle basis with respect to which the column state has a conserved particle number $N_{d+1}n_{\eta_{-}}$. In Eq. (F5) the state is expressed in a basis where the CM $G_{\chi}^{\text{col}}$ has a nonvanishing off-diagonal block $Q$ corresponding to superconducting correlations, such that only the parity of the particle number is conserved [cf. Eq. (A1)]. Below Eq. (F5) we note that half the entries of $G_{\chi}^{\text{col}}$ vanish, namely, those with $(\eta_{\chi})_{\mu} \ne (\eta_{\chi})_{\mu'}$. A new basis, in which the complex CM has only a diagonal block $R$ and hence a fixed particle number, is created as follows: For those modes $1 \le m \le n_{\text{modes}}$ with $(\eta_{\chi})_m = -1$, we exchange annihilation and creation operators since the latter have $(\eta_{\chi})_{m+n_{\text{modes}}} = 1$. This is a particle-hole transformation corresponding to the mapping

$$a_{\tau, (x, x_{d+1})} \mapsto a_{\tau, (x, x_{d+1})}^{\dagger} \qquad (G6)$$

for all physical modes with $\eta_{\tau,x} = -1$ while leaving modes with $\eta_{\tau,x} = 1$ unchanged, and similarly for the virtual modes.

The number of particles in the column state is given by the expectation value from Eq. (F8) of the generator of its global U(1) symmetry. Note that this is not the physical particle number, but rather the number of particles in the system composed of the physical legs and the virtual legs in the first $d$

spatial dimensions. Under the transformation of Eq. (G6), the particle-number expectation value from Eq. (F8) transforms as explained in Eq. (13). We thus find that in the new basis of modes the particle number is $N_{d+1}n_{\eta_{-}}$ as claimed above.

### a. Chern PEPS

Let us apply these results to the Chern PEPS derived from the SSH pumping MPS. For simplicity, we restrict ourselves to a column $A^{\text{col}}$ of A sites with $n_{\text{modes}} = 3$ DOFs per site and $\eta_A = 1$ and $\eta_{L,A} = \eta_{R,A} = -1$ [cf. the discussion below Eq. (13)]. Therefore, $n_{\eta_{-}} = 2$ and Eq. (G5) becomes

$$n_{\mathcal{A}_L}^{\lambda=1} \ge \max\{3L - N_y, 0\}. \qquad (G7)$$

This shows that Eq. (25) gives the maximal number of entangled modes compatible with the U(1) symmetry of the SSH model MPS.

For the PEPS defined by the trivial cycle $\phi_{\text{triv}}$ from Eq. (17), the discussion above can be refined: since $\beta = 0$ throughout the interpolation, all right virtual modes decouple from the column tensor such that there trivially are $N_y$ entanglement levels $|\lambda| = 1$. Let us investigate if the U(1) symmetry causes additional decoupled levels in the system of the coupled physical and left virtual particles. This system has only two DOFs per site, namely, the physical leg with $\eta_A = 1$ and the left virtual leg with $\eta_{L,A} = -1$. In this case, Eq. (G5) gives a trivial lower bound for the number of decoupled modes with levels $\lambda = 1$:

$$n_{\mathcal{A}_L}^{\lambda=1} \ge \max\{2L - N_y, 0\} = 0 \quad \text{for } L \le N_y/2. \qquad (G8)$$

Hence, the U(1) symmetry does not cause any additional decoupled modes, and the number of entangled modes is given by $2L$ as discussed in the main text.

### b. Chiral hinge PEPS

The discussion for the three-dimensional chiral hinge PEPS from Sec. IV C is analogous, where a column of sites on the sublattice 1 has $n_{\text{modes}} = 5$ DOFs per site with $\eta_1 = 1$ and $\eta_{L,1} = \eta_{U,1} = \eta_{R,1} = \eta_{D,1} = -1$. Therefore, $n_{\eta_{-}} = 4$ and Eq. (G5) implies $n_{\mathcal{A}_L}^{\lambda=1} \ge \max\{5L - N_z, 0\}$.

For the mirror-symmetric case with couplings $\alpha = \alpha_x = \alpha_y$ and $\beta = \beta_x = \beta_y$, this bound can be refined. Indeed, by defining the linear combinations $b_{\pm,1}^{LD} = (b_{L,1} \pm b_{D,1})/\sqrt{2}$ and $b_{\pm,1}^{UR} = (b_{R,1} \pm b_{R,1})/\sqrt{2}$, the local fiducial state $|Q^{[1]}\rangle$ from Eq. (H1) below can be written as

$$|Q^{[1]}\rangle = [\gamma - \sqrt{2}\beta a_1^\dagger(b_{-,1}^{LD})^\dagger - \sqrt{2}\alpha a_1^\dagger(b_{-,1}^{UR})^\dagger]|\Omega\rangle. \qquad (G9)$$

Therefore, when considering only one column of sites on the sublattice 1, two virtual fermionic modes decouple from the local tensor on each site. Effectively, the remaining coupled system therefore has only $n_{\text{modes}} = 3$ DOFs per site with $\eta_1 = 1$ and $\eta_{-1}^{LD} = \eta_{-1}^{UR} = -1$, such that the number of additional disentangled modes due to the U(1) symmetry can be estimated as above for the Chern PEPS. Hence, with the identification $N_y \mapsto N_z$, the bound on the number of disentangled modes is given by Eqs. (G7) and (G8) for the PEPSs derived from $\phi_{\text{pump}}$ and $\phi_{\text{triv}}$, respectively.
## APPENDIX H: QUADRUPOLINE PEPS AS GfTNS

In this Appendix we apply the formalism of GfTNSs to the quadrupole model pumping PEPS of Eq. (29). In Appendix H1, we show how the state can be expressed as a GfTNS. This allows us to compute its Bloch CM and a Bloch parent Hamiltonian on the torus in Appendix H2. These sections are completely analogous to Appendix E for the SSH pumping MPS. We therefore refer the reader to this Appendix for a detailed explanation of each step. Finally, in Appendix H3 we discuss the ES of the PEPS when the parameters are chosen such that the state represents the OAI dimerized phase of the quadrupole model.

### 1. Expression as GfTNS

The quadrupole model pumping PEPS is defined in Eq. (29) of the main text in the language of local tensors. Here, we want to reexpress this state in the formalism of GfTNSs. We recall that each unit cell consists of 2 × 2 lattice sites, and that the PEPS has physical dimension 2 and bond dimension 2. This corresponds to $f = 1$ physical fermion per lattice site and $\xi = 1$ virtual fermion per nearest-neighbor bond and lattice site, represented in Fig. 6(b) by blue and red circles, respectively.

As explained in Appendix B, the local PEPS tensors $A^{[\tau]}$ on the four sublattices $\tau = 1, 2, 3, 4$ from the main text correspond to local fiducial states $|Q^{[\tau]}\rangle$, whose basis coefficients are given by the local tensors [see Eq. (B10)]. We write $a_\tau$ for the annihilation operator of the physical fermion and $b_{L,\tau}$, $b_{U,\tau}$, $b_{R,\tau}$, $b_{D,\tau}$ for the annihilation operators of the left, up, right, down virtual fermions on the sublattice $\tau$ (we dropped the unit-cell index due to the translation invariance). Let $|\Omega\rangle$ denote the vacuum annihilated by all these operators. Applying Eq. (B10), we see that the local fiducial states on the four sublattices derived from the local tensors of Eq. (29) are

$$|Q^{[1]}\rangle = [\gamma - \beta_x a_1^\dagger b_{L,1}^\dagger - \alpha_y a_1^\dagger b_{U,1}^\dagger + \alpha_x a_1^\dagger b_{R,1}^\dagger \\ + \beta_y a_1^\dagger b_{D,1}^\dagger] |\Omega\rangle, \qquad (H1a)$$

$$|Q^{[2]}\rangle = [\gamma - \alpha_x a_2^\dagger b_{L,2}^\dagger + \beta_y a_2^\dagger b_{U,2}^\dagger + \beta_x a_2^\dagger b_{R,2}^\dagger \\ - \alpha_y a_2^\dagger b_{D,2}^\dagger] |\Omega\rangle, \qquad (H1b)$$

$$|Q^{[3]}\rangle = [\gamma + \alpha_x a_3^\dagger b_{L,3}^\dagger + \alpha_y a_3^\dagger b_{U,3}^\dagger + \beta_x a_3^\dagger b_{R,3}^\dagger \\ + \beta_y a_3^\dagger b_{D,3}^\dagger] |\Omega\rangle, \qquad (H1c)$$

$$|Q^{[4]}\rangle = [\gamma + \beta_x a_4^\dagger b_{L,4}^\dagger + \beta_y a_4^\dagger b_{U,4}^\dagger + \alpha_x a_4^\dagger b_{R,4}^\dagger \\ + \alpha_y a_4^\dagger b_{D,4}^\dagger] |\Omega\rangle. \qquad (H1d)$$

These fiducial states are of the form of Eq. (A11) with only zero- and second-order terms in the creation operators. Hence, the fiducial states are Gaussian and can be

parametrized as in Eqs. (A10) and (B12) with antisymmetric coefficient matrices

$$M_1 = \frac{1}{2 \times 2^{1/4}} \begin{pmatrix} 0 & -b_x & -a_y & a_x & b_y \\ b_x & 0 & 0 & 0 & 0 \\ a_y & 0 & 0 & 0 & 0 \\ -a_x & 0 & 0 & 0 & 0 \\ -b_y & 0 & 0 & 0 & 0 \end{pmatrix}, \qquad (H2a)$$

$$M_2 = \frac{1}{2 \times 2^{1/4}} \begin{pmatrix} 0 & -a_x & b_y & b_x & -a_y \\ a_x & 0 & 0 & 0 & 0 \\ -b_y & 0 & 0 & 0 & 0 \\ -b_x & 0 & 0 & 0 & 0 \\ a_y & 0 & 0 & 0 & 0 \end{pmatrix}, \qquad (H2b)$$

$$M_3 = \frac{1}{2 \times 2^{1/4}} \begin{pmatrix} 0 & a_x & a_y & b_x & b_y \\ -a_x & 0 & 0 & 0 & 0 \\ -a_y & 0 & 0 & 0 & 0 \\ -b_x & 0 & 0 & 0 & 0 \\ -b_y & 0 & 0 & 0 & 0 \end{pmatrix}, \qquad (H2c)$$

$$M_4 = \frac{1}{2 \times 2^{1/4}} \begin{pmatrix} 0 & b_x & b_y & a_x & a_y \\ -b_x & 0 & 0 & 0 & 0 \\ -b_y & 0 & 0 & 0 & 0 \\ -a_x & 0 & 0 & 0 & 0 \\ -a_y & 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (H2d)$$

Here, we defined the quotients $a_x = \alpha_x/\gamma$, $a_y = \alpha_y/\gamma$, $b_x = \beta_x/\gamma$, and $b_y = \beta_y/\gamma$ of the parameters from Eq. (29). We absorbed the remaining factor $\gamma$ into the normalization constant $N$ in Eq. (A11).

We have thus successfully written the PEPS from Eq. (29) as a GfTNS. Before proceeding, we choose the orientation of the virtual bonds as follows: all horizontal bonds are oriented from left to right and all vertical bonds are oriented from top to bottom.

### 2. Bloch CM and parent Hamiltonian

Having expressed the quadrupole model pumping PEPS as a GfTNS, we now want to use the formalism from Appendix C to compute the Bloch CM and a Bloch parent Hamiltonian for the state on a torus. We will compute the Bloch CM by evaluating Eq. (C8), and from there obtain the parent Hamiltonian via Eq. (C10).

#### a. CM of unit cell

The Majorana CM $\Gamma_{Q_x}$ in Eq. (C8) refers to the fiducial state of a unit cell, not that of a single site. We therefore need to compute $\Gamma_{Q_x}$ from the fiducial states for each individual sublattice given in Eq. (H1). This is done by contracting the four virtual bonds within one unit cell. In the language of fiducial states, we project the tensor product $|Q^{[1]}\rangle \otimes |Q^{[2]}\rangle \otimes |Q^{[3]}\rangle \otimes |Q^{[4]}\rangle$ of the fiducial states for a unit cell on the product of the maximally entangled states of the four virtual bonds connecting the lattice sites within this unit cell. This is analogous to the computation leading to Eq. (E13), and we refer the reader to Appendix E 2 a for details. Here, we will only state the result for the Majorana CM $\Gamma_{Q_x}$.
<figure><img src="image_10.png" /><figcaption>FIG. 10. Unit cell of the quadrupole model pumping PEPS after projection of the virtual fermions corresponding to the bonds within the unit cell (marked by green ellipses) onto maximally entangled states. This results in a fiducial state with four physical fermions (blue circles) collected into the block p, four horizontal virtual fermions in the block $v_h$, and four virtual fermions in the block $v_v$ (red circles for both). The order of the individual fermions within these blocks is indicated by the labels next to each circle.</figcaption></figure>

The fiducial state of one unit cell, after contraction of the virtual bonds within the plaquette, describes four physical fermions and two virtual fermions per direction left, up, right, down (see Fig. 10). We collect the four physical, horizontal virtual and vertical virtual fermions into blocks labeled p, $v_h$, and $v_v$. The order of the individual fermions within these blocks is indicated in Fig. 10. Since the Majorana CM $\Gamma_{Q_x}$ of the fiducial state of a unit cell is antisymmetric (see Appendix A 1), we can then separate it into 8 × 8 blocks as

$$ \Gamma_{Q_x} = \begin{pmatrix} \mathbf{A}_p & \mathbf{B}_{pv_h} & \mathbf{B}_{pv_v} \\ -\mathbf{B}_{pv_h}^T & \mathbf{D}_{v_h} & \mathbf{B}_{v_h v_v} \\ -\mathbf{B}_{pv_v}^T & -\mathbf{B}_{v_h v_v}^T & \mathbf{D}_{v_v} \end{pmatrix}. \quad (H3) $$

The blocks $A_p$, $D_{v_h}$, and $D_{v_v}$ describe the reduced fiducial state of the physical, horizontal virtual and vertical virtual subsystem, respectively. The off-diagonal blocks describe the

coupling between these three subsystems. Equation (H3) is a generalization of Eq. (C6) with

$$ \mathbf{A} = \mathbf{A}_p, \quad (H4a) $$

$$ \mathbf{B} = (\mathbf{B}_{pv_h} \ \mathbf{B}_{pv_v}), \quad (H4b) $$

$$ \mathbf{D} = \begin{pmatrix} \mathbf{D}_{v_h} & \mathbf{B}_{v_h v_v} \\ -\mathbf{B}_{v_h v_v}^T & \mathbf{D}_{v_v} \end{pmatrix} \quad (H4c) $$

which we find useful since we want to treat horizontal and vertical virtual fermions separately.

We are now ready to give the expression for the Majorana CM $\Gamma_{Q_x}$ of the fiducial state of a unit cell. For simplicity, we set $\alpha_x = \alpha_y = \alpha$ and $\beta_x = \beta_y = \beta$ from now on, corresponding to parameters $a_x = a_y = a$ and $b_x = b_y = b$ in Eq. (H2). The diagonal blocks $A_p$, $D_{v_h}$, and $D_{v_v}$ of $\Gamma_{Q_x}$ take the form

$$ Z^{(2)}(r, s, u) = \begin{pmatrix} 0 & -r & 0 & -u & 0 & -s & 0 & 0 \\ r & 0 & -u & 0 & -s & 0 & 0 & 0 \\ 0 & u^* & 0 & -r & 0 & 0 & 0 & -s \\ u^* & 0 & r & 0 & 0 & 0 & -s & 0 \\ 0 & s^* & 0 & 0 & 0 & -r & 0 & u \\ s^* & 0 & 0 & 0 & r & 0 & u & 0 \\ 0 & 0 & 0 & s^* & 0 & -u^* & 0 & -r \\ 0 & 0 & s^* & 0 & -u^* & 0 & r & 0 \end{pmatrix} \quad (H5) $$

with parameters $r \in \mathbb{R}$ and $s, u \in \mathbb{C}$. Specifically,

$$ \mathbf{A}_p = Z^{(2)}(r_p, -s_p, s_p), \quad (H6a) $$

$$ \mathbf{D}_{v_h} = Z^{(2)}(r_v, s_v, s_v), \quad (H6b) $$

$$ \mathbf{D}_{v_v} = Z^{(2)}(r_v, -s_v, s_v) \quad (H6c) $$

with parameters

$$ r_p = \frac{\sqrt{2}(2a^4 + b^4 - 1)}{4a^2 + 2\sqrt{2}a^4 + \sqrt{2}(1 + b^4)}, \quad (H7a) $$

$$ s_p = \frac{2b^2}{4a^2 + 2\sqrt{2}a^4 + \sqrt{2}(1 + b^4)}, \quad (H7b) $$

$$ r_v = \frac{2a^2 + \sqrt{2}(b^4 + 1)}{4a^2 + 2\sqrt{2}a^4 + \sqrt{2}(1 + b^4)}, \quad (H7c) $$

$$ s_v = \frac{\sqrt{2}a^2b^2}{4a^2 + 2\sqrt{2}a^4 + \sqrt{2}(1 + b^4)}, \quad (H7d) $$

where the denominator stems from the matrix inversion in the Schur complements used to evaluate the projection on the virtual bonds within one unit cell. The off-diagonal blocks of $\Gamma_{Q_x}$ are

$$ \mathbf{B}_{v_h v_v} = \begin{pmatrix} 0 & 1 - r_v & 0 & s_v & 0 & -s_v & 0 & 0 \\ r_v - 1 & 0 & s_v & 0 & -s_v & 0 & 0 & 0 \\ 0 & s_v & 0 & 0 & 0 & 1 - r_v & 0 & -s_v \\ s_v & 0 & 0 & 0 & r_v - 1 & 0 & -s_v & 0 \\ 0 & s_v & 0 & r_v - 1 & 0 & 0 & 0 & s_v \\ s_v & 0 & 1 - r_v & 0 & 0 & 0 & s_v & 0 \\ 0 & 0 & 0 & s_v & 0 & s_v & 0 & 1 - r_v \\ 0 & 0 & s_v & 0 & s_v & 0 & r_v - 1 & 0 \end{pmatrix}, \quad (H8a) $$
$$
\mathbf{B}_{pv_h} = \frac{a}{2^{1/4}}
\begin{pmatrix}
0 & r_p+1 & 0 & -s_p & 0 & -s_p & 0 & 0 \\
r_p+1 & 0 & s_p & 0 & s_p & 0 & 0 & 0 \\
0 & s_p & 0 & 0 & 0 & r_p+1 & 0 & s_p \\
-s_p & 0 & 0 & 0 & r_p+1 & 0 & -s_p & 0 \\
0 & -s_p & 0 & -r_p-1 & 0 & 0 & 0 & s_p \\
s_p & 0 & -r_p-1 & 0 & 0 & 0 & -s_p & 0 \\
0 & 0 & 0 & s_p & 0 & -s_p & 0 & r_p+1 \\
0 & 0 & -s_p & 0 & s_p & 0 & r_p+1 & 0
\end{pmatrix},
\quad \text{(H8b)}
$$

$$
\mathbf{B}_{pv_v} = \frac{a}{2^{1/4}}
\begin{pmatrix}
0 & r_p+1 & 0 & s_p & 0 & -s_p & 0 & 0 \\
r_p+1 & 0 & -s_p & 0 & s_p & 0 & 0 & 0 \\
0 & s_p & 0 & -r_p-1 & 0 & 0 & 0 & s_p \\
-s_p & 0 & -r_p-1 & 0 & 0 & 0 & -s_p & 0 \\
0 & -s_p & 0 & 0 & 0 & -r_p-1 & 0 & s_p \\
s_p & 0 & 0 & 0 & -r_p-1 & 0 & -s_p & 0 \\
0 & 0 & 0 & s_p & 0 & s_p & 0 & r_p+1 \\
0 & 0 & -s_p & 0 & -s_p & 0 & r_p+1 & 0
\end{pmatrix}.
\quad \text{(H8c)}
$$

b. Bloch CM

We can now compute the Bloch CM $\tilde{\Gamma}_{|\psi)}(\mathbf{k})$ of the physical state on a torus by evaluating Eq. (C8). Here, the Majorana CM of the local fiducial state of the unit cell with its blocks **A**, **B**, and **D** is given in the previous subsection.

In the same basis as **D**, the Bloch CM $\tilde{\Gamma}_{\omega}(\mathbf{k})$ of the total virtual maximally entangled state of size $16 \times 16$ is

$$
\tilde{\Gamma}_{\omega}(\mathbf{k}) = \begin{pmatrix} \tilde{\Gamma}_{\omega, v_h}(\mathbf{k}) & 0 \\ 0 & \tilde{\Gamma}_{\omega, v_v}(\mathbf{k}) \end{pmatrix}, \quad \text{(H9a)}
$$

where the blocks $\tilde{\Gamma}_{\omega, v_h}(\mathbf{k})$ and $\tilde{\Gamma}_{\omega, v_v}(\mathbf{k})$ refer to the horizontal and vertical virtual fermions, respectively. Since the horizontal fermions on the sublattice 1 couple only to the horizontal fermions on sublattice 3 (and similarly for 2 and 4), the CM $\tilde{\Gamma}_{\omega, v_h}(\mathbf{k})$ is given by a sum of two copies of the CM from Eq. (C7) for an MPS with $\xi = 1$. Hence,

$$
\tilde{\Gamma}_{\omega, v_h}(\mathbf{k}) = \begin{pmatrix} 0 & 0 & -\sigma_1 e^{-ik_x} & 0 \\ 0 & 0 & 0 & -\sigma_1 e^{-ik_x} \\ \sigma_1 e^{ik_x} & 0 & 0 & 0 \\ 0 & \sigma_1 e^{ik_x} & 0 & 0 \end{pmatrix}. \quad \text{(H9b)}
$$

The expression for $\tilde{\Gamma}_{\omega, v_v}(\mathbf{k})$ is analogous with $k_x \to -k_y$ (where the negative sign indicates that the bonds are oriented from top to bottom and hence point in the direction of negative $y$). Note that $\tilde{\Gamma}_{\omega, v_h}(\mathbf{k}) = Z^{(2)}(0, e^{-ik_x}, 0)$ and $\tilde{\Gamma}_{\omega, v_v}(\mathbf{k}) = Z^{(2)}(0, e^{ik_x}, 0)$ are of the form of Eq. (H5).

For the quadrupole model pumping PEPS, the evaluation of Eq. (C8) requires the inversion of the matrix $\mathbf{D} + \tilde{\Gamma}_{\omega}(\mathbf{k})$ of size $16 \times 16$. This can be done analytically using the special representation from Eq. (H5). Indeed, one can show that the matrix $Z^{(2)}$ has the properties

$$
Z^{(2)}(r, s, u) + Z^{(2)}(r', s', u') = Z^{(2)}(r + r', s + s', u + u'), \quad \text{(H10a)}
$$

$$
\det[Z^{(2)}(r, s, u)] = (r^2 + ss^* + uu^*)^4, \quad \text{(H10b)}
$$

$$
[Z^{(2)}(r, s, u)]^{-1} = - \frac{Z^{(2)}(r, s, u)}{[\det Z^{(2)}(r, s, u)]^{1/4}}. \quad \text{(H10c)}
$$

Moreover, conjugation with the matrix $\mathbf{B}_{v_h v_v}$ from Eq. (H8a), which gives the off-diagonal block of **D**, returns a matrix of the form of Eq. (H5),

$$
\mathbf{B}_{v_h v_v} Z^{(2)}(r, s, u) \mathbf{B}_{v_h v_v}^T = Z^{(2)}(r', s', u'), \quad \text{(H11)}
$$

where $r' = r[(1 - r_v)^2 - 2s_v^2] + 2(1 - r_v)s_v \text{Re}(s - u)$, $s' = -s(1 - r_v)^2 + 2rs_v(r_v - 1) - 2s_v^2[\text{Re}(u) + i\text{Im}(s)]$, and $u' = u(1 - r_v)^2 + 2rs_v(r_v - 1) + 2s_v^2[\text{Re}(s) + i\text{Im}(u)]$. Using these identities, the inverse of $\mathbf{D} + \tilde{\Gamma}_{\omega}(\mathbf{k})$ can be evaluated blockwise, and we compute the determinant in Eq. (C9) as

$$
q(\mathbf{k}) = \frac{2^8(1 + a^4 + b^4 + a^2b^2 \cos k_x + a^2b^2 \cos k_y)^4}{(1 + 2\sqrt{2}a^2 + 2a^4 + b^4)^4}. \quad \text{(H12)}
$$

This is strictly positive unless $\gamma = 0$ and $|\alpha| = |\beta|$, such that the PEPS is well defined everywhere except for these parameter values.

Proceeding thus, we find that the Bloch Majorana CM $\tilde{\Gamma}_{|\psi)}(\mathbf{k})$ of the physical state is also of the form of Eq. (H5),

$$
\tilde{\Gamma}_{|\psi\rangle} = Z^{(2)}(r(\mathbf{k}), s(\mathbf{k}), u(\mathbf{k})), \quad \text{(H13a)}
$$

where the parameters are

$$
r(\mathbf{k}) = \frac{1 - a^4 - b^4 - a^2 b^2 \cos k_x - a^2 b^2 \cos k_y}{1 + a^4 + b^4 + a^2 b^2 \cos k_x + a^2 b^2 \cos k_y}, \quad \text{(H13b)}
$$

$$
s(\mathbf{k}) = \frac{-\sqrt{2}(b^2 + a^2 e^{-ik_x})}{1 + a^4 + b^4 + a^2 b^2 \cos k_x + a^2 b^2 \cos k_y}, \quad \text{(H13c)}
$$

$$
u(\mathbf{k}) = \frac{\sqrt{2}(b^2 + a^2 e^{ik_y})}{1 + a^4 + b^4 + a^2 b^2 \cos k_x + a^2 b^2 \cos k_y}. \quad \text{(H13d)}
$$

In the denominator, we recognize the fourth root $q(\mathbf{k})^{1/4}$ coming from the matrix inversion according to Eq. (H10c).

c. Parent Hamiltonian

We are now in a position to find a parent Hamiltonian for the physical state on a torus, using Eq. (C10) and the result for the Bloch Majorana CM $\tilde{\Gamma}_{|\psi)}$ from Eq. (H13). To
get an intuitive understanding of the Hamiltonian, we find it useful to express it in terms of the original fermionic modes before the particle-hole transformation of Eq. (28). Their FT is related to the FT of the new modes as $\hat{a}_{\tau, k} = a_{\tau, k}$ for $\tau = 1, 2$ and $\hat{a}_{\tau, k} = a_{\tau, -k}^\dagger$ for $\tau = 3, 4$. With respect to the original fermionic modes, Eq. (C10) is given by

$$H_{\epsilon} = \sum_{\mathbf{k}} \epsilon(\mathbf{k}) \times \hat{a}_{\tau, \mathbf{k}}^{\dagger} \begin{pmatrix} r(\mathbf{k}) & 0 & -s(\mathbf{k})^{*} & u(\mathbf{k}) \\ 0 & r(\mathbf{k}) & -u(\mathbf{k})^{*} & -s(\mathbf{k}) \\ -s(\mathbf{k}) & -u(\mathbf{k}) & -r(\mathbf{k}) & 0 \\ u(\mathbf{k})^{*} & -s(\mathbf{k})^{*} & 0 & -r(\mathbf{k}) \end{pmatrix}_{\tau \tau'} \hat{a}_{\tau', \mathbf{k}}, \quad (H14)$$

where a summation over $\tau, \tau'$ is implied in the second line.

As discussed in Appendix H2c, the properties of the parent Hamiltonian are determined by our choice of dispersion relation $\epsilon(\mathbf{k})$. By setting $\epsilon(\mathbf{k}) = q(\mathbf{k})$, we would obtain a gapped parent Hamiltonian with coupling between up to fourth-nearest-neighbor unit cells. However, we can find a more short-ranged parent Hamiltonian due to the special form of the fiducial state CM using Eq. (H5). Indeed, by choosing the dispersion function

$$\epsilon(\mathbf{k}) = \frac{1 + a^4 + b^4 + a^2 b^2 \cos k_x + a^2 b^2 \cos k_y}{a^4 + b^4 + 1} \quad (H15)$$

which is proportional to the denominator of $r, s$, and $u$ from Eq. (H13), we obtain a second-nearest-neighbor parent Hamiltonian with Bloch representation

$$H_{\epsilon}(\mathbf{k}) = \frac{1}{a^4 + b^4 + 1} \times [(1 - a^4 - b^4 - a^2 b^2 \cos k_x - a^2 b^2 \cos k_y)\sigma_3 \otimes \sigma_0 + \sqrt{2}(b^2 + a^2 \cos k_x)\sigma_1 \otimes \sigma_0 + \sqrt{2}a^2 \sin k_x(-\sigma_2 \otimes \sigma_3) + \sqrt{2}(b^2 + a^2 \cos k_y)(-\sigma_2 \otimes \sigma_2) + \sqrt{2}a^2 \sin k_y(-\sigma_2 \otimes \sigma_1)], \quad (H16)$$

where $\sigma_0$ is the identity matrix of dimension two.

### 3. ES in topological quadrupole phase

When $b=0$, the system described by the PEPS of Eq. (H2) splits into decoupled four-site plaquettes shifted from the

unit cell by one site in both directions, corresponding to the OAI dimerized phase of the quadrupole model if $a=1$. The Majorana CM $\Gamma_{\text{Plaquette}}$ describing the physical state of one such decoupled plaquette takes the form of Eq. (H5) with

$$r = \frac{1 - a^4}{1 + a^4}, \quad (H17a)$$

$$u = s = - \frac{\sqrt{2} a^2}{1 + a^4}. \quad (H17b)$$

As an application, we will derive the ES contributions from the edges and corners given in Eqs. (34a) and (34b), respectively.

We begin with the four corners. The Majorana CM $\Gamma_{\text{corner}, \tau}$ for a single corner site on the sublattice $\tau$ is given by the corresponding block of dimension 2 on the diagonal of $\Gamma_{\text{Plaquette}}$. Specifically,

$$\Gamma_{\text{corner}, \tau} = \begin{pmatrix} 0 & -r \\ r & 0 \end{pmatrix} \quad (H18)$$

for $\tau = 1, 2, 3, 4$. We now transform $\Gamma_{\text{corner}, \tau}$ to the basis of the original complex fermionic modes before the particle-hole transformation of Eq. (28). Then, the CM of the corner site has a vanishing off-diagonal block $\hat{Q}_{\text{corner}, \tau}^* = 0$ and a diagonal block $\hat{R}_{\text{corner}, \tau}^* = \mp \frac{i}{2} \lambda_{\text{corner}}$ with $\lambda_{\text{corner}} = r$, where the negative and positive signs hold for the sublattices $\tau = 1, 2$ and $\tau = 3, 4$, respectively. Using the expression for $r$ from Eq. (H17a) and $a = \alpha/\gamma$, we obtain the formula for the corner ES level given in Eq. (34b) with one level per corner.

Similarly, the Majorana CM $\Gamma_{\text{edge}, \tau_1, \tau_2}$ for two decoupled edge sites on the sublattices $\tau_1$ and $\tau_2$ is given by a block of dimension 4 of $\Gamma_{\text{Plaquette}}$. Concretely,

$$\Gamma_{\text{edge}, \tau_1 \tau_2} = \begin{pmatrix} 0 & -r & 0 & -s \\ r & 0 & -s & 0 \\ 0 & s & 0 & -r \\ s & 0 & r & 0 \end{pmatrix} \quad (H19)$$

for $(\tau_1, \tau_2) \in \{(3, 1), (4, 2), (4, 1), (2, 3)\}$. In the basis of the original complex fermionic modes before the particle-hole transformation of Eq. (28), $\Gamma_{\text{edge}, \tau_1, \tau_2}$ takes the following form: it has a vanishing off-diagonal block $\hat{Q}_{\text{edge}, \tau_1 \tau_2}^* = 0$ and a nonzero diagonal block $\hat{R}_{\text{edge}, \tau_1 \tau_2}^*$ with doubly degenerate eigenvalues $\pm \frac{i}{2} \lambda_{\text{edge}}$ and

$$\lambda_{\text{edge}} = \sqrt{r^2 + s^2} = \frac{\sqrt{1 + a^8}}{1 + a^4}. \quad (H20)$$

This corresponds to the formula given in Eq. (34a).

[1] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science **357**, 61 (2017).

[2] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv. **4**, eaat0346 (2018).

[3] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B **96**, 245115 (2017).

[4] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. **119**, 246401 (2017).

[5] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. **82**, 3045 (2010).

[6] F. D. M. Haldane, Phys. Rev. Lett. **61**, 2015 (1988).

[7] B. Bradlyn, L. Elcoro, J. Cano, M. Vergniory, Z. Wang, C. Felser, M. Aroyo, and B. A. Bernevig, Nature (London) **547**, 298 (2017).

[8] B. J. Wieder and B. A. Bernevig, arXiv:1810.02373.

[9] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. B **22**, 2099 (1980).
[10] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nat. Phys. **14**, 918 (2018).

[11] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature (London) **555**, 342 (2018).

[12] H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Nat. Mater. **18**, 108 (2019).

[13] X. Ni, M. Weiner, A. Alù, and A. B. Khanikaev, Nat. Mater. **18**, 113 (2019).

[14] S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Nat. Photonics **13**, 692 (2019).

[15] A. E. Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, and M. Bourennane, Nat. Photonics **13**, 697 (2019).

[16] B.-Y. Xie, G.-X. Su, H.-F. Wang, H. Su, X.-P. Shen, P. Zhan, M.-H. Lu, Z.-L. Wang, and Y.-F. Chen, Phys. Rev. Lett. **122**, 233903 (2019).

[17] Y. Yang, Z. Jia, Y. Wu, R.-C. Xiao, Z. H. Hang, H. Jiang, and X. C. Xie, Sci. Bulletin (2020), doi:10.1016/j.scib.2020.01.024.

[18] C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature (London) **555**, 346 (2018).

[19] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Nat. Phys. **14**, 925 (2018).

[20] M. Serra-Garcia, R. Süsstrunk, and S. D. Huber, Phys. Rev. B **99**, 020304(R) (2019).

[21] B. A. Bernevig and T. L. Hughes, *Topological Insulators and Topological Superconductors* (Princeton University Press, Princeton, 2013).

[22] L. Fidkowski, Phys. Rev. Lett. **104**, 130502 (2010).

[23] H. Li and F. D. M. Haldane, Phys. Rev. Lett. **101**, 010504 (2008).

[24] Q. Wang, D. Wang, and Q.-H. Wang, Europhys. Lett. **124**, 50005 (2018).

[25] P. Zhu, K. Loehr, and T. L. Hughes, arXiv:1910.10180.

[26] R. Takahashi, Y. Tanaka, and S. Murakami, Phys. Rev. Research **2**, 013300 (2020).

[27] J. I. Cirac and F. Verstraete, J. Phys. A: Math. Theor. **42**, 504004 (2009).

[28] J. I. Cirac, D. Poilblanc, N. Schuch, and F. Verstraete, Phys. Rev. B **83**, 245134 (2011).

[29] M. B. Hastings, Phys. Rev. B **73**, 085115 (2006).

[30] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Phys. Rev. Lett. **100**, 070502 (2008).

[31] N. Schuch, D. Pérez-García, and I. Cirac, Phys. Rev. B **84**, 165139 (2011).

[32] N. Schuch, I. Cirac, and D. Pérez-García, Ann. Phys. **325**, 2153 (2010).

[33] C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A **81**, 052338 (2010).

[34] S. Yang, T. B. Wahl, H.-H. Tu, N. Schuch, and J. I. Cirac, Phys. Rev. Lett. **114**, 106803 (2015).

[35] A. Kitaev, Ann. Phys. **303**, 2 (2003).

[36] M. A. Levin and X.-G. Wen, Phys. Rev. B **71**, 045110 (2005).

[37] F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066.

[38] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Phys. Rev. Lett. **96**, 220601 (2006).

[39] Z.-C. Gu, M. Levin, B. Swingle, and X.-G. Wen, Phys. Rev. B **79**, 085118 (2009).

[40] T. B. Wahl, H.-H. Tu, N. Schuch, and J. I. Cirac, Phys. Rev. Lett. **111**, 236805 (2013).

[41] J. Dubail and N. Read, Phys. Rev. B **92**, 205307 (2015).

[42] N. Read, Phys. Rev. B **95**, 115309 (2017).

[43] M. P. Zaletel and R. S. K. Mong, Phys. Rev. B **86**, 245305 (2012).

[44] D. J. Thouless, Phys. Rev. B **27**, 6083 (1983).

[45] A. J. Károly, *A Short Course on Topological Insulators: Band-structure Topology and Edge States in One and Two Dimensions* (Springer, Berlin, 2016).

[46] M. Fannes, B. Nachtergaele, and R. F. Werner, Commun. Math. Phys. **144**, 443 (1992).

[47] T. B. Wahl, S. T. Haßler, H.-H. Tu, J. I. Cirac, and N. Schuch, Phys. Rev. B **90**, 115133 (2014).

[48] A. Deitmar and S. Echterhoff, *Principles of Harmonic Analysis* (Springer, Berlin, 2014).

[49] N. Bultinck, D. J. Williamson, J. Haegeman, and F. Verstraete, Phys. Rev. B **95**, 075108 (2017).

[50] D. Pérez-García, M. M. Wolf, M. Sanz, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. **100**, 167202 (2008).

[51] S. Bravyi, Quantum Inf. Comput. **5**, 216 (2005).

[52] I. Peschel, J. Phys. A: Math. Gen. **36**, L205 (2003).

[53] N. Schuch and B. Bauer, Phys. Rev. B **100**, 245121 (2019).

[54] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Inf. Comput. **7**, 401 (2007).

[55] M. T. Fishman and S. R. White, Phys. Rev. B **92**, 075132 (2015).

[56] B. J. Wieder, Z. Wang, J. Cano, X. Dai, L. M. Schoop, B. Bradlyn, and B. A. Bernevig, Strong and “fragile” topological dirac semimetals with higher-order fermi arcs, Nat. Commun. **11**, 627 (2020).

[57] N. Bultinck, D. J. Williamson, J. Haegeman, and F. Verstraete, J. Phys. A: Math. Theor. **51**, 025202 (2017).

[58] Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. B **84**, 075128 (2011).

[59] Y. Zhang and A. Vishwanath, Phys. Rev. B **87**, 161113(R) (2013).

[60] S. A. Parameswaran, R. Roy, and S. L. Sondhi, C. R. Phys. **14**, 816 (2013), topological insulators / Isolants topologiques.

[61] Y.-H. Wu, L. Wang, and H.-H. Tu, arXiv:1910.11011.

[62] I. Peschel and V. Eisler, J. Phys. A: Math. Theor. **42**, 504003 (2009).

[63] A. Kitaev, Ann. Phys. **321**, 2 (2006).

[64] A. Botero and B. Reznik, Phys. Rev. A **67**, 052311 (2003).

[65] G. Giedke and J. Ignacio Cirac, Phys. Rev. A **66**, 032316 (2002).

[66] J. Eisert, S. Scheel, and M. B. Plenio, Phys. Rev. Lett. **89**, 137903 (2002).

[67] J. Fiurášek, Phys. Rev. Lett. **89**, 137904 (2002).

[68] E. Zohar, M. Burrello, T. B. Wahl, and J. I. Cirac, Ann. Phys. **363**, 385 (2015).

[69] A. Alexandrinata, T. L. Hughes, and B. A. Bernevig, Phys. Rev. B **84**, 195103 (2011).