Stabilizing multiple topological fermions on a quantum computer

In classical and single-particle settings, non-trivial band topology always gives rise to robust boundary modes. For quantum many-body systems, however, multiple topological fermions are not always able to coexist, since Pauli exclusion prevents additional fermions from occupying the limited number of available topological modes. In this work, we show, through IBM quantum computers, how one can robustly stabilize more fermions than the number of topological modes through specially designed 2-fermion interactions. Our demonstration hinges on the realization of BDI- and D-class topological Hamiltonians on transmon-based quantum hardware, and relied on a tensor network-aided circuit recompilation approach. We also achieved the full reconstruction of multiple-fermion topological band structures through iterative quantum phase estimation (IQPE). All in all, our work showcases how advances in quantum algorithm implementation enable noisy intermediate-scale quantum (NISQ) devices to be exploited for topological stabilization beyond the context of single-particle topological invariants.


INTRODUCTION
Many-body quantum effects like particle statistics and Hubbard interactions underscore some of the most exciting condensed matter phenomena, such as superconductivity and fractionalization [1][2][3] . Their interplay with single-particle properties like band topology is particularly fascinating, with rich fractionalized quasiparticles emerging out of Coulomb repulsion within dispersionless Landau levels, for instance. But unlike in purely singlebody settings, particle statistics often play a central role, determining all Fermi surface properties like electrical conductivity. In topological systems, Pauli exclusion also implies that topological robustness is only conferred upon the few electrons occupying the limited number of topological modes 4 .
Compared to single-body topological phenomena realizable in classical metamaterials, novel many-body phases are traditionally much harder to engineer and physically demonstrate. This is due to challenges in accessing and manipulating intrinsically fragile quantum states, unlike the macroscopic classical degrees of freedom of photonic, acoustic and electrical circuit lattices [5][6][7][8] . Fortunately, there is considerable recent experimental progress in various quantum systems such as ultracold atomic lattices, 9 , photonic systems 10 , silicon 11 and trapped ion systems 12 , which make Richard Feynman's vision 13 of utilizing quantum systems to simulate quantum Hamiltonians a reality.
Of particular versatility are universal quantum simulators, also known as quantum computers, which allow arbitrary quantum systems to be simulated, thereby in principle enabling any quantum phenomenon to be physically realized. Through the appropriate design of sequences of quantum operations, known collectively as quantum algorithms, quantum computers are capable of performing a vast array of computational tasks, some at polynomial or exponential resource advantage over classical algorithms. Indeed, even in the current noisy intermediate-scale quantum (NISQ) era, quantum computers have already shown great promise, with demonstrations of quantum advantage in limited settings 14,15 , mapping topology in parameter space 16 , the achievement of chemical accuracy in intermediate-scale electronic structure calculations 17 , and in neutron scattering and exotic magnetic phenomena simulations 18,19 .
This work shall employ computations made with IBM quantum computers, which not only rank favourably with other cuttingedge platforms 20,21 in hardware performance, such as gate fidelities and decoherence times, but is also fully accessible via the cloud. IBM Quantum (IBM Q) currently offers access to up to 65-qubit machines based on superconducting transmon qubits, and provides an open-source software development kit called Qiskit 22,23 . Thus far, IBM Q has been successfully utilized in simulating spin models 24 , global quantum quenches 25 , quantum chemistry problems 26 , topological phenomena [27][28][29][30] , machine learning 31 and various other applications 32 .
We emphasize that quantum computation in the current NISQera is still plagued with significant limitations. Main bottlenecks include low gate fidelity, decoherence, limited qubit connectivity and limited number of qubits 33 , which together impose constraints on circuit depth and structure. Amidst qubit noise and readout error, typical depths of Oð10 1 Þ entangling gate layers are presently feasible for precision results. Of utmost current priority is hence the development of error mitigation and circuit optimization approaches that maximize computational capability within hardware bounds, thereby enabling the practical use of quantum simulators in contexts where classical simulators, for example, purpose-designed electrical circuits 34 , are inadmissible.
In this work, we stabilize BDI-and D-class topological boundary states, as well as demonstrate a full band structure reconstruction of the fermionic extended Kitaev chain (KC) on IBM quantum computers. Compared to existing quantum computer realizations of other topological states [27][28][29] , some of which are performed in parameter space, ours was performed on a longer (12-qubit) chain with physical open boundaries that host topological modes. Furthermore, our extended KC Hamiltonian contains multiple 1 non-local couplings which presented heavy demands on circuit depth and complexity, necessitating the use of tensor networkaided circuit recompilation techniques beyond traditional trotterization. Crucially, by exploiting the quantum nature of the IBM machines, we engineered effective interactions to stabilize more fermions than originally allowed by the number of topological zero modes, hence physically realizing a few-body phenomenon that has not been possible in existing classical realizations.

