Quantum dynamical characterization and simulation of topological phases with high-order band inversion surfaces

Xiang-Long Yu1,∗ Wentao Ji2,3,∗ Lin Zhang4,5,∗ Ya Wang2,3,† Jiansheng Wu1,6,‡ and Xiong-Jun Liu4,5,1§ 1 Department of Physics and Shenzhen Institute for Quantum Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, P. R. China 2 Hefei National Laboratory for Physical Sciences at the Microscale and Department of Modern Physics, University of Science and Technology of China, Hefei 230026, P. R. China 3 CAS Key Laboratory of Microscale Magnetic Resonance, University of Science and Technology of China, Hefei 230026, P. R. China 4 International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, P. R. China 5 Collaborative Innovation Center of Quantum Matter, Beijing 100871, P. R. China and 6 Guangdong Provincial Key Laboratory of Quantum Science and Engineering, Shenzhen Institute for Quantum Science and Engineering, Southern University of Science and Technology, Shenzhen 518055, Guangdong, P. R. China


I. INTRODUCTION
The bulk-boundary correspondence (BBC) is a fundamental mechanism in topological quantum phases, such as in quantum Hall effect, topological insulators, and topological superconductors, in which the topological number of the bulk links to the number of the robust gapless states on the boundary [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. In consequence, the BBC is widely utilized to identify the topological quantum states and measure the topological invariants, e.g. by transport and angle-resolved photoelectron spectroscopy (ARPES) experiments [18][19][20][21][22][23][24]. Recently, a new development has been made for topological systems, called the higher-order topological phases [25][26][27][28][29][30][31]. In such quantum systems, a d-dimensional (dD) bulk has nontrivial topology in the presence of a certain spatial symmetry, while the corresponding gapless modes do not exist on the (d − 1)D boundary but survive on a lower (d − n)D boundary with n > 1, due to spatial symmetry breaking on the former but preserving on the later. This gives a high-order BBC in the real space.
While the BBC is well defined in the real space, it seems surement of the topological phases are also confirmed in experiments based on dBSC compared with traditional methods built on equilibrium theories [45][46][47][48][49][50][51].
In this work, we propose the concept of high-order BISs for the dynamical characterization and simulation of topological quantum phases by extending the dBSC to high order, and report the experimental observation. The nth-order BIS is a (d−n)D interface on the (n−1)th-order BIS. From a dimension reduction approach, we show that through quantum quenches the bulk topology of a dD phase can be uniquely characterized by the dynamical topology emerging on any high-order BISs, manifesting the high-order dBSC. The prediction provides the optimal schemes and considerably expand the freedom for dynamical characterization and simulation of topological phases. Experimentally, we build up a quantum simulator using nitrogen-vacancy (NV) center to investigate the high-order dBSC in a 3D chiral topological phase through emulating each momentum one by one, with the great advantages of the minimal measurement strategy in the simulation being demonstrated.
The paper is organized as follows. In Sec. II, we introduce the concept of high-order BISs and topological invariants from dimension reduction. In Secs. III and IV, we propose two dynamical schemes for the characterization and simulation of topological phases. In Sec. V, we experimentally build up a quantum simulator using NV center for the dynamical characterization of topological phases. Finally, we conclude in Sev. VI with an outlook. Technical details are provided in Appendices.

II. HIGH-ORDER BAND INVERSION SURFACES AND TOPOLOGICAL CHARACTERIZATION
We start with the dD topological phases classified by integer invariants in the Altland-Zirnbauer (AZ) symmetry classes, which are winding numbers in odd dimensions and Chern numbers in even dimensions [7,53,54]. The Hamiltonian is described in the form [55,56], where the Clifford algebra matrices γ i satisfies anticommunication relation, γ i , γ j = 2δ i j , and has matrix dimension n d = 2 d/2 (or 2 (d+1)/2 ) when the spatial dimension d is even (or odd). For 1D and 2D systems, the γ matrices are the Pauli ones and Eq. (1) describes a two-band model, such as the Su-Schrieffer-Heeger chain [57] and the Haldane model of the integer quantum Hall effect [58]. For 3D and 4D systems, the γ matrices are the Dirac ones and the corresponding Hamiltonian describes a four-band model, such as 3D DIII class superconductors and AIII class topological insulators [53], and 4D quantum Hall insulators [59]. In general, the γ matrices can be regarded as the spin or pseudospin operators. The bulk topology defined in the whole Brillouin zone (BZ) is quantified by a dD invariant, and can be reduced to the lower-dimensional topology on the (d − 1)D BIS which is obtained by choosing one component, e.g. the i 0 th component of h(k) to satisfy h i 0 (k) = 0, with i 0 ∈ {0 ∼ d}, and is a (d − 1)D closed surface in the BZ [32]. Here we name it as the first-order BIS, 1-BIS = {k ∈ BZ|h i 0 (k) = 0}, which can be interpreted in a physically transparent way as follow. We take that the h i 0 -term characterizes the energy dispersions of the bulk bands, while the remaining terms of the Hamiltonian, denoted by h (1) = {h i 1 , h i 2 , · · · , h i d }, characterize the coupling (pseudospin-orbit coupling) between different bands. Then the 1-BIS is simply a momentum subspace, across which the energy of the half of these bulk bands is inverted with respect to that of the other half before taking into account the interband couplings. The interband coupling term h (1) opens a gap on the BIS, leading to the gapped topological phase, and moreover, the vector field h (1) exhibits nontrivial topology on the 1-BIS, which is captured by a (d − 1)D topological invariant (winding or Chern number) on the 1-BIS. This renders the bulk-surface duality of the free-fermion topological phases with integer invariants [32].
The central concept we propose here is the high-order BIS which can be introduced based on a dimension reduction method. A key observation is that the interband coupling term h (1) can be further treated as a (d − 1)D gapped Hamiltonian with integer invariant defined on the 1-BIS. Correspondingly, we can define a higher-order (second-order) BIS on the 1-BIS via another component h i 1 of h (1) , namely 2-BIS = {k ∈ 1-BIS|h i 1 (k) = 0}, which equals to {k ∈ BZ|h i 0 (k) = h i 1 (k) = 0} and is a (d −2)D closed surface. In a similar process we can show that the bulk topology is reduced to the (d − 2)D invariant of the (d − 1)D vector field h (2) = {h i 2 , h i 3 , · · · , h i d } on the 2-BIS (see details in Appendix A). Repeating this dimension reduction we obtain the nth-order BIS by where h i α 's are n components of h. Obviously, the n-BIS is a (d − n)D closed surface and formed by the k points satisfying α=1∼(n−1) h i α (k) = 0. Across the n-BIS the band energies switch sign in the absence of h (n) = {h i n , h i n+1 , · · · , h i d }. With the above process the bulk topology then reduces to the (d − n)D invariant obtained by h (n) on the n-BIS that where ] m is the topological invariant of the vector field f on the mD manifold M (see details of proof in Appendix A) and characterizes the integer times thatf covers over the corresponding spherical surface S m when k runs over M. Eq. (3) shows a correspondence that maps the classification of bulk topology to the characterization on the n-BIS. This correspondence is valid for any highorder BISs. In addition, two points are worthwhile to mention for the dimension reduction process. Firstly, the configurations of the high-order BISs are sharply different if choosing different components (h i ) of the Hamiltonian for the definition in Eq. (2), with two typical cases being illustrated in Fig. 1 (a) and (b), respectively, where we particularly sketch the dimension reduction from 3D to 0D in the last three steps. Secondly, whenever the BIS of any certain order does not exist in the BZ, the corresponding system is topologically trivial.
The highest-order d-BIS is of zero dimension and is obtained for n = d, and consists of a finite number of momentum points. Thus the corresponding topological invariant T d-BIS [h i d ] only depends on a finite number of momentum points where j corresponds to different sectors of the d-BIS. The subscripts L j and R j are the left-and right-hand points of the jth sector, respectively. They appear in pairs due to the fermion doubling theorem [60]. The Left and the Right depend on the direction of the integral path on the 1D (d−1)-BIS. This formula is similar to the Brouwer degree of the map h(k) from k space to h space [61] and it maps the classification of bulk topology to the characterization in zero dimension. In this case, the computation of the topological invariant can be greatly optimized. In the following we further apply the dimension reduction based on high-order BISs to the dynamical characterization of topological phases.

III. DYNAMICAL CHARACTERIZATION AND SIMULATION: SCHEME I
We turn to dynamical characterization of topological phases by showing the high-order dBSC and present quantum simulation with concrete models. The quantum dynamics are induced by quenching an initially fully polarized trivial phase to topological regime. The pre-and post-quench Hamiltonians are H pre = H k + δm j γ j with δm j H k and H post = H k , respectively, where j denotes the quench axis. Hence, the initial phase is the trivial ground state of δm j γ j . The spin precesses with respect to the vector field h of H post after quenching. The characterization is captured by three main steps: (i) Choose a certain direction to perform a deep quench and measure the spin polarization along the same direction. Then the 1-BIS is formed by the momenta with vanishing time-averaged spin polarization (TASP), given by where i and j correspond to the spin polarization and quench axes, respectively, and ρ j (0) is the density matrix of the initial state at the time t = 0. The Eq. (5) can be simplified as γ i j = −h i h j / d i =0 h 2 i , and the 1-BIS is determined by γ j j = 0 which corresponds to h j = 0. On the 1-BIS the vector field h (1) is perpendicular to the initial spin polarization, leading to a spin precession within the plane perpendicular to h (1) . (ii) Repeating the quench process with respect to a new axis γ j , one measures another 1-BIS defined by γ j j = 0. The intersection between the two 1-BISs gives the 2-BIS. With this process we can perform sequential dimension reduction of the BISs, and determine the n-BIS by n-BIS = {k| γ 0 0 = γ 1 1 = · · · = γ n−1 n−1 = 0}, where without losing generality we have set the quench directions i = 0 to n − 1 in turn, and in each quench only the ith spin component needs to be measured. Note that to determine the j-BIS, the spin dynamics are need to be measured only on the subspace of ( j − 1)-BIS, so the number of momentum points need to be measured decreases rapidly. (iii) One can verify that at the momentum k away from BISs the TASP is finite. Measure spin polarizations along the remaining d − n + 1 directions near the n-BIS to obtain a dynamical spin texture (DST) on the n-BIS. Its components are given by Here Υ k is a normalization factor and k ⊥ denotes the momentum perpendicular to the n-BIS within the (n − 1)-BIS and pointing to h n−1 > 0. Physically, the DST describes the variation of the remaining TASP components γ l n−1 across the n-BIS. On the n-BIS we further find that the DST g nth (k) = h (n) (k), which is a key result implying that the topology on n-BIS is captured by the DST. Then the topological invariant is intuitively described by the coverage of the DST g nth (k) over the corresponding (d − n)D spherical surface S d−n formed by the unit vectorĥ (n) (k) when k runs over the n-BIS. In particular, for n = d we obtain T d- We present the quantum simulation with two realistic models to show explicitly how the high-order dBSC can provide dynamical characterization of topological phases based on dimension reduction. First we consider a 2D quantum anomalous Hall system with Hamiltonian with h x = m x + t x so sin k x , h y = m y + t y so sin k y , h z = m z − t 0 cos k x − t 0 cos k y , and t x,y so = m z = t 0 . This model has been widely realized in experiments [45,46,62,63]. We first deeply quench the h x axis from a large m x to 0 with m y = 0, and measure TASP σ x x in the whole BZ. Then the 1-BIS with σ x x = 0 can be obtained, which are two lines, k x = 0 and π, as shown in Fig. 2 (a-1). On the 1-BIS we further set m x = 0, deeply quench the h y axis, and measure σ y y . Then the 2-BIS which is formed by four momentum points can be obtained by 2-BIS = {k| σ x x = 0, σ y y = 0}, and is marked in Fig. 2 (b-1). With the same quench we could measure the TASP σ z y around the four k points and calculate the corresponding DST [defined in Eq. (6)] on the 2-BIS, as shown in Fig. 2 (b-2). We can see that the DSTs are anti-parallel on k x = 0, while they are parallel on k x = π. Accordingly, the former is nontrivial and the latter is trivial. On the whole, the 2D system is topologically nontrivial with the topological invariant T 2-BIS = 1. In order to verify the correctness of the above result of the high-order BIS, we also investigate the DSTs on the 1-BIS, as presented in Fig. 2 (a-2). While the DST exhibits a non-zero 1D winding along k x = 0, there is no winding along k x = π. The total non-zero winding number that characterizes the topologically nontrivial phase is consistent with that on the 0D 2-BIS. The validity of the present scheme can be confirmed with a different quench order. We first quench the h z axis and a new 1D BIS, which is a closed loop, is obtained through measuring σ z z , as shown in Fig. 3. Then on the BIS we quench the h x axis and measure σ x x . A 0D BIS is formed by two k points (B and D) with σ x x = 0. Through analyzing DSTs on the two BISs, we can conclude that the system is topologically nontrivial with T 2-BIS = 1, in agreement with the former result.
The characterization with high-order dBSC necessitates far less measurements of spin evolutions than that with low-order dBSC. As pointed out, in determining the information of n-BIS, measurement is required only on the (n − 1)-BIS, not in the whole BZ. This greatly simplifies the measurement strategy in real experiments.

IV. DYNAMICAL CHARACTERIZATION AND SIMULATION: SCHEME II
We proceed to study the second scheme to implement the dBSC and quantum simulation. Unlike the above dynamical scheme, the second one necessitates to measure fewer or even only a single spin component to determine the complete topological information of the system, which further optimizes the experimental studies.
We present the dynamical characterization with two steps. First, we determine the high-order BISs by measuring only the spin component γ 0 . In particular, we perform the quench along the γ 1 axis and measure γ 0 1 , which vanishes on two types of 1-BIS. One is identical to that determined by γ 0 0 = 0, and another is identical to that by γ 1 1 = 0, as introduced in the characterization scheme I. We keep only the later one. Then we apply the quench along γ 2 axis and measure γ 0 2 , which also leads to two types of 1-BIS, respectively identical to γ 0 0 = 0 and γ 2 2 = 0. We again keep the later one, which intersects with the 1-BIS identical to γ 1 1 and gives the 2-BIS. With this dimension reduction process, we can obtain the general nth order BIS by n-BIS = {k| γ 0 i = 0, i = 1, 2, . . . , n} − {k| γ 0 0 = 0}. Then, we measure the remaining spin components γ l to determine the DST as Compared with Eq. (6), it shows qualitatively the same physics that is the variation of the remaining TASP components, including the 0th one, across the n-BIS. With the results of g nth (k) the topological invariant on the n-BIS is simply given by T n-BIS [g]. In this characterization the d − n + 1 components of spin-polarization need to be measured. Similarly, for the case n = d, we determine the topology by g 0 (k) on the d-BIS, T d-BIS [g] = 1 2 d-BIS j sgn g dth 0,L j − sgn g dth 0,R j , corresponding to the coverage of the DST g dth 0 (k) over the 0D spherical surface S 0 , i.e. the 0D d-BIS. Only the γ 0 component needs to be measured in determining the BISs and dynamical topological invariant in all d times of quenches.
We present dynamical simulation with this approach for the 3D system described in Eq. (8), with only the γ 0 -component being measured in each quench. Fig. 5 (a) shows the common 1-BIS 0 with γ 0 0 = 0, as well as the other one identical to γ 1 1 = 0. They emerge in quenching the h 1 axis and measuring γ 0 1 . Keeping the latter yields the 1-BIS, which are two planes with k x = 0 and π [ Fig. 5 (b)]. Further, by respectively quenching h 2 and h 3 axes, and measuring γ 0 2 and γ 0 3 in the same way, we eventually obtain the 2-BIS and 3-BIS for the 3D system. The result shows that the two BISs are re- Advantages of characterization with high-order BISs.-We make comments on the above dynamical schemes. For the characterization based on the nth-order BIS, the total number of measurements is d+1 to determine the BIS and corresponding topological invariant in the both dynamical schemes. Nevertheless, the measurement strategy is largely simplified (optimized) for characterization with higher-order (highest-order) BISs. We compare the two extreme cases with n = 1 and n = d. In the former case one takes a single quench and measurement to determine 1-BIS, but measures TASPs for the d different spin components near the (d − 1)D 1-BIS to obtain the corresponding DST g 1st (k). In the latter one takes d times quenches to determine d-BIS, but measures TASPs along a single direction [within the (d − 1)-BIS] near several discrete points of d-BIS, which is much simplified. Moreover essentially, during the dimension reduction for n = d the number of momentum points necessitating measurement decreases very rapidly to extract the high-order BIS information, so that actual detection is remarkably simplified compared with that for n = 1. For the second dynamical scheme, only the γ 0component along a single direction is measured for n = d to determine the dynamical topological invariant, which might be useful for real experiments. In general quantum simulation of the dynamical scheme with d-BISs is optimized.

V. EXPERIMENTAL DEMONSTRATION WITH QUANTUM SIMULATOR
We present now the experimental study to verify the advantages of high-order BISs in topological characterization. As the study is carried out in momentum space, we can measure the high-order dynamical bulk-surface correspondence by emulating the BZ with one set of parameters for one momentum. We build the quantum simulator from solid-state spins of NV center in diamond [50,65] and simulate the Scheme I. A triplet ground state (S = 1) of the electrons around the center and the intrinsic nitrogen-14 ( 14 N) nuclear spin (I = 1) form a coupled system [see Fig. 6 (a)]. The Hamiltonian reads whereŝ z andî z respectively denote electron and nuclear spin operators. A magnetic field of 446 Gauss is applied along the NV's symmetry axis. It allows optical nuclear polarization [66] and leads to an electron Zeeman splitting α e = 1250 MHz and an nuclear one α n = 137 kHz. Besides, the other parameters β e = 2.87 GHz, β n = −4.95 MHz, and λ e−n = −2.16 MHz are the electronic zero-field splitting, the nuclear quadrupolar interaction ,and the hyperfine interaction, respectively. The two-qubit subsystem is formed by the subspace of {m S = 0, −1} ⊗ {m I = +1, 0}, relabeled as {|0 , |1 } S ⊗ {|0 , |1 } I , on which the Pauli operators σ and τ are defined. We apply a microwave pulse to produce a driving field of strength ω 0 and obtain an effective Hamiltonian H RWA = 2π( λ e−n 4 σ z ⊗ τ z + ω 0 cos ϕσ x ⊗ 1 − ω 0 sin ϕσ y ⊗ 1) under the rotating-wave approximation, where ϕ is the phase of the microwave pulse. A τ y -rotation of angle θ driven by a radio-frequency pulse transforms the Hamiltonian into where the magnitudes of all the coefficients can be rescaled by an overall factor based on the evolution time (see details in Appendix B). From the above experimental operations, we can map the momentum-space parameters k to the experimental parameters (ω 0 , ϕ, θ) and successfully emulate with H eff the 3D chiral topological insulator given by the Hamiltonian H 3D [Eq. (8)]. We measure the quench dynamics to determine the bulk topology. The quantum circuit for the simulation experiment is plotted in Fig. 6 (b) and simulated parameters are set as m 0 = 1.4t 0 and t so = 0.2t 0 in Eq. (8). There are three main steps in the whole experimental procedure. (i) We prepare an initial state |00 and rotate the polarization direction to "0" by a unitary control. Then, quench m 0 and the state evolves under the Hamiltonian H eff , emulating H 3D . Meanwhile, measure the TASP γ 0 0 in the BZ to obtain the 1-BIS, which is a closed surface [50]. (ii) Repeat the first step to obtain the 2-BIS and the 3-BIS through quenching and measuring along "1" and "2" directions, but the measurements are only taken on the 1-BIS and the 2-BIS, respectively, much less than whole BZ. In Fig. 6 (c-e) we can see that the former with γ 1 1 = 0 is a closed loop in the k x = 0 plane and the latter with γ 2 2 = 0 includes only two k points. (iii) The topological invariant is finally determined by the DST g 3rd 3 (k) given from the variation of γ 3 2 across the two points of the 3-BIS [Eq. (6)]. The experimental results are in good agreement with theoretical calculations in Fig. 6 (f) and the simulated system is topologically nontrivial due to the anti-parallel DSTs on the 1D 2-BIS. Compared with previous studies where one needs to measure TASPs for three different spin components near the 2D 1-BIS to obtain the corresponding DST [50,51], which contains a large number of momentum points, here one only needs to measure TASPs along a single direction near two single points of 3-BIS. Even though we also need another step to reduce the dimension of BISs, the total number of actual measurements is much less and thus the strategy is greatly simplified.
The above experimental study verifies that the characterization and simulation of topological phases using high-order BISs are of clear advantages and high efficiency. While this verification does not reply on the construction of a real topological band for the single-particle regime, the spin-qubit simulator can be further applied to investigate the interacting effects if taking into account more qubits. The interacting effects on the quench dynamics can be simulated by simultaneous control of the different spin qubits during the evolution. The topological characterization based on high-order BISs could provide optimal schemes to explore the correlated topological quench dynamics in experiment, which shall be considered in our next work.

VI. CONCLUSION AND OUTLOOK
We have proposed the concept of high-order BISs, with which we developed a new dynamical theory to characterize and simulate topological quantum phases, and experimentally built up a quantum simulator using spin qubits to verify the feasibility of the new theory in the dynamical characterization through emulating different momenta one by one. Especially, for the highest-order BIS with zero dimension, the topological invariants can be determined with the minimal measurement strategy, showing the optimal scheme with fundamental advantages in characterizing and simulating topological phases, as verified in our simulation experiment.
Theoretically, the present dimension reduction through high-order BISs is performed by reducing the degree of freedom corresponding to the Clifford algebra space, which is similar to, but different from that in typical classification theories like K-theory [67] for the topological phases. The dimension reduction in the K-theory is performed directly in the real or momentum space. While so far the dynamical theory is developed for topological phases with integer invariants, the dimension reduction approach proposed here shows important insight into the feasibility of establishing the dynamical classification of the complete set of topological phases [7,53,54] in the AZ ten-fold symmetry classes [68] and further crystalline topological phases [7], which is stimulating and shall be presented in our next works. 2016ZT06D348), Guangdong Provincial Key Laboratory Topological invariants of band insulators are usually defined in the whole or half BZ [4? , 5]. The bulk-surface duality [32] shows that the integer topological invariant of a d-dimensional (dD) topological phase can be generically mapped to the topological number defined on the (d − 1)D BIS [32]. Here we propose the concept of high-order BISs by introducing a dimension reduction approach, and show that the topological invariant of a dD topological phase can be reduced to the lower-dimensional topological numbers on arbitrary high-order BISs, rendering the high-order bulk-surface duality.
We start from a general dD gapped phase within the Z classification. Its corresponding Hamiltonian is given in Eq. (1). Here, for the convenience of following description, we choose any component of the Hamiltonian as h 0 (k) and the other components can be regarded as a spin-orbit field h so,i (k) with i = 1 ∼ d. As well known, the topology of a system remains unchanged under any continuous deformation that does not close the energy gap. Thus, we can perform a general analysis and only consider the spin-orbit field near the BIS with zero field in the other region. For this purpose, we transform the Hamiltonian into such a form, whereĥ so,i (k) = h so,i (k) / d j=1 h 2 so, j (k) makes the Hamiltonian be normalized and α k is a step function, as illustrated in Fig. 7 (a). α k = π 2 , − π 2 , and 0, when k points are inside, outside, and on the BIS, respectively. For the Z classification, when the system is odd (even)-dimensional, the corresponding topological invariant, uniformly defined as T d , is winding number (Chern number) [56, 69? ]. Before performing a continuous dimension reduction of the BIS, we first investigate and calculate the topological invariants in general odd and even systems, respectively. ±Y represents rotation of nuclear spin about ±y axis by an angle θ. e −iHt corresponds to evolution of the system under the rotating-wave approximation for a time duration of t eff . The net effect is equivalent to the evolution under H 3D for a duration of t on the right-hand side. (b-c) Sequences to measure populations. The operation labeled γ i (pink) is a deep quench along i direction, which requires unitary operations to transform the state |00 > to an eigenstate of the pre-quench Hamiltonian. e −iH 3D t (purple) corresponds to the post-quench dynamical evolution, whoes detailed quantum cirtuit is plotted in (a). The operation γ j (azure) corresponds to measurement of the γ j component, which requires an operation to transform γ j to the z basis. π MW and π RF respectively represent the microwave and ratio-frequency pulses, which are applied to rotate spins. PLM (green) corresponds to a photoluminescence measurement. Idle corresponds to a waiting time equal to the total time of the γ i , H 3D and γ j steps.
For a (2n+1)D system, its topological invariant is the winding number, defined as where the Hamiltonian has a chiral symmetry γ = i n+1 2n+1 i=0 γ i and Tr γH(dH) 2n+1 (A3) Here 2n+1 l=0 dh l = dh 0 ∧ · · · ∧ dh 2n+1 with l i. The winding number can be simplified as Substitute the above deformed Hamiltonian Eq. (A1) into the wedge product term, where k ⊥ denotes the momentum perpendicular to the BIS, j corresponds to different sectors of the BIS, and Finally, the winding number of the (2n + 1)D system, defined on the 2nD BIS, reads We have changed the range of the index i from 1 ∼ 2n + 1 to 0 ∼ 2n. One can see that the integrated area of the winding number has been changed from the (2n + 1)D BZ to the 2nD BIS. It manifests a bulk-surface correspondence. For a 2nD system, its topological invariant is the nth-order Chern number, defined as where Tr H(dH) 2n = (2n)!(−2i) n The Chern number can be calculated in the same way as the winding number. We can obtain Similarly, its integrated area is also changed from the 2nD BZ to the (2n − 1)D BIS with reducing one dimension. Interestingly, through comparing the right-hand sides of Eqs. (A6) and (A9), which are respectively the winding number defined on the BIS of the (2n + 1)D system and the Chern number defined in the BZ of the 2nD system, we find that they have a similar form. (i) Both of them have the same constant factor. (ii) In the former the summation over i includes 2n + 1 components of the (2n + 1)D system, while in the latter it includes all the components of the 2nD system. It is important that the number of the components is equal to each other. (iii) In the former the integrated area is the BIS which are 2nD, while in the latter it is the BZ which is also 2nD. Therefore, when the BIS of the (2n + 1)D system is regarded as the BZ of a new 2nD system, the two equations are completely equivalent.
Similarly, we compare Eq. (A10) and the winding number of a generic (2n − 1)D system, which can be obtained from Eq. (A4) and reads We can see that they are also equivalent. Based on the above analysis and calculations of general systems, a continuous dimension reduction of the BIS can be performed, as illustrated in Fig. 7 (c). In each process of the dimension reduction, we regard the obtained BIS as a new BZ for next dimension reduction. Repeating this process, the BIS can be reduced to (d − n) dimensions for a dD system. It is called the nth-order BIS, defined as where h i α 's are n arbitrary components chosen as h 0 during the dimension reduction. The corresponding topological invariant, which is equal to the original invariant T d , is given as where Ω 1 = α n + β n , Ω 2 = 0, Ω 3 = −α e + β e + α n + β n − λ e−n , and Ω 4 = −α e + β e .
To simulate γ 1 and γ 2 terms in H 3D [Eq. (8)], a microwave pulse of frequency ω 0 = −α e + β e − λ e−n /2 is applied to couple basis states, including |00 ↔ |10 and |01 ↔ |11 . Its corresponding interacting Hamiltonian reads (B2) Then we transform the total Hamiltonian H tot = H 0 + H int to the rotating frame defined by the microwave field. Under proper rotating-wave approximation the system Hamiltonian can be written as (B3) We further apply a rotation U rot = e −iθτ y /2 to the system Hamiltonian, and then obtain the effective Hamiltonian Eq. (11). At the same time, we imply θ = arctan(h 3 /h 1 ), and define the effective time as a rescale of the simulation time t, i.e., t eff = κt. Our purpose is to reproduce the same effect as U 3D = e −iH 3D t with the simulated evolution U eff = e −iH eff t eff . Thus, H 3D = κH eff . The parameters satisfy To show the mapping relation between the momentum space parameters and experimental ones, we substitute h 0,1,2,3 in the above equations with momentum space parameters, giving θ = arctan |t so sin k z | m 0 −t 0( cos k x +cos k y +cos k z ) , ω 0 = |λ e−n | 4 t 2 so( sin 2 k x +sin 2 k y ) [m0−t0(cos k x +cos k y +cos k z )] 2 +t 2 so sin 2 k z , ϕ = arctan sin k y sin k x , For a given set of momentum space parameters in H 3D , we can obtain corresponding experimental parameters to emulate the same system evolution over time t eff .
In the quantum simulation experiment the initial state of the system is an eigenstate of γ 0 = σ z τ z , |00 , which is obtained by applying a green laser pulse. It corresponds to a deep quench along γ 0 . Deep quenches along other axes are realized by either applying a microwave or radio-frequency pulse to prepare the system onto the eigenstates of γ 1,2,3 , as illustrated in Fig. 6 (b). For post-quench evolution, we plot the corresponding experimental circuit in Fig. 8 (a). There are three steps in the post-quench operations: (i) Rotate the nuclear spin along −y axis for an angle θ. (ii) Apply a microwave pulse with a driving strength of ω 0 and a phase of ϕ during the evolution time of t eff . (iii) Rotate back the nuclear spin along y axis for the same angle θ. The net effect of this whole process is identical to the evolution of the system under H 3D during an evolution time of t.

