Realization of higher-order topological lattices on a quantum computer

Programmable quantum simulators may one day outperform classical computers at certain tasks. But at present, the range of viable applications with noisy intermediate-scale quantum (NISQ) devices remains limited by gate errors and the number of high-quality qubits. Here, we develop an approach that places digital NISQ hardware as a versatile platform for simulating multi-dimensional condensed matter systems. Our method encodes a high-dimensional lattice in terms of many-body interactions on a reduced-dimension model, thereby taking full advantage of the exponentially large Hilbert space of the host quantum system. With circuit optimization and error mitigation techniques, we measured on IBM superconducting quantum processors the topological state dynamics and protected mid-gap spectra of higher-order topological lattices, in up to four dimensions, with high accuracy. Our projected resource requirements scale favorably with system size and lattice dimensionality compared to classical computation, suggesting a possible route to useful quantum advantage in the longer term.

In this work, we develop a new approach that establishes NISQ hardware as particularly suitable platforms for the simulation of generic multi-dimensional condensed-matter lattice models.As a demonstration, we utilize our method to realize, on transmon-based su-perconducting quantum devices, higher-order topological (HOT) phases in high-dimensional lattices of unprecedented size and complexity.Unlike previous quantum simulator studies that implemented topological models through synthetic dimensions [26][27][28][29][30][31][32][33][34][35][36], we realize HOT lattices in real space, in up to d = 4 dimensions (tesseract).Central to our approach is a mapping procedure that encodes single-particle degrees of freedom of a highdimensional lattice within the many-body Fock space of an interacting one-dimensional (1D) model.This enables us to take full advantage of the exponentially large Hilbert space innately accessible by a quantum computer, while drastically reducing the number of qubits needed for direct simulation.
We remark that classical simulation of highdimensional HOT lattices is expensive, and the resource complexity of our quantum simulation approach can be favorable over classical numerical methods-e.g.exact diagonalization (ED).Although the scalable realization of larger HOT systems requires hardware exceeding present capabilities, our approach could present an avenue toward quantum advantage as quantum simulator platforms continue to rapidly improve.
Our unique approach, especially, lends the reinterpretation of HOT robustness as an interactionmediated phenomenon in 1D, thereby introducing a new class of 1D many-body systems with topologically protected clustering properties.HOT, in particular, becomes an alternative mechanism for enforcing d-body clustering or repulsion [63][64][65][66] beyond the scope of known mechanisms such as Coulomb repulsion.Promoting these HOT-protected boundary states to spatial corners further opens avenues for applications, such as the realization of a topological qubit using HOT superconductors that host elusive Majorana corner modes [67][68][69].This complements the search for HOT phases in quantum materials [49,50,[70][71][72][73][74][75][76], which is presently still in its infancy.

II. RESULTS & DISCUSSION
A. Mapping higher-dimensional lattices to 1D quantum chains While small quasi-1D and 2D systems have been simulated on digital quantum computers [77,78], the explicit simulation of higher-dimensional lattices remains elusive.Directly simulating a d-dimensional lattice of width L along each dimension requires ∼L d qubits 1 .For large dimensionality d or lattice size L, this quickly becomes infeasible on NISQ devices, which are significantly limited by the number of usable qubits, qubit connectivity, gate errors, and decoherence times.
To overcome these hardware limitations, we devise an approach to exploit the exponentially large many-body Hilbert space of an interacting qubit chain.The key inspiration is that most local lattice models only access a small portion of the full Hilbert space (particularly noninteracting models and models with symmetries), and a L d -site lattice can be consistently represented with far fewer than L d qubits.To do so, we introduce an exact mapping that reduces d-dimensional lattices to 1D chains hosting d-particle interactions, which is naturally simulable on a quantum computer that accesses and operates on the many-body Hilbert space of a register of qubits. 1 The prefactor depends on spin degrees of freedom of the particles inhabiting the lattice.In particular, m = 1 qubits representing each site suffices for spinless particles, but m > 1 is required for spinful fermions and bosons.In real space, where we have associated the band degrees of freedom to a sublattice structure γ, and h γγ ′ rr ′ = 0 for |r − r ′ | outside the coupling range of the model, i.e. adjacent sites for a nearest-neighbor model, next-adjacent for next-nearestneighbor, etc.The operator c rγ annihilates particle excitations on sublattice γ of site r.
To take advantage of the degrees of freedom in the many-body Hilbert space, our mapping is defined such that the hopping of a single particle on the original ddimensional lattice from (r ′ , γ ′ ) to (r, γ) becomes the simultaneous hoppings of d particles, each of a distinct species, from locations (r ′ 1 , . . ., r ′ d ) to (r 1 , . . ., r d ) and sublattice γ ′ to γ on a one-dimensional interacting chain.Explicitly, this map is given by where r α is the α th component of r, and {ω α ℓγ } d α=1 represent d excitation species hosted on sublattice γ of site ℓ on the interacting chain, yielding In the single-particle context, exchange statistics is unimportant and {ω α } can be taken to be commuting.This mapping framework accommodates any lattice dimension and geometry, and any number of bands or sublattice degrees of freedom.As the mapping is performed at the second-quantized level, any one-body Hamiltonian expressed in second-quantized form can be treated, which encompasses a wide variety of single-body topological phenomena of interest.With slight modifications, this mapping can be extended to admit interaction terms in the original d-dimensional lattice Hamiltonian, although we do not explore them further in this work.

Mapping HOT lattices onto 1D qubit chains
For concreteness, we specialize our Hamiltonian to higher-order topological (HOT) systems henceforth, and shall detail how our mapping enables them to be encoded on quantum processors.The simplest square lattice with HOT corner modes [47,51,79] may be constructed from the paradigmatic one-dimensional Su-Schrieffer Heeger (SSH) model [55,80].To allow for sufficient degrees of freedom for topological localization, we minimally require a 2D mesh of two different types of SSH chains in each direction, arranged in an alternating fashion where c xy is the annihilation operator acting on site (x, y) of the lattice and u α r1r2 takes values of either v α r1r2 for intra-cell hopping (odd r 2 ) or w α r1r2 for inter-cell hopping (even r 2 ), α ∈ {x, y}.Conceptually, we recognize that the 2D lattice momentum space can be equivalently interpreted as the joint configuration momentum space of two particles, specifically, the (1 + 1)-body sector of a corresponding one-dimensional interacting chain.We map c xy → µ x ν y , where µ ℓ and ν ℓ annihilate hardcore bosons of two different species at site ℓ on the chain.In the notation of Eq. 2, we identify ω 1 ℓ = µ ℓ and ω 2 ℓ = ν ℓ , and the sublattice structure has been absorbed into the (parity of) spatial coordinates.This yields an effective one-dimensional, two-boson chain described by where n ω ℓ is the number operator for species ω at site ℓ of the chain.As written, each term in H 2D chain represents an effective SSH model for one particular species µ or ν, with the other species not participating in hopping but merely present (hence its number operator).These twobody interactions arising in H 2D chain appear convoluted, but can be readily accommodated on a quantum computer, taking advantage of the quantum nature of the platform.To realize H 2D chain on a quantum computer, we utilize 2 qubits to represent each site of the chain, associating the unoccupied, µ-occupied, ν-occupied and both µ, ν-occupied boson states to qubit states |00⟩, |01⟩, |10⟩, and |11⟩ respectively.Thus 2L qubits are needed for the simulation, a significant reduction from L 2 qubits without the mapping, especially for large lattice sizes.We present simulation results on IBM quantum computers for lattice size L ∼ O (10) in Section II C.
Our methodology naturally generalizes to higher dimensions.Specifically, a d-dimensional HOT lattice maps onto a d-species interacting one-dimensional chain, and d qubits are employed to represent each site of the chain, providing sufficient many-body degrees of freedom to encode the 2 d occupancy basis states of each site.We write where α enumerates the directions along which hoppings occur and êα is the unit vector along α.As before, the hopping coefficients alternate between interand intra-cell values that can be different in each direction.Compactly, for parity function π, and r α are spatial coordinates in non-α directions-see Supplementary where ω α ℓ annihilates a hardcore boson of species α at site ℓ of the chain and n α ℓ is the number operator of species α.In the d = 2 square lattice above, we had r = (x, y) and {ω α } = {µ, ν}.The highest dimensional HOT lattice we shall examine is the d = 4 tesseract, for which r = (x, y, z, w) and {ω α } = {µ, ν, η, ξ}.In total, a d-dimensional lattice Hamiltonian has d × 2 d distinct hopping coefficients, since there are d different lattice directions and 2 d−1 distinct edges along each direction, each comprising 2 distinct hopping amplitudes for inter-and intra-cell hopping.Appropriately tuning these coefficients allow the manifestation of robust HOT modes along the boundaries (corners, edges, etc.) of the lattices-see Figures 3-6 for schematics of the various lattice configurations investigated in our experiments.
Accordingly, the equivalent interacting 1D chain requires dL qubits to realize, an overwhelming reduction from the L d otherwise needed in a direct simulation of H dD lattice without the mapping.We remark that such a significant compression is possible because HOT is inherently a single-particle phenomenon.See Methods for further details and optimizations of our mapping scheme.

B. Simulation on quantum hardware
With our mapping, a d-dimensional HOT lattice H dD lattice with L d sites is mapped onto an interacting 1D chain H dD chain with dL number of qubits, which can be feasibly realized on existing NISQ devices for L ∼ O(10) and d ≤ 4. While the resultant interactions in H dD chain are inevitably complicated, below we describe how H dD chain can be viably simulated on quantum hardware.
A high-level overview of our general framework for simulating HOT time-evolution is illustrated in Figure 2. To evolve an initial state |ψ 0 ⟩, it is necessary to implement the unitary propagator U (t) = exp −iH dD chain t as a quantum circuit, such that the circuit yields |ψ(t)⟩ = U (t) |ψ 0 ⟩ and desired observables can be measured upon termination.A standard method to implement U (t) is trotterization, which decomposes H dD chain in the spin-1/2 basis and splits time-evolution into small steps (see Methods for details).However, while straightforward, such an approach yields deep circuits unsuitable for presentgeneration NISQ hardware.To compress the circuits, we utilize a tensor network-aided recompilation technique [81][82][83][84][85].We exploit the number-conserving symmetries of H dD chain , arising from H dD lattice and the nature of our mapping (see Methods), to enhance circuit construction performance and quality at large circuit breadths (up to 32 qubits).Moreover, to improve data quality amidst hardware noise, we employ a suite of error mitigation techniques, in particular readout mitigation (RO) that approximately corrects bit-flip errors during measurement [19,86], a post-selection (PS) technique that discards results in unphysical Fock-space sectors [19,82,87], and averaging across machines and qubit chains (see Methods).
After acting on |ψ 0 ⟩ by the quantum circuit that effects U (t), terminal computational-basis measurements are performed on the simulation qubits.We retrieve the site-resolved occupancy densities ρ(r) = ⟨c † r c r ⟩ = ⟨ d α=1 n α rα ⟩ on the d-dimensional lattice, and the extent of evolution of |ψ(t)⟩ away from |ψ 0 ⟩, whose occupancy densities are ρ 0 (r), is assessed via the occupancy fidelity Compared to the state fidelity F = |⟨ψ 0 |ψ⟩|2 , the occupancy fidelity F ρ is considerably more resource-efficient to measure on quantum hardware.
In addition to time-evolution, we can also directly probe the energy spectrum of our simulated Hamiltonian H dD chain through iterative quantum phase estimation (IQPE) [88,89]-see Methods.Specifically, to characterize the topology of HOT systems, we use IQPE to probe the existence of midgap HOT modes at exponentially suppressed (effectively zero for L ≫ 1) energies.In contrast to quantum phase estimation [90,91], IQPE circuits are shallower and require fewer qubits, and are thus preferable for implementation on NISQ hardware.As our interest is in HOT modes, we initiate IQPE with maximally localized boundary states that are easily constructed a priori, which exhibit good overlap (> 80% state fidelity) with HOT eigenstates, and examine whether IQPE converges consistently towards zero energy.These states are listed in Supplementary Table 2.

C. Two-dimensional HOT square lattice
As the lowest-dimensional incarnation of HOT lattices, the d = 2 staggered square lattice harbors only one type of HOT modes-zero-dimensional corner modes (Figure 1a).Previously, such HOT corner modes on 2D lattices have been realized in various metamaterials [92][93][94][95][96][97][98][99][100][101], but not in a quantum setting to-date.Our equivalent 1D hardcore boson chain is also the first demonstration of interaction-induced higher-order topology, where the topological localization is mediated not due to physical SSH-like couplings or band polarization, but due to the combined exclusion effects from all its interaction terms.We emphasize that our physically realized 1D chain contains highly non-trivial interaction terms involving multiple sites-the illustrative example in Figure 3f for an L = 6 chain already contains a multitude of interactions, even though it is much smaller than the L = 10 and L = 16 systems we simulated on quantum hardware.As evident, the 2 × 2 d = 8 unique types of interactions, corresponding to the 8 different couplings on the lattice, are mostly non-local; but this does not prohibit their implementation on quantum circuits.Indeed, the versatility of digital quantum simulators in realizing effectively arbitrary interactions allows the implementation of complex interacting Hamiltonian terms, and is critical in enabling our quantum device simulations.
In our experiments, we consider three different scenarios: C0, having no topological corner modes; C2, having two corner modes at corners (x, y) = (1, 1) and (L, 1); and C4, having corner modes on all four corners.These scenarios can be obtained by appropriately tuning the 8 coupling parameters in the Hamiltonian (Eq.4)-see Supplementary Table 1 for parameter values [102].
We first show that the correct degeneracy of midgap HOT modes can be measured on each of the configuration C0, C2 and C4 on IBM transmon-based quantum computers, as presented in Figure 3a.For a start, we used a 20-qubit chain, which logically encodes a 10 × 10 HOT lattice, with an additional ancillary qubit for IQPE readout.The number of topological corner modes in each case is accurately obtained through the degeneracy of midgap states of exponentially suppressed energy (red), as measured through IQPE executed on quantum hardwaresee Methods for details.That these midgap modes are indeed corner-localized is verified via numerical (classical) diagonalization, as in the insets of Figure 3a.
Next, we demonstrate highly accurate dynamical state evolution on larger 32-qubit chains on quantum hardware.We time-evolve various initial states on 16 × 16 HOT lattices in the C0, C2, and C4 configurations, and measure their final 2 site-resolved occupancy densities ρ(x, y).The resultant occupancy fidelity plots (Figures 3b-3e) conform to the expectation that states localized on topological corners survive the longest, and are also in excellent agreement with reference data from exact diagonalization (ED).For instance, a localized state at the corner (x 0 , y 0 ) = (1, 1) is robust on C2 and C4 lattice configurations (Figure 3b), whereas one localized on the (x 0 , y 0 ) = (1, L) corner is robust only on the C4 configuration (Figure 3c).These fidelity decay trends are corroborated with the measured site-resolved occupancy density ρ(x, y): low occupancy fidelity is always accompanied by a diffused ρ(x, y) away from the initial state, whereas strongly localized states have high occupancy fidelity.In general, the heavy overlap between an initial state and a HOT eigenstate confers topological robustness, resulting in significantly slowed decay; this is apparent from the occupancy fidelities F ρ , which remain near unity over time.In comparison, states that do not enjoy topological protection, such as the (1, L)-localized state on the C2 configuration and all initial states on the C0 configuration, rapidly delocalize and decay quickly.
Our experimental runs remain accurate even for initial states that are situated away from the lattice corners, such that they cannot enjoy full topological protection.
In Figure 3d, the initial state at (x 0 , y 0 ) = (2, 1), which neighbors the corner (1, 1), loses its fidelity much sooner than the corner initial state of Figure 3b, even for the C2 and C4 topological corner configurations.That said, its fidelity evolution still agrees well with ED reference data.In a similar vein, an initial state that is somewhat delocalized at a corner (Figure 3e) is still conferred a degree of stability when the corner is topological.

D. Three-dimensional HOT cubic lattice
Next, we extend our investigation to the staggered cubic lattice in 3D, which hosts third-order HOT corner modes (Figure 1a).These elusive corner modes have to date only been realized in classical platforms [103][104][105][106][107] or in synthetic electronic lattices [108].Compared to the 2D cases, the implementation of the 3D HOT lattice (Eq.6) as a 1D interacting chain (Eq.7) on quantum hardware is more sophisticated.The larger dimensionality of the staggered cubic lattice, in comparison to the square lat-tice, is reflected by a larger density of multi-site interaction terms on the interacting chain.This is illustrated in Figure 4b for the minimal 4 × 4 × 4 lattice, where the combination of the various d = 3-body interactions gives rise to emergent corner robustness (which appears as up to 3-body boundary clustering as seen on the 1D chain).
On quantum hardware, we implemented 18-qubit chains representing 6×6×6 cubic lattices in four configurations, specifically, the trivial lattice (C0), two geometrically inequivalent configurations hosting four topological corners (C4a, C4b), and a configuration with all 2 3 = 8 topological corners (C8).Similar to the 2D HOT lattice, we first present the degeneracy of zero-energy topological modes (header row of Figure 4a) with low-energy spectral data (red diamonds) accurately obtained via IQPE.
From the first row of Figure 4a, it is apparent that initial states localized on topological corners enjoy significant robustness.Namely, the measured site-resolved occupancy densities ρ(x, y, z) (four right columns) indicate that the localization of (x 0 , y 0 , z 0 ) = (1, 1, 1) corner initial states on C4a, C4b and C8 configurations are maintained, and measured occupancy fidelities remain near unity.In comparison, an initial corner-localized state on the C0 configuration, which hosts no topological corner modes, delocalizes quickly.Moving away from the corners, an edge-localized state adjacent to a topological corner is conferred slight, but nonetheless present, stability (second row of Figure 4a), as observed from the slower decay of the (x 0 , y 0 , z 0 ) = (2, 1, 1) state on C4a, C4b and C8 configurations in comparison to the C0 topologically trivial lattice.This conferred robustness is diminished for states localized further from topological corners, for instance surface-localized states (third row), and is virtually unnoticeable for states localized in the bulk (fourth row), which decay rapidly for all topological configurations.Initial states that are slightly delocalized near a corner enjoy some protection when the corner is topological, but is unstable when the corner is trivial (fifth row of Figure 4a).We again highlight the quantitative agreement of our quantum hardware simulation results with theoretical ED predictions.

E. Four-dimensional tesseract hyperlattice-HOT corner and edge modes
We now turn to our key results-the NISQ quantum hardware simulation of four-dimensional staggered tesseract HOT lattices.A true 4D lattice is difficult to simulate on most experimental platforms, and with a few exceptions [109,110], most works to date have relied on using synthetic dimensions [27,58,111,112].In comparison, utilizing our exact mapping (Eqs.6 and 7) that exploits the exponentially large many-body Hilbert space accessible by a quantum computer, a tesseract lattice can be directly simulated on a physical 1D spin chain, with the number of spatial dimensions only limited by the number of qubits.The tesseract unit cell can be visual-ized as two interlinked three-dimensional cubes (spanned by x, y, z axes) living in adjacent w-slices (Figure 5).The full tesseract lattice of side length L is then represented as successive cubes with different w coordinates, stacked successively from inside out, with the inner and outer wireframe cubes being w = 1 and w = 6 slices.Being more sophisticated, the 4D HOT lattice features various types of HOT corner, edge, and surface modes (Figure 1a); we presently focus on the fourth-order (hexadecapolar) HOT corner modes, as well as the third-order (octopolar) HOT edge modes.
To start, we realized a dL = 4 × 6 = 24-qubit chain on the quantum processor, which encodes our 6 × 6 × 6 × 6 HOT tesseract.The 4-body (8-operator) interactions now come in d × 2 d = 64 types-half of them are illustrated in Figure 5d, which depicts only the minimal L = 4 case.As previously discussed in Section II A, these interactions are each a product of d − 1 density terms and a hopping process, the latter acting on the particle species that encodes the coupling direction on the HOT tesseract.In generic models with non-axially aligned hopping, these interactions could be a product of up to d hopping processes.As we shortly illustrate, despite the complexity of the interactions, the signal-tonoise ratio in our hardware simulations (Figure 5a) remain reasonably good.
In Figure 5, we consider the configurations C0, C4, C8 and C16, which correspond respectively to the topologically trivial scenario and lattice configurations hosting four, eight, and all sixteen HOT corner modes, as schematically sketched in the header row.Similar to the 2D and 3D HOT lattices, site-resolved occupancy density ρ(x, y, z, w) and occupancy fidelities measured on quantum hardware reveal strong robustness for initial states localized at topological corners, as illustrated by the strongly localized final states in the C4, C8, and C16 cases (Figure 5a).However, their stability is now slightly lower, partly due to the more spacious 4D configuration space into which the state can diffuse, as seen from the colored clouds of partly occupied sites after time evolution.Evidently, the stability diminishes as we proceed to the edge-and surface-localized initial states (Figure 5b).
Next, we investigate a lattice configuration that supports HOT edge modes (or commonly referred to as topological hinge states in literature [48]).So far we have seen topological robustness only from topological corner sites (Figure 5); but with appropriate parameter tuning (see Supplementary Table 1), topological modes can be made to lie along entire edges.This is illustrated in the header row of Figure 6, where topological modes lie along the y-edges.As our HOT lattices are constructed from a mesh of alternating SSH chains, we expect the topological edges to have wavefunction support (nonzero occupancy) only on alternate sites, consistent with the cumulative occupancy densities of the midgap zero-energy modes.This is corroborated by site-resolved occupancy densities and occupancy fidelities measured on quantum hardware, which demonstrate that initial states localized on sites with topological wavefunction support are significantly more robust (Figures 6a-6b), i.e. (x 0 , y 0 , z 0 , w 0 ) = (1, 3, 1, L) overlaps with the topological mode on (1, y, 1, L), y ∈ {1, 3, 5} sites and is hence robust, but (1, 2, 1, L) is not.Stability of the initial state is reduced as we move farther from the corner, as can be seen, for instance, by comparing occupancy fidelities and the size of the final occupancy cloud for (1, 1, 1, L) and (1, 3, 1, L) in Figures 6a-6b, which is expected from the decaying y-profile of the topological edge mode.Finally, our measurements verify that surface-localized states do not enjoy topological protection (Figure 6c) as they are localized far away from the topological edges.It is noteworthy that such measurements into the interior of the 4D lattice can be made without additional difficulty on our 1D qubit chain, but doing so can present significant challenges on other platforms, even electrical (topolectrical) circuits.

F. Resource scaling
Our approach of mapping a d-dimensional HOT lattice onto an interacting 1D chain enabled a drastic reduction in the number of qubits required for simulation, and served a pivotal role in enabling the hardware realizations presented in this work.Here, we further illustrate that employing this mapping for simulation on quantum computers can provide a resource advantage over exact diagonalization3 on classical computers, particularly at large lattice dimensionality d and linear size L.
To be concrete, we consider simulation tasks of the following broad type:given an initial state |ψ 0 ⟩, we wish to perform time-evolution to |ψ(t)⟩, and extract the expectation value of an observable O that is local, that is, ⟨O⟩ is dependent on O l d number of sites on the lattice for a fixed neighborhood of radius l independent of L. State preparation resources for |ψ 0 ⟩ are excluded from our considerations, as there can be significant variations in costs depending on the choice of specification of the state 4 .Measurement costs for computing ⟨O⟩, however, are considered.To ensure a meaningful comparison, we assume first-order Pauli-basis trotterization for the construction of quantum circuits, such that circuit preparation is algorithmically straightforward given a lattice Hamiltonian.As a baseline, exact diagonalization (ED) of a d-dimensional, length L system with a single particle generally requires O(L 3d ) run-time and O(L 2d ) dense classical storage to complete a task of such a type [113].with our approach on a quantum computer (QC) scales with L 2 , rather than the higher power of L 3d through classical exact diagonalization (ED).For (b) fixed L and varying d, our approach scales promisingly, scaling like 4 d instead of (L 3 ) d for ED.We assume conventional trotterization for circuit construction, and at large L and d, our mapping and quantum simulation approach can provide a resource advantage over classical numerical methods (e.g.exact diagonalization).
A direct implementation of a generic Hamiltonian using our mapping gives O dL d • 2 d Pauli strings per trotter step (see Methods), where hoppings along each edge of the lattice, extensive in number, are allowed to be independently tuned.However, physically relevant lattices typically host only a systematic subset of hopping processes, described by a sub-extensive number of parameters.In particular, in the HOT lattices we considered, the hopping amplitude u α r along each axis α is dependent only on α and the parities of coordinates r.Noting the sub-extensive number of distinct hoppings, the lattice Hamiltonian can be written in a more favorable factorized form, yielding O dL • 2 2d Pauli strings per trotter step (see Methods).Decomposing into a hardware gate set, the number of gates in a time-evolution circuit scales as O d 2 L 2 • 2 2d /ϵ in the worst-case for simulation precision ϵ, assuming all-to-all connectivity between qubits.Importantly, imposing nearest-neighbor connectivity on the qubit chain, as in the actual IBM processors (see Figure 8), does not alter this bound.Crucially, there is no exponential scaling of d in L (of form ∼L d ), unlike classical ED.
For large L and d, the circuit preparation and execution time can be lower than the O L 3d run-time of classical ED.We illustrate this in Figure 7, which shows a qualitative comparison of run-time scaling between the quantum simulation approach and ED.We have assumed execution time on hardware to scale as the number of gates on the circuit O d 2 L 2 • 2 2d /ϵ , which neglects speed-ups afforded by parallelization of gates acting on disjoint qubits [114,115].Classical memory usage is similarly bounded during circuit construction, straightforwardly reducible to O(dL) by constructing and executing gates in a streaming fashion [116,117], and worstcase O 2 ld during readout to compute ⟨O⟩, reducible to a constant supposing basis changes to map components of O onto the computational basis of a fixed number of measured qubits can be implemented on the quantum circuits [118,119].
The favorable resource scaling (run-time and memory), in combination with the modest dL qubits required, suggests promising scalability of our mapped quantum simulation approach, especially in realizing larger and higher-dimensional HOT lattices.We re-iterate, however, that trotterized circuits without additional optimization remain largely too deep for present-generation NISQ hardware to execute feasibly.The use of qudit hardware architectures in place of qubits can allow shallower circuits [120][121][122]; in particular, using a qudit of local Hilbert space dimension ≥ 2 d instead of a group of d qubits avoids, to a degree, decomposition of long-range multi-site gates, assuming the ability to efficiently and accurately perform single-and two-qudit operations [123,124].Nonetheless, for the quantum simulation of sophisticated topological lattices as described to be achieved in its full potential, fault-tolerant quantum computation, at the least quantum devices with vastly improved error characteristics and decoherence times, will likely be needed.

G. Discussion and outlook
Although there exists a plethora of metamaterials and classical systems demonstrating high-dimensional topological or HOT phases [48,57,92,97,100,101,103,[125][126][127][128][129], these paradigmatic topological states have not been directly simulated on NISQ quantum platforms.Recent demonstrations of various enigmatic condensedmatter phenomena on digital quantum computers, ranging from discrete time crystals [130,131] to topological phases [81,[132][133][134][135][136], illustrate promising capabilities despite hardware constraints, such as limited gate fidelities and decoherence times.Despite the rapid pace of quantum hardware development and algorithmic advancements, the simulation of HOT phases on quantum computers remain inherently difficult, given their prohibitive qubit number requirements.
By fully exploiting the exponentially large many-body Hilbert space accessible on a quantum computer, we achieved the first simulation of HOT Hamiltonians in up to four dimensions in a quantum setting.We not only accurately measured the density evolution of initial states and demonstrated their robustness arising from higher-order topology, but also conclusively measured their midgap topological energy spectra, which report on the number of degenerate HOT modes hosted by the high-dimensional lattice.In overcoming the challenges associated with NISQ quantum simulation, we pave the way for further experimental investigation of HOT phenomena, valuable especially given the paucity of quantum material HOT candidates.
Though outside of the present HOT scope, we remark that a similar mapping scheme apply also to interacting lattices, at an expense of a larger local Hilbert space dimension for each quantum chain site and more complicated couplings.Namely, to simulate an interacting d-dimensional lattice containing up to m particles, one can utilize m copies of the d species of particles on the chain, and in the worst-case, hopping processes on the lattice translate into interactions involving the md particles on the chain.In such a scheme, O(mdL) qubits suffice to implement the chain, with possible further reductions depending on the lattice.
III. ACKNOWLEDGEMENTS J.M.K. and T.T. thank Wei En Ng and Yong Han Phee of the National University of Singapore for helpful discussions.The authors acknowledge the use of IBM Quantum services for this work.The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.

IV. AUTHOR CONTRIBUTIONS
C.H.L. initiated and supervised the project, and wrote parts of the manuscript.J.M.K. developed most of the quantum simulation codebase, ran experiments on emulators and hardware, and wrote parts of the manuscript.T.T. contributed to the codebase, ran experiments, and wrote parts of the manuscript.C.H.L. acknowledges the National Research Foundation Singapore grant under its QEP2.0 programme (NRF2021-QEP2-02-P09).All authors analyzed computational and experiment results.The manuscript reflects the contributions of all authors.

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

VI. CODE AVAILABILITY
The code used in this study is available from the corresponding author upon reasonable request.

METHODS
Quantum hardware.We utilized IBM transmon-based superconducting quantum devices in our experiments.For time-evolution on the square lattice, we used QV-32 devices ibmq_manhattan (65 qubits) and ibmq_brooklyn (65 qubits); for time-evolution on the cubic and tesseract lattices, we used QV-32 devices ibmq_sydney (27 qubits) and ibmq_toronto (27 qubits), and QV-128 devices ibmq_mumbai (27 qubits) and ibmq_montreal (27 qubits).The latter group of QV-128 devices were also used for iterative quantum phase estimation (IQPE).The quantum volume (QV) reflects an approximate measure of the aggregate capability of the machine-number of qubits, gate erorr rates, and decoherence times [137,138].To provide ballpark measures of performance, the relaxation T1 and dephasing T2 times range 80 µs ≤ T1 ≈ T2 ≤ 130 µs on average on the devices.Typical single-qubit ( √ X) gate errors are O 10 −4 , and typical two-qubit (CX) gate errors are O 10 −2 .Illustration of qubit layouts for 27-qubit and 65-qubit devices are provided in Figure 8.
Hardcore bosons.Our prescribed mapping transforms ddimensional non-interacting lattices to one-dimensional interacting chains hosting hardcore bosons {ω α } d α=1 , which satisfy the standard hardcore bosonic mixed commutation relations for sites ℓ, ℓ ′ on the chain.The anti-commutation relations forbid double-occupancy of any site ℓ on the chain by two bosons of the same species.We remark that, in the singleparticle context of HOT lattices, particle statistics can be neglected since there is no second particle on the d-dimensional lattice to exchange with.The choice of hardcore bosons over ordinary bosons in our mapping can be viewed as a form of Fock space truncation supporting simulation on quantum computers; occupancy-dependent phase factors are also avoided compared to a fermionic mapping.We may likewise choose {ω α } d α=1 to be mutually commuting.Restricted sector.The symmetries of our Hamiltonians impose constraints on their physically allowed states.The number-conserving properties of H dD lattice and H dD chain in each species, and the nature of our mapping in encoding coordinates along axes α as the location of a species-α boson along the chain, implies that exactly one boson of each species occupies the chain at any time.States lying outside this Fockspace sector, which we call the restricted sector, are unphysical.As described later, we take advantage of this constraint in our recompilation and post-selection procedures.
Representation in terms of qubits.By definition, where |m 1 ℓ , . . ., m d ℓ ⟩ are Fock basis states at site ℓ on the chain, with m β ℓ ∈ {0, 1} indicating the occupancy of species β.Using d qubits to represent each site of the chain, we associate where |0⟩ q , |1⟩ q are the computational basis states of a qubit.This fixes the representation for Pauli raising and lowering operators σ ± = (σ x ± iσ y )/2, and σ γ ℓα denotes σ γ acting on the α th qubit of the ℓ th group, representing the ℓ th site, each group comprising d qubits.This representation satisfies the mixed hardcore bosonic commutation relations of ω α and (ω α ) † .The number operator for projector I z ℓα = (I ℓα − σ z ℓα )/2 onto the |1⟩ q subspace of the α th qubit of the ℓ th group.Terminal computational-basis (σ z ) measurements on each qubit on the chain thus suffices to produce site-resolved occupancy density estimates throughout the lattice.Specifically, ρ(r is the probability of detecting |1⟩ q outcomes on all qubits representing site r of the lattice.We adopt the convention that a measurement records bit m ∈ {0, 1} for a |m⟩ q outcome, thus the outcomes for dL qubits are dL-length bitstrings.Then, ρ(r) is the fraction of bitstrings with 1s on all qubits {rαα} d α=1 .The statistical reliability of ρ(r) improves with the number of measurement samples; in our experiments, we perform ≥ 32000 shots for each simulation circuit.
Trotterization.Given a generic Hamiltonian H, conventional trotterization decomposes H in the spin-1/2 basis and performs time-evolution through a series of small time steps.Writing H = γ Aγ σ γ where Aγ ∈ C and σ γ are Pauli strings, the first-order trotterization scheme is for M time steps (trotter steps) each covering ∆t = t/M .The error term arises from the non-commutation of σ γ strings in H, and implies that a number of steps scaling as 1/ϵ are required to achieve a simulation precision ϵ.Each of the e −iAγ σ γ ∆t terms can be implemented via standard quantum circuit primitives, thus allowing a straightforward transcription of the time-evolution circuit.There exist higher-order formulae of error O 1/M 2 that produce deeper circuits.In the context of HOT lattices, substituting our representation of hardcore bosonic operators (Eqs.12 and 13) into H dD chain (Eq.7) gives the spin-1/2 decomposition Imposing nearest-neighbor (NN) qubit connectivity, however, necessitates SWAP gates to effect CX operations between non-adjacent qubits; then the maximal distance between qubits acted upon by the Pauli string determines circuit depth.Here, the maximal is bounded by the number of qubits dL on the chain, thus each exponentiated Pauli string can be implemented with O(dL) gates, using a cascading SWAP scheme.Thus, a time-evolution circuit depth of O(d 2 L d • 2 d /ϵ) and O(d 2 L d+1 • 2 d /ϵ) are expected for allto-all and NN qubit connectivities respectively, at simulation precision ϵ.
In the above decomposition (Eq.15) and analysis for generic Hamiltonians, we have not assumed any particular structure to the hoppings u α r ; thus the results are applicable to any d-dimensional lattice hosting NN hoppings.In our HOT lattices, however, u α r depends only on α and the parities of coordinates r, alternating between intra-and inter-cell counterparts.This structure allows a more compact decomposition within the physically relevant sector, where u α π is now labeled by the parity vector π of lattice coordinates, and the parity operator Recompilation.Circuit recompilation starts from a circuit ansatz, whose parameters are dynamically optimized [81][82][83][84][85]139].We use a circuit ansatz comprising an initial layer of U3 general single-qubit rotation gates on all qubits followed by K ansatz layers, each comprising a layer of CX gates entangling adjacent qubits and a layer of U3 rotations (see Figure 2e).This ansatz provides sufficient entangling power to accommodate the couplings in H dD chain we investigated, at modest K ≤ 10 depth.Each U3 gate in the ansatz is associated with rotation angles (θ, ϕ, λ); we collate all of them into parameter vector ϑ.Then, given a target circuit unitary V and an initial state |ψ0⟩, the optimization problem can be numerically treated, where V ϑ is the circuit ansatz unitary with parameters ϑ.The recompiled circuit is then the ansatz with optimal parameters ϑ * fixing the U3 gates.
To enhance recompilation performance and the quality of recompiled circuits, we note that we require access only to the restricted sector of H dD chain , as established above; thus we may focus optimization within this sector.Concretely, we treat where fidelity F R is taken over the restricted sector, such that loss evaluation is made vastly cheaper.Note that postselection enforces the symmetries of H dD chain that forbid access to states outside the restricted sector.A complication is that F R no longer has fixed normalization; to stabilize the optimization, we impose a minimum normalization baseline Λ0, by regularizing the sector-specific fidelity with a softmax Ω(x, Λ) = Λ + log [1 + e k(x−Λ) ]/δ of reasonable sharpness δ.In our implementation, we estimate V θ through autodifferentiable tensor network-based ansatz simulation [140], and use L-BFGS-B with basin-hopping to perform the optimization [141,142], for ≤ 10 2 hops and ≤ 10 3 iterations per hop, terminated at estimated F R ≥ 99.9%.Recompilation is tuned to never exceed a few minutes per circuit.
Iterative quantum phase estimation (IQPE).Given U |ψ⟩ = e 2πiϕ |ψ⟩ for a unitary U and an eigenstate |ψ⟩, quantum phase estimation (QPE) measures eigenphase ϕ ∈ [0, 1), in principle to arbitrary precision [90,91,143,144].Setting U = e −iHt allows the inference of eigenenergy E = −2πϕ/t of |ψ⟩, enabling the probing of the spectrum of H.In standard circuit constructions of QPE, each bit of ϕ in binary is measured by an ancillary qubit.In comparison, IQPE uses a single ancillary qubit and does not require multi-qubit quantum Fourier transforms [88,89], and is thus preferable on NISQ devices.
Truncating ϕ = 0.ϕ1ϕ2 . . .ϕm to m bits, the IQPE algorithm iterates k = m to k = 1.In the circuit for the k = m iteration, a controlled-U 2 m−1 block is applied, and the ancilla qubit is measured to determine ϕm.The probability p0 of measuring an ancilla state of |0⟩ q is cos 2 [(0.ϕm)π] = 1 − ϕm in the absence of noise [88]; amidst noise, or when |ψ⟩ is not an exact eigenstate, the ancilla measurement is no longer deterministic, but a thresholded inference of ϕm = 0 if p0 > 1/2 and ϕm = 1 otherwise can still be applied.Subsequently, in iteration k, a controlled-U 2 k−1 block is applied, and a feedback rotation of angle ω k = −2π(0.0ϕk+1 ϕ k+2 . . .ϕm) is applied on the ancilla to remove the phase due to the previous bits, before likewise inferring ϕ k .The circuit schematic for iteration k is illustrated in Figure 2c.The energy resolution of IQPE is determined by the number of iterations m executed.
To probe corner HOT modes, we select |ψ⟩ to be perfectly localized on lattice corners.More precisely, we pick |ψ⟩ to be simple superpositions of corner-localized states, to emulate the generic profiles of topological eigenstates on finite-size lattices.Despite their simplicity, such |ψ⟩ closely resemble exact HOT modes (≥ 80% fidelities).Readout error mitigation is applied to all qubits, and post-selection is applied to the simulation qubits to select for the physically relevant restricted sector.
Readout error mitigation (RO).Readout error refers to the chance of measuring |1⟩ q when a qubit is in |0⟩ q or vice versa, which is non-negligible on NISQ devices.A standard method to mitigate these errors is to characterize the measurement bit-flip probabilities of each qubit; then linear inversion can be performed on raw measurement counts in experiments to recover approximate true measurement counts [3,19,86,145].We first describe complete readout error mitigation.For terminal measurements on N qubits, where N = dL for time-evolution and N = dL + 1 for IQPE, the measurement bitstrings that can result are S = {0, 1} N .Given a calibration matrix M with entry M ss ′ recording the probability of measuring s ∈ S when the true result should be s ′ ∈ S, and denoting the raw measurement count cs and corrected count c ′ s for bitstring s ∈ S, we may write where c = c0 c1 . . .c 2 N −1 T and M + is the pseudoinverse of M .The estimated c ′ may carry negative entries due to the approximate M ; we zero these entries as they are unphysical.A least-squares fit with non-negative constraints can alternatively be used [3,146], but we chose the pseudoinverse to reduce post-processing costs.The matrix M is constructed at the start of each run by running calibration circuits, which prepare the N qubits in state |s⟩ for each s ∈ S and collect measurement probabilities.Accordingly, 2 N calibration circuits are required.This is prohibitively expensive for large N , which motivates tensored readout error mitigation.
In the tensored scheme, we partition the N qubits into k sub-registers, containing N1, . . ., N k qubits.Then calibration matrices M 1 , . . ., M k can be acquired separately for each sub-register, and c ′ = (M 1 ⊗ . . .⊗ M k ) + c likewise applies.In practice, for efficiency, we avoid the explicit tensor product by block-wise operating on c with the pseudoinverses of the calibration matrices.The scheme requires 2 max(N 1 ,...,N k ) < 2 N calibration circuits, since circuits operating on disjoint subsets of qubits can be merged.This reduction in cost, however, is at the expense of neglecting correlations in readout errors between qubits in different sub-registers.In our experiments, we use up to k = 4 sub-registers with up to 8 qubits in each.
Post-selection (PS).We exploit the symmetry constraints of H dD chain to mitigate gate-round errors incurred in our simulation circuits, to a limited extent.Specifically, recall that only the restricted sector of H dD chain is physically allowed.Occupancies outside of the sector are unphysical and should not arise during simulation-measurements that detect so thus indicate erroneous results.At no additional cost, we may piggyback on the terminal computational-basis measurements on simulation qubits, which provide site-resolved occupancy densities, to filter and discard such noise-afflicted outcomes.Explicitly, we accept counts c ′ s of measurement bitstring s ∈ S ′ , where s ℓα denotes the measurement outcome on the α th qubit of the ℓ th group in bitstring s.We zero the counts c ′ s for s / ∈ S ′ , which violate symmetries and are unphysical.
Qubit selection and averaging.The quantum devices we utilize provide more qubits than needed for our experiments, and there are significant variations in the quality of qubits within the same device.We hence use an exhaustive search scheme to select qubits of the lowest estimated error rates to use in our simulations.Given the available CX couplings between qubits, we identify all distinct qubit chains of length N , and compute the following fitness function for each chain, where ε M i is the calibrated readout error for qubit i, and ε CX i,j is the calibrated CX gate error between qubits i and j.By convention, ε CX 1,0 = ε CX N,N +1 = 0.The fitness function above is designed to emulate the structure of our recompilation circuit ansatz, with K ∼ 10 ansatz layers.We pick the qubit chains with highest Q for our simulations.This selection is only approximate in nature, since the error rates are obtained from the scheduled calibration (∼ daily) of quantum devices and are subject to drift and fluctuations over time.Illustrations of qubit chains selected through this procedure are shown in Figure 8 on 65-and 27-qubit processors.Moreover, to average out fluctuations in data and further reduce the effects of hardware noise, we collate error-mitigated measurement bitstring counts over multiple qubit chains per device as selected by our search, and multiple devices (as listed above).
Parameter selection.The HOT lattice models used can be built from an alternating d-dimensional patchwork of alternating Su-Schrieffer Heeger (SSH) chains [80,102].Each SSH-type chain is characterized by a topological invariant, its winding number.That is, each constituent SSH chain in our lattice is topological when its corresponding hopping strengths satisfy |v| < |w|.
All of our HOT corner states are categorized as type-I [102], with support only on a single sublattice and exponentially decaying amplitudes in all orthogonal directions.This is interpreted as a higher-dimensional manifestation of the 1D topology of the SSH model, such that the winding numbers of the edges that intersect a topological corner are all 1.We provide a complete list of hopping parameters v α π , w α π for our models in Supplementary Table 1.For a d-dimensional lattice with E number of edges, appropriate choices of parameters can be systematically found using simple offsets {δj} E j=1 such that the SSH-like coupling parameters satisfy v = λ(1 − δ), w = λ(1 + δ), with λ = −1 for a subset of edges chosen to generate a π flux through each face of the lattice, thereby gapping the spectrum [47].A positive (negative) choice of δ for an edge results in it being topologically non-trivial (trivial).
On the 2D lattice, to achieve midgap 1D edge modes, one can tune the hopping parameters away from a configuration with four corner modes (C4) towards one with only two neighboring corner modes.In this transition, two of the corner modes vanish as the bands become gapless along one momentum direction, leading to the emergence of zero-energy 1D edge modes.This approach straightforwardly generalizes to engineer zero-energy (higher-order) edge modes in higherdimensional lattices, for example on the tesseract lattice.Supplementary Table 1.Values of intra-and inter-cell hopping coefficients, respectively v α π(rα) and w α π(rα) along direction α ∈ {x, y, z, w} and at parity π(rα) on the remaining axes, used in simulations of the square, cubic, and tesseract HOT lattices.Numerical values are explicitly listed for clarity; a systematic procedure to assign coefficients for corner-mode configurations on square and cubic lattices is available in Ref. [1], generalizable to higher dimensions.

FIG. 1 .
FIG.1.Higher-order topological (HOT) states in higher-dimensional lattices.(a) In the different number of dimensions d, n th order topological states manifest as robust codimension-n corner, edge, surface, and volume states protected by higher-order topological invariants.HOT states refer to topological states with n > 1.(b) Schematics of HOT corner modes on hyperlattices in different dimensions.The illustrated lattice sizes are small for visual clarity; the HOT states we realized on quantum hardware are on considerably larger 16 × 16, 6 × 6 × 6, and 6 × 6 × 6 × 6 lattices.

FIG. 2 .
FIG.2.High-level schematic of our approach on simulating high-dimensional lattice models on quantum hardware.(a)-(b) Mapping of a higher-dimensional lattice to a one-dimensional interacting chain to facilitate quantum simulation on near-term devices.Concretely, a two-dimensional single-particle lattice can be represented by a two-species interacting chain; a three-dimensional lattice can be represented by a three-species chain with three-body interactions.(c) Overview of quantum simulation methodology: higher-dimensional lattices are first mapped onto interacting chains, then onto qubits; various techniques, such as (d) trotterization and (e) ansatz-based recompilation, enable the construction of quantum circuits for dynamical time-evolution, or iterative quantum phase estimation (IQPE) for probing the spectrum.The quantum circuits are executed on the quantum processor, and results are post-processed with readout error mitigation (RO) and post-selection (PS) to reduce effects of hardware noise.See Methods for elaborations on the mapping procedure, and quantum circuit construction and optimization.

FIG. 3 .
FIG. 3. Quantum processor measurements of 2D HOT zero modes and their roles in preserving state fidelity.(a) Ordered eigenenergies on a 10 × 10 lattice for the topologically trivial C0 and nontrivial C2 and C4 configurations.They correspond to 0, 2 and 4 midgap zero modes (red diamonds), as measured via IQPE on a 20-qubit quantum chain plus an additional ancillary qubit; the shaded red band indicates the IQPE energy resolution.The corner state profiles (right insets) and other eigenenergies (black and gray dots) are numerically obtained via ED.Time-evolution of four initial states on a 16 × 16 lattice mapped onto a 32-qubit chain-(b)-(c) localized at corners to highlight topological distinction, (d) localized along an edge, and (e) delocalized in the vicinity of a corner.Left plots show occupancy fidelity for the various lattice configurations, obtained from exact diagonalization (ED) and hardware (HW), with insets showing the site-resolved occupancy density ρ(x, y) of the initial states (darker shading represents higher density).Right grid shows occupancy density measured on hardware at two later times.States with good overlap with robust corners exhibit minimal evolution.Error bars represent standard deviation across repetitions on different qubit chains and devices.In general, the heavy overlap between an initial state and a HOT eigenstate confers topological robustness, resulting in significantly slowed decay.(f ) Schematic of the interacting chain Hamiltonian, mapped from the parent 2D lattice, illustrated for a smaller 6 × 6 square lattice.The physical sites of the interacting boson chain are colored black, with their many-body interactions represented by colored vertices.

FIG. 4 .
FIG. 4. Quantum processor measurements of HOT corner modes on the staggered cubic lattice.(a) Header row displays energy spectra for the topologically trivial C0 and inequivalent nontrivial C4a, C4b and C8 configurations.The configurations host 0, 4 and 8 midgap zero modes (red diamonds), as measured via IQPE on an 18-qubit chain plus an ancillary qubit; the shaded red band indicates the IQPE energy resolution.Schematics illustrating the locations of topologically robust corners are shown on the right.Subsequent rows depict time-evolution of five initial states on a 6 × 6 × 6 lattice mapped onto an 18-qubit chain-localized at a corner, on an edge, on a face, and in the bulk of the cube, and delocalized in the vicinity of a corner.Leftmost column plots occupancy fidelity for the various lattice configurations, obtained from exact diagonalization (ED) and hardware (HW), with insets showing the site-resolved occupancy density ρ(x, y, z) of the initial state (darker shading represents higher density).Central grid shows occupancy density measured on hardware at a later time (t = 0.6), for the corresponding initial state (row) and lattice configuration (column).Error bars represent standard deviation across repetitions on different qubit chains and devices.Again, initial states localized close to topological corners exhibit higher occupational fidelity.(b) Hamiltonian schematic of the interacting chain realizing a minimal 4 × 4 × 4 cubic lattice.Sites on the chain are colored black; colored vertices connecting to multiple sites on the chain denote interaction terms.

TimeFIG. 5 .
FIG. 5.Quantum processor measurements of HOT corner modes on a 4D tesseract lattice.A L = 6 tesseract lattice is illustrated as six cube-slices indexed by w and highlighted on a color map.Header row displays energy spectra computed numerically for the topologically trivial C0 and nontrivial C4, C8, and C16 configurations.The configurations host 0, 4, 8, and 16 midgap zero modes (black circles).Schematics on the right illustrate the locations of the topologically robust corners.Subsequent rows depict time-evolution of three initial states on a 6 × 6 × 6 × 6 lattice mapped onto a 24-qubit chain-localized on (a) a corner, (b) an edge, and (c) a face.Leftmost column plots occupancy fidelity for the various lattice configurations, obtained from exact diagonalization (ED) and quantum hardware (HW), with insets showing the site-resolved occupancy density ρ(x, y, z, w) of the initial state.Central grid shows occupancy density measured on hardware at the final simulation time (t = 0.6), for the corresponding initial state (row) and lattice configuration (column).Color of individual sites (spheres) denote their w-coordinate and color saturation denote occupancy of the site; unoccupied sites are translucent.Error bars represent standard deviation across repetitions on different qubit chains and devices.Initial states with less overlap with topological corners exhibit slightly lower stability than their lower dimensional counterparts, as these states diffuse into the more spacious 4D configuration space.(d) Hamiltonian schematic of the interacting chain realizing a minimal 4 × 4 × 4 × 4 tesseract lattice.Sites on the chain are colored black; colored vertices connecting to multiple sites on the chain denote interaction terms.To limit visual clutter, only v α π couplings are shown; a corresponding set of w α π couplings are present in the Hamiltonian but have been omitted from the diagram.

(FIG. 6 .
FIG. 6.Quantum processor measurements of robustness from HOT edge (or hinge) modes on a 4D tesseract lattice.Our mapping facilitates the realization of any desired HOT modes, beyond the aforementioned corner mode examples.Header row on the left displays energy spectrum for a configuration of the tesseract harboring topologically non-trivial edges (midgap mode energies in black).Accompanying schematic highlights alternating sites with topological edge wavefunction support.Subsequent columns present site-resolved occupancy density ρ(x, y, z, w) for a 6 × 6 × 6 × 6 lattice mapped onto a 24-qubit chain, measured on quantum hardware at t = 0 (first row) and final simulation time t = 0.6 (second row), for three different experiments.(a) A corner-localized state along a topological edge is robust, compared to one along a non-topological edge.(b) On a topologically non-trivial edge, a state localized on a site with topological wavefunction support is robust, compared to one localized on a site without support.(c) A surface-localized state far away from the topological edges diffuses into a large occupancy cloud.Bottom leftmost summarizes occupancy fidelities for the various initial states, obtained from exact diagonalization (ED) and hardware (HW).Error bars represent standard deviation across repetitions on different qubit chains and devices.

FIG. 7 .
FIG.7.Favorable resource scaling.Comparison of asymptotic computational time required for the dynamical simulation of d-dimensional, size-L lattice Hamiltonians of similar complexity as our HOT lattices.With (a) fixed lattice dimension d and increasing lattice size L, the time taken with our approach on a quantum computer (QC) scales with L 2 , rather than the higher power of L 3d through classical exact diagonalization (ED).For (b) fixed L and varying d, our approach scales promisingly, scaling like 4 d instead of (L 3 ) d for ED.We assume conventional trotterization for circuit construction, and at large L and d, our mapping and quantum simulation approach can provide a resource advantage over classical numerical methods (e.g.exact diagonalization).

17 FIG. 8 .
FIG. 8. Qubit layout on quantum hardware.Diagram showing qubit layout and connectivity on (a) 65-qubit and (b) 27-qubit quantum processors.Native CX gates are available on pairs of qubits joined by an edge.An example of a 32-qubit chain in (a) and an 18-qubit chain in (b), used for 16 × 16 square and 6 × 6 × 6 cubic HOT lattice simulations respectively, are highlighted in gray.The selection of qubit chains is performed via a search scheme that minimizes an estimated error metric (see Methods).
so that naively counting, there are O(dL d ) coupling terms each hosting O(2 d ) Pauli strings, totalling O(dL d • 2 d ) Pauli strings comprising H dD chain and thus manifesting in each trotter step.These Pauli strings have ≤ d weight, that is, they act non-trivially on ≤ d qubits, so their exponentiated forms can each be implemented with O(d) gates, assuming all-to-all qubit connectivity on hardware.
P β τ |ψ⟩ = |ψ⟩ for |ψ⟩ harboring a particle on a parity-τ position along axis β and P β τ |ψ⟩ = 0 otherwise (parity 1−τ ).Instead of having O(dL d ) unique coupling terms, there are O(d • 2 d ) coupling terms each hosting O(L • 2 d ) Pauli strings, thus totaling O(dL • 2 2d ) Pauli strings comprising H dD chain .The strings are of weight O(dL), so for both all-to-all and NN qubit connectivity, a time-evolution circuit depth of O(d 2 L 2 • 2 2d /ϵ) is expected.In comparison to the O d 2 L d+1 • 2 d /ϵ depth for the naive trotterization above with NN qubit connectivity, the present scheme turns the factor L d → L•2 d , which is significantly advantageous for L ≫ 1.