RESULTS AND DISCUSSION Physical motivation and models
Topological robustness is a highly sought-after property exemplified by the extraordinarily long survival duration of boundary states in specially designed topological lattices 35,36 . This robustness has inspired various potential applications like topological lasers and sensors 37,38 , and originates from the integer quantization of topological invariants characterizing the lattice band structure 4 . While non-trivial topology has been demonstrated in a wide variety of classical metamaterials, true quantization of the response can only be observed in quantum settings (but see refs. 39,40 ).
Yet, a fully quantum topological fermionic system is also limited by the fact that there can only be as many topological fermions as available topological modes. Since the latter is determined by the topological invariant, which is typically an integer of order Oð10 0 Þ, it will be of great scientific and practical interest to probe how quantum interactions can also enhance the number of robustly surviving fermions beyond what is dictated by the topological invariant. We shall first introduce our topological lattice models, and later discuss how specific interactions preserve the fidelity of multiple fermions by hosting many-body states that are adiabatically connected to non-interacting topological states.
In this work we simulate on a quantum computer the extended Kitaev model 41,42 , which is a 1D open chain containing nextnearest neighbour (NNN) couplings in addition to its underlying nearest-neighbour (NN) structure, with two sites per unit cell: þ h:c:; (1) for N unit cells, chemical potential μ, tunnelling coefficients v 1,2 , superconducting pairing Δ 1,2 , relative phase ϕ, and fermionic operators c j , d j respectively acting on the A-and B-sublattices of the j-th unit cell. Our model H KC is an extension of the standard Kitaev chain model that has been intensely studied for hosting Majorana zero modes 43 , which is of particular interest to faulttolerant topological quantum computing 44 . Its extra degrees of freedom in Eq. (1) allow for smooth tuning across the BDI-and D-symmetry protected topological (SPT) classes elaborated below, and will ultimately also be useful in the design of interactions that preserve the robustness of multiple fermions. When ϕ = 0, H KC preserves time-reversal symmetry (TRS) in addition to parity and charge conjugation symmetry, and belongs to the BDI class of the tenfold-way topological classification 45 , characterized by a winding number ν 2 Z. Our particular model possesses ν = 0, 1, 2 regimes, the first trivial and the latter two respectively exhibiting twofold and fourfold degeneracy of midgap topological modes localized at either boundary. Mathematically, 2πν corresponds to the winding of its Berry phase across one Brillouin zone period (see Supplementary Note 1).
When ϕ ≠ 0, TRS is broken and H KC falls into the D-class characterized by Z 2 invariant, with its topologically trivial/non-trivial phases respectively labelled by Berry phase winding γ = 0, π which encapsulates the relative configurations of the k and −k paths 41,42 . This D-class phase minimally requires NNN couplings, not present in the BDI class which contains the extremely well-known Su-Schrieffer-Heeger (SSH) model 46,47 , and is thus much less frequently investigated let alone realized in quantum settings 48 . We note that the SSH model, which we shall also simulate, is a special case of the extended Kitaev model upon unitary rotation: with v, w intra-cell and inter-cell hopping coefficients. We remark that H SSH and H KC are quadratic in the fermionic basis and are, in principle, efficiently simulable classically 49 ; but with the addition of quartic interactions, both models cannot be efficiently classically simulated.