Time-averaged spin polarization
In the experiment all spin polarizations are finally measured in z basis of electron and nuclear spins, where a population measurement p |i, j with i, j = 0, 1 can be obtained through the optical readout. For the case of γ 0 = σ z τ z , which is already in the z basis, the corresponding spin polarization is essentially p |00 − p |01 − p |10 + p |11 . When we measure the spin polarization γ 1(2) = σ x(y) ⊗ 1, we apply a π/2 rotation on the electron spin about −y and x axes to map the σ x and σ y components to σ z , respectively. The corresponding spin polarization is given by p |00 + p |01 − p |10 − p |11 . Similarly, for the γ 3 = σ z τ x readout, a π/2 rotation on the nuclear spin about −y axis transforms the γ 3 readout to the γ 0 readout. These operations are depicted in Fig. 6 (b).
The population readout needs to record the photoluminescence (PL) photon count of the spin state, which is the average of all four levels weighted by their populations (n tot = n 1 p |00> + n 2 p |01> + n 3 p |10> + n 4 p |11> ). In order to obtain the populations, we apply different RF and MW pulses to produce equations of their different linear combinations and solve them. In Fig. 8 (b) and (c) we show the sequences of the population measurement, and the corresponding equations are given as n 1 n 2 n 3 n 4 n 3 n 4 n 1 n 2 n 2 n 1 n 3 n 4 n 3 n 4 n 2 n 1 (B6) We measure and average the spin polarization γ i (k) under quenching along the j direction over a series of time to obtained the time-averaged spin polarization γ i (k) j . For different polarizations and quench directions, we choose the same simulation time to keep consistency in experiments. The time range is chosen from 0 to the maximum The shot noise of the optical readout leads to the dominant error in our quantum simulation, and yields a normal distribution of the photon counts with a mean of N and a standard deviation of √ N. We set N between 2000 and 3000 for a fixed 20000 repetitions of each sequence.