Persistent boundary modes from topological protection
As the first experiments on the IBM quantum computer, we demonstrate that initial states at the boundary survive much longer when the Hamiltonian is topological. We emphasize that this is true for a rather large parameter space of initial states, not just exact topological eigenstates (see Supplementary Note 4). For a start, various perfectly localized 1-fermion initial states defined in Table 1 are evolved via H SSH and H KC . We utilize a tensor network-aided recompilation technique (see Methods) to construct quantum circuits for time-evolution, in order to overcome circuit depth limitations on NISQ hardware. Our approach is based on prior literature [50][51][52] , with added sectorspecificity to enhance performance at larger fermion numbers and qubit counts. Furthermore, we employ a suite of error mitigation

Model
Regime of stability The list includes idealized 1-particle edge states localized at certain sites near the boundary, as well as subscripted 'M' states at the middle of the chain, used for comparison purposes. Some of these states are adiabatically connected to topological eigenstates, as elaborated in Supplementary Note 1. Multi-fermion idealized states are constructed from the respective fermionic creation operators c[α] † -for example 2 LR techniques to improve data quality, in particular readout mitigation (RO) 25,53 , post-selection (PS) 25,54 , and averaging across machines and qubit chains (see Methods). Through computational-basis measurements on simulation qubits, the along the chain is retrieved, and the extent of evolution away from ψ 0 j i, whose occupancy is O 0 , can be assessed via the occupancy fidelity We present 1-particle time-evolution results for H SSH in Fig. 1a, comparing raw data, data with RO and PS mitigation, and data with both. The effectiveness of the error mitigation methods is clear, with the occupancy fidelity of all initial states approaching close to exact diagonalization (ED) results with RO and PS applied. The stability, i.e. persistence, of the initial state is quantified via the occupancy fidelity F O ; the slower the decay in F O , the more stable the state. As expected, the decay of the initial states 1 L , 1 R localized at the left and right boundary sites are slow in the topological regime, compared to that of the middle-localized state 1 M , which overlaps negligibly with topological boundary modes; by contrast, in the trivial regime topological protection does not exist and all states decay quickly. We remark that edge-localized Fig. 1 Topologically robust edge states and band structures simulated on quantum computers. a 1-particle time-evolution results for the SSH model H SSH , with well-preserved occupancy fidelity for perfectly boundary-localized 1 L;R states when ν = 1 (topological), and rapidly decaying occupancy fidelity for ν = 0 (trivial) cases and the middle-localized 1 M state. While raw hardware data already agrees qualitatively with exact diagonalization (ED) results, the agreement becomes excellent with RO and PS error mitigation. b Exact topological state wavefunctions at the left boundary, which constitute the dominant components of our 1-site, 2-site and 4-site boundary-localized initial states detailed in Table 1. c 1-particle time-evolved occupancy fidelity and spatial state density for BDI-class H KC in topological (ν = 2) and trivial (ν = 0) phases; and D-class H KC with ϕ = π/4 in topological (γ = 0) and trivial phases. Again, boundary localized states persist much longer when the system is topological, despite not being exact eigenstates. Middle initial states recur periodically due to boundary reflection from finite system length. d IQPE-obtained single-particle band structures for H SSH and H KC , with topological phase transitions marked by changes in 2ν, the number of degenerate midgap (E = 0) states. For H KC , parameters are held at v 1 = v 2 = Δ 1 = Δ 2 ≡ v, with ϕ = 0 for BDI-class and ϕ = π/4 for D-class. Time-evolution results are computed on chains with N = 6 unit cells (n = 12 qubits), with error bars obtained from repeated experiments across different qubit chains and machines (ibmq_paris, ibmq_toronto, ibmq_manhattan, and ibmq_boeblingen); band structure results are on n = 10 chains with an additional ancilla qubit, with error bars smaller than plot markers. Detailed parameter specifications are given in Table 2.
states can be stable even with significant perturbations, which is a hallmark of topological robustness. For instance, the stability of 1 L , in the sense of slow decay of F O , is preserved even when mixed with non-SPT states localized on the neighbouring and next-neighbouring sites (see Supplementary Note 4).
For H KC with NNN couplings, even relatively delocalized initial states can exhibit slow decay. We determine 'idealized' maximally localized boundary modes that are adiabatically connected to topological eigenstates (Supplementary Note 1), that exhibits negligible decay, such as to facilitate the design of more general persistent states. In the BDI class, they are 1 L 1 and 1 L 2 as defined in Table 1, which are localized on the first and second unit cells respectively. In Fig. 1c ) are robust only for the ν = 2 sector, which supports two topological boundary modes at each end. For D-class H KC systems with nonzero ϕ and NNN couplings, the four idealized boundary states are each localized on the four boundary sites with ϕ-dependent relative phases between site orbitals, as defined in Table 1. Indeed, almost negligible decay was observed (Fig. 1c) , which stays near unity as long as α j j 2 β j j 2 ( 1, that is, when a dominant eigenstate component exists. This can be checked by comparing all our initial states with exact topological eigenstates (Fig. 1b). In realistic settings where the initial boundary state generically overlaps with arbitrarily many eigenstates, the energy gap between the dominant topological eigenstate and all other states is crucial to stability. It introduces a separation of frequency scales that avoids overwhelming destructive interference. Experimentally, this energy gap can be verified from a full reconstruction of the topological band structure. Below, we first describe our approach for mapping the band structure, and then present its results for both single-fermion and two-fermion systems.
To probe and reconstruct the band structure on quantum hardware, we perform iterative quantum phase estimation (IQPE) 55,56 , which supports the estimation of eigenenergies of a simulated Hamiltonian, in principle to arbitrary precision (see Methods). Compared to quantum phase estimation (QPE) 57,58 , IQPE circuits are shallower and require fewer qubits, and are thus suited for implementation on current-generation NISQ hardware. IQPE results for H SSH and H KC in the 1-particle sector are shown in Fig. 1d. The BDI topological phase transition in H SSH at w = v is apparent, with a pair of midgap states separating from the bulk and approaching degeneracy at E = 0 as w/v increases. The H KC model possesses richer behaviour-in the BDI case, transitions from the ν = 1 phase into the trivial (ν = 0) and then into ν = 2 phase occur as v/μ increases, for illustrative v ≡ v 1 = v 2 = Δ 1 = Δ 2 . The ν = 1 phase exhibits twofold degeneracy of midgap states like in the SSH model, and the ν = 2 phase exhibits fourfold degeneracy at E = 0 due to having two zero modes at each boundary. In the D class, twofold midgap topological degeneracy occurs in the γ = π but not γ = 0 phase. For the latter, fourfold degeneracy is broken as the eigenenergies split into pairs of positive energy and negative energy states, which can be regarded as remnants of topological midgap states. In all cases, the reconstructed band structure from hardware execution of IQPE closely match ED results.
Topological stability for multiple fermions-non-interacting case Compared to existing classical realizations of BDI-and D-class topological states, our IBM Q realizations possess the advantage of granting natural access to quantum many-body effects like fermionic statistics and interactions. We shall first discuss the former, specifically on the various ways whereby up to four fermions can simultaneously enjoy topological robustness.
We first consider H SSH and the two-fermion sector, in which twofermion boundary states can be constructed either as 2 LR ¼ c½1 L y c½1 R y vac j i with 1 fermion at each boundary, or alternatively as 2 L AA or 2 L AB with one fermion in the 1 L orbital and the other in the nearest A or B site. Of these, only 2 LR has both particles overlapping significantly with topological boundary modes, so only it would be conferred stability in the topological phase (Fig. 2b). The 2 L AA and 2 L AB states, like 1 M , are unstable whether topological modes exist or not, despite one of the fermions being in the topologically stable orbital 1 L . The upshot is that since H SSH possesses only one topological state at each end of the chain, it is impossible to construct a 2-particle topological state localized at only one boundary. Generalization to multiple fermions follows.
To realize 2-fermion topological states localized at a single boundary, the fourfold midgap topological modes of H KC , with two at each boundary, can be utilized. In particular, the two fermions can occupy 1 L 1 and 1 L 2 near the left boundary, resulting in 2 L

12
, or analogously 2 R 12 for the right boundary. Relaxing the requirement of localization on a single edge, robust states such as j2 LR ij i ¼ c½1 L i y c½1 R j y jvaci are all possible, where i, j ∈ {1, 2}. Indeed, hardware results indicate these states are conferred stability in the topological regime; the contrast against the non topologically protected middle-localized 2 M is drastic (Fig. 2c). The occupancy density maps concur that these 2-fermion edge localized states persist almost perfectly in the topological phase, but diffuses in the trivial phase.
As a generalization to larger numbers of effective particles, we demonstrate that up to 4 fermions can enjoy stability in the D-class H KC system, with each occupying a different topological mode. Indeed, the 4-fermion state 4 LLRR is robust compared to the middle-localized 4 M state in the topological phase (Fig. 2c).

Stability for multiple fermions-interacting case
Quantum interactions can present new phase transitions 59 and avenues of stability with no non-interacting analogues 60,61 . The understanding of strongly correlated topological models is also crucial in understanding the phenomenology of real materials 62,63 . While weak interactions may only perturb the stability of SPT boundary modes, strong interactions can drastically break the symmetry protection altogether, leading to new preferred states. Our time-evolution and IQPE methods, hinging on circuit recompilation, readily supports the study of strongly correlated models, facilitating digital quantum computers as alternative experimental platforms from existing ultracold atomic lattices 64 or photonic crystals 65 . We consider two-body Hubbard interactions which have been commonly used to model strong density-density interactions 64,66 . Among the simplest 1D topological model with interactions is the SSH-Hubbard model, which has been explored in various contexts 67,68 . We are interested in Hubbard interactions that directly compete with the boundarylocalization of SSH topological modes. Specifically, we consider interactions between all successive sites (U 2 ), sites within the same unit cell (U 4 ), and same sublattice sites across adjacent unit cells (U 5 ). Our interacting SSH model is then obtained by adding these interaction terms to Eq. (2), where c y j (d y j ) creates a fermion in the A (B) sublattice of the j-th unit cell. We emphasize that these interactions are homogeneous across all unit cells, and do not induce boundary effects on their own, without the interplay with topological boundary localization. The schematic in Fig. 3a illustrates the interaction terms in Eq. (4). The U 4 and U 5 terms are intentionally chosen to induce asymmetry across the two sublattices-they impact states with inhomogeneous polarizations differently, such as SSH topological modes which occupy one sublattice exclusively. We dynamically evolve on IBM Q various boundary-localized 2-fermion states 2 LR , 2 L AB , 2 L AA , 2 R AB , 2 R BB , schematically shown in Fig. 3b. As before, the superscripts indicate the edge localization (left/right) of the 2 fermions, while the subscripts indicate the occupied sublattices; we remind the reader that the chain starts with an A-site on the left. Our chosen 2-particle interactions have enabled us to achieve stable multi-fermion boundary states due to either topological localization or interaction-induced polarization, and often the interplay of both. As such, we can for instance produce stable states localized on a single boundary in H SSH full , reminiscent of topological states, but yet existing even in the topologically trivial regime of the non-interacting part H SSH . The conferred stability generally increases with interaction strength at diminishing returns; beyond U/w~10, little additional stability is gained from stronger interactions. Yet, the impact of the topological character of the non-interacting terms in our systems remains significant even in the face of arbitrarily strong interactions.
In Fig. 3c-d, we present the evolution of various boundarylocalized 2-fermion initial states subject to strong interactions U 2 , U 4 or U 5 of relative strength 10. While such strong interactions may seem clearly dominant compared to the non-interacting H SSH part, the topological character of H SSH still affects the evolution of most of these states significantly. When U 2 is switched on (Fig. 3c) , not shown) retain near-perfect fidelity, even when the chain is topologically trivial. For the 2 LR initial state localized at both boundaries, topology enables a longer survival time. To obtain stable 2-fermion states localized at a single boundary, we may switch on U 4 and U 5 ; when U 5 > 0, both types of left boundary modes ( 2 L AB and 2 L AA ) are stabilized. On the other hand, to obtain stable right boundary states, we require U 4 > 0 as well, such as to balance with U 5 . This is reflected from the stability plot in parameter space (Fig. 3e)-the time-averaged occupancy fidelity F O for the left-localized 2-particle boundary states is minimal along the U 4 = U 5 diagonal, being greatly destabilized relative to the rightlocalized modes. The preferential stabilization at either boundary by U 4 and U 5 is a consequence of their coupling of specific sublattice pairs on finite chains; when the interaction couples all sublattices like for U 2 , two-fermion boundary modes are stable on both edges, and the preferential stability is not observed.
Earlier, we showed that the stability of initial boundary states can be maintained even when they differ rather significantly from topological ansatz states (Supplementary Note 4). We now extend the same remark to boundary states that are not topologically protected, but which are conferred robustness under the combination of topology and interaction-induced effects. In Fig. 3f, we mix 2 L AA with 2 L BB perturbation, neither of which are SPT, and examine their evolution under varying U 2 interactions. When U 2 = 0, as expected, F O is relatively low; but for U 2 > 0 an increase in F O is observed, reaching ≳ 90% for U 2 ≳ 5/2.  Table 2.
This stabilization holds across a wide range of mixing amplitudes. Most interestingly, the stabilization occurs only in the topological phase of H SSH -indeed in the trivial regime, U 2 has virtually no effect on state stability. Hence, neither topology nor interactions alone suffices to enable the observed perturbation-robust stabilization of non-SPT states.
Similarly, one can study Hubbard interactions on the extended Kitaev model. We again focus on interaction terms that induce sublattice asymmetry, and thus left/right boundary asymmetry for topological modes. The extended Kitaev chain (Eq. (1)), however, has longer-range couplings compared to the SSH model (Eq. (2)). The NNN hoppings and pairings were crucial in the noninteracting model to demonstrate a richer topology, exhibiting more than one pair of topological modes in the ν = 2 phase of the BDI class. We thus consider interactions that compete with this topology, specifically between occupied A-sites of NN (U 2 ) an NNN (U 3 ) unit cells, with These interactions are schematically illustrated in Fig. 3a. In the BDI class of the non-interacting part H KC , we dynamically evolve the various 2-fermion states used earlier, 2 L AA and 2 L BB , which now have large overlaps with the topological modes of the noninteracting Kitaev chain H KC (see Fig. 3b for schematics). Recall that the BDI-class of H KC can host up to two pairs of topological modes at the two boundaries. As shown in Fig. 4a, the presence of either interaction (U 2 or U 3 ) worsens the fidelity of all 2-particle boundary modes when ν = 2. These interactions hence disrupt the symmetry protection conferred in the BDI class. When ν = 1, there is one fewer pair of topological boundary modes, but the states 2 L A0A and 2 R B0B -which do not overlap significantly with any non-interacting topological mode-now enjoy enhanced robustness in the presence of interactions. Also, all edge-localized modes can be made robust in the topologically trivial case ν = 0.
In the D class of H KC , we demonstrate that U 3 interactions can confer stability to 3-and 4-fermion states (Fig. 4b), which carry a greater number of fermions on a boundary than available topological modes, and are hence topologically unprotected. The stability of SPT 2 L 12 is not destroyed when U 3 is imposed, and simultaneously the non-SPT 3 LLL and 4 LLLR are stabilized, in both the γ = 0 topological and trivial phases of H KC , the latter more pronounced. Though not shown in the main text, U 2 has a similar but weaker effect (see Supplementary Note 4). We are thus able to engineer, in both H SSH and (both BDI-and D-classes of) H KC , enhanced stability of edge modes that goes beyond what can be achieved from standard topology alone. These results demonstrate how NNN Hubbard interactions can lead to avenues of interesting dynamical behaviour, which is not wholly surprising, given that they are known to lead to new phases like topological Mott phases in other contexts imbued with appropriate lattice structure, for example, the honeycomb lattice 69 .  Table 2 for detailed parameter specifications.
Finally, we present the full two-fermion band structure of the interacting SSH model (Eq. (4)) reconstructed on quantum hardware using IQPE (Fig. 3g). Similar results are obtained against U 5 instead of U 2 (Supplementary Note 4). With strong interactions, an additional band of energies whose large gap scales linearly with interaction strength appears. However, the interpretation of this additional set of bands is not straightforward. The interacting Hamiltonian does not share the same eigenstates as its non-interacting counterpart; in fact differences in eigenstate wavefunctions are drastic even for modest U~v, w comparable to the non-interacting H SSH . One therefore cannot hope to understand the separation of the band purely using the non-interacting eigenstates. Moreover, ED reveals that the band does not necessarily contain 2-fermion boundarylocalized modes, even though initial 2-fermion boundarylocalized states can be stabilized in some cases. The near-perfect fidelity of these stabilized 2-particle boundary modes, some representatives having been previously demonstrated (Fig. 3), is due to strong overlap with one of the many eigenstates ψ j i of the interacting Hamiltonian; but the location of ψ j i in the band structure may be in either set of bands.

Discussion and outlook
By realizing the interacting extended Kitaev Chain on IBM quantum computers, we have demonstrated how various boundary states can be robustly preserved by topological mechanisms of BDI-and D-class symmetries. Importantly, this topological protection is not limited to topological eigenstates, and interplays non-trivially with 2-body interactions and Pauli exclusion statistics in multi-fermion settings. Most spectacularly, we discovered avenues where interactions allow more fermions to be stabilized at one boundary than suggested by the number of available topological modes alone. Our work also illustrates how tensor network-aided circuit recompilation techniques beyond traditional trotterization enable the simulation and full band structure reconstruction of complex topological Hamiltonians on quantum circuits. Our approach can be further extended to more sophisticated many-body interacting Hamiltonians, presenting new opportunities and raising the state-of-the-art in the quantum simulation of strongly correlated topological systems.

Quantum simulation of state evolution
Given an initial state ψ 0 j i and a time-independent Hamiltonian H ∈ {H SSH , H KC }, the time-evolved state is ψðtÞ j i¼UðtÞ ψ 0 j i, with propagator U(t) = e −iHt . A schematic of the quantum circuit performing time-evolution is given in Fig. 5a. Traditional implementation of U(t) on quantum circuits entails expanding the Hamiltonian H in the spin-1/2 basis and employing trotterization 70,71 ; but acceptable trotterization error requires numerous time steps, yielding infeasibly deep circuits. Furthermore, each trotter step, comprising terms e Àiβ σ σΔt for Pauli strings σ and coefficients β σ (see Supplementary Note 3), requires layers of entangling gates scaling with the weight of σ, thus impairing its usefulness for models with longer-range couplings.
To transcend these limitations, we employ an implementation strategy known as circuit recompilation, built upon prior studies [50][51][52] . A circuit Table 2. Summary of parameter specification used in all figures.

Figure
Regime Parameters This regime contains eigenstates that are adiabatically connected to the ν = 0 phase as we approach BDI class (ϕ → 0).  Table 2 for detailed parameter specifications. ansatz (Fig. 5c) comprising n L ≤ 8 repetitive layers of general U 3 rotation and CX entangling gates is iteratively optimized, through tensor networkbased quantum simulation 72 , to approach the intended unitary. Specifically, the U 3 angles χ = (θ, ϕ, λ) are fine-tuned with L-BFGS-B with basinhopping (see Supplementary Note 3). We design the circuit ansatz to contain no longer-range CXs, in order to conform to qubit connectivity on hardware, and the CXs are placed in a brickwork pattern to maximize entangling power, to accommodate the intended unitary in as few layers as possible. This technique yields order-of-magnitude shallower circuits than trotterization at comparable error rates, and is critical in our acquisition of high-quality experiment data on NISQ hardware. For larger fermion numbers (≥3), we improve the optimization procedure to focus on the relevant sectors, referred to as sector-specific recompilation, so as to overcome the impairment in performance brought about by the increased entanglement in evolved states.

KC BDI-Class
We employ error mitigation techniques to improve the quality of hardware results-first, readout error mitigation (RO) reverses bit-flips during measurement based on prior calibration of qubits 25,53 ; second, post-selection (PS) is performed on particle number 25,54 . Indeed, since H SSH full and H KC full are number-conserving, time-evolved states that fall outside the particle number sector of ψ 0 j i are unphysical and can be discarded; no additional circuit depth is incurred since the measurements for O suffice. To feasibly accommodate the number of qubits n used in our simulations (7 ≤ n ≤ 13), we adopt a tensored RO scheme that differs from the open-source implementation in Qiskit 22 . See Supplementary Note 2 for an introductory overview of quantum computation, and Supplementary Note 3 for technical details of our abovementioned quantum simulation methods.
We end the section with a conceptual remark pertaining to the mapping between the fermionic models and the spin-1/2 qubits on the quantum computers. In traditional trotterization, the Jordan-Wigner or Bravyi-Kitaev transforms 70,71 map the fermionic Hamiltonian into spin chains; the constituent Pauli strings are then naturally expressible on quantum circuits. In comparison, using circuit recompilation, instead of assembling quantum circuits from a fixed set of formulae, an approximate mapping between the fermionic propagator and spin-1/2 quantum gates is dynamically determined during optimization.

Iterative quantum phase estimation
Given U ψ j i ¼ e 2πiϕ ψ j i for unitary U and eigenstate ψ j i, IQPE estimates the eigenphase ϕ ∈ [0, 1), in principle to arbitrary precision. Setting U = e −iHt allows the inference of eigenenergy E = −2πϕ/t of ψ j i. As mentioned, compared to QPE 57,58 , IQPE circuits are shallower and require fewer qubits -only a single ancilla qubit and a single controlled-unitary block is required. There is no need for multi-qubit inverse Fourier transforms. Truncating the binary expansion ϕ = 0. ϕ 1 ϕ 2 …ϕ m to m bits, IQPE iterates from k = m to k = 1; in iteration k, a controlled-U 2 kÀ1 block and a feedback R z (ω k ) = −2π(0.0ϕ k+1 ϕ k+2 …ϕ m ) rotation are applied, and the ancilla qubit is measured to determine ϕ k . An IQPE circuit diagram is shown in Fig. 5b.
To minimize circuit depth, the initialization of ψ j i and the controlledunitary block are implemented via recompilation; RO mitigation is applied to all qubits, and PS is applied to the simulation qubits to select for specific particle number sectors. Indeed, by performing IQPE over numerous ψ j iwhich need not be exact eigenstates, since superposition states collapse at measurement and yield expected eigenenergies nonetheless-the band structure of the Hamiltonian H can be recovered at arbitrarily high resolution. See Supplementary Note 3 for further implementation details.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.