Simulating Chern insulators on a superconducting quantum processor

The quantum Hall effect, fundamental in modern condensed matter physics, continuously inspires new theories and predicts emergent phases of matter. Here we experimentally demonstrate three types of Chern insulators with synthetic dimensions on a programable 30-qubit-ladder superconducting processor. We directly measure the band structures of the 2D Chern insulator along synthetic dimensions with various configurations of Aubry-André-Harper chains and observe dynamical localisation of edge excitations. With these two signatures of topology, our experiments implement the bulk-edge correspondence in the synthetic 2D Chern insulator. Moreover, we simulate two different bilayer Chern insulators on the ladder-type superconducting processor. With the same and opposite periodically modulated on-site potentials for two coupled chains, we simulate topologically nontrivial edge states with zero Hall conductivity and a Chern insulator with higher Chern numbers, respectively. Our work shows the potential of using superconducting qubits for investigating different intriguing topological phases of quantum matter.


Introduction
Topological phases of matter [1,2], classified beyond Landau's symmetry-breaking theory, have been attracting growing interest in recent decades. It started with the discovery of the two-dimensional (2D) integer quantum Hall effect (QHE) [3], which arises from the topological nature of Bloch bands characterised by the Chern number [4]. Topological band theory provides a direct link between theory and experiments, which are successful for identifying salient characteristics of topological states and predicting new classes of topological phases [5]. The existence of robust edge states is deeply related to the topology of gapped bulk band structures, which is the ubiquitous bulk-edge correspondence in topological systems.
Since exploring higher-dimensional physics requires an exponentially growing number of qubits, there is growing interest in creating synthetic dimensions to construct a higher-dimensional lattice in a lower-dimensional system with its internal degrees of freedom [6][7][8][9]. Moreover, a lower-dimensional topological charge pump shares the same topological origin as higher-dimensional topological physics, e.g., a one-dimensional (1D) Thouless pump for the 2D integer QHE [10] and a 2D topological pump for a four-dimensional quantum Hall system [11,12].
Here we observe several topological signatures of 2D and bilayer Chern insulators with synthetic dimensions on a programmable 30-qubit-ladder superconducting processor. We experimentally measure the band structure of the 2D Chern insulator along a synthetic dimension by analysing the temporal frequency of the system's response to local perturbations. By monitoring quantum walks (QWs) of a single excitation initially prepared at an edge qubit, we demonstrate dynamical localisation of the topologically protected edge states. The measured band structure and the dynamical signatures of topological edge states together demonstrate the bulk-edge correspondence. Furthermore, we synthesise two different bilayer Chern insulators on the 30-qubit-ladder processor. Given the same periodically modulated on-site potentials of two coupled chains, we probe the nontrivial topological edge states with zero Hall conductivity. When on-site potentials of two coupled chains have opposite signs, a Chern insulator with higher Chern numbers is probed. Our results show that superconducting simulation platforms are capable for studying different intriguing topological phases of quantum matter.
With a dimensional reduction procedure [31], a 2D integer quantum Hall system can be mapped to a 1D model with a periodic parameter as a synthetic dimension. In this context, we experimentally simulate a 2D integer QHE associated with Chern insulator using 15 qubits on one leg of the ladder (i.e., qubits labelled by Q 1,↑ to Q 15,↑ ), where the on-site potential V j,↑ of each qubit is periodically modulated. This system can be described by the 1D Aubry-André-Harper (AAH) model [32,33] with a tight-binding Hamiltonian where the second index ↑ is omitted. Here, b determines the modulation periodicity, and the modulation phase ϕ corresponds to the momentum in a synthetic dimension [31], see Fig. 1b. Note that there exists hopping between next-nearest neighbour (NNN) qubits H ′ = J ′ ∥ 13 j=1 (ĉ † jĉ j+2 + H.c.) on the same leg with a strength J ′ ∥ ≃ 0.1J ∥ , see the Supplementary Note 1. This model is topologically equivalent to the lattice model of the 2D integer QHE, proposed by Hofstadter [34], where electrons hop within the 2D lattice, subjected to a perpendicular magnetic field with b being the magnetic flux (normalised to the magnetic flux quantum) threading each unit cell. In our experiments, we fix b = 1 3 , and vary ϕ from 0 to 2π to obtain various instances of AAH models with potentials by tuning the qubits' frequencies as ω j = ω 0 + ∆ 15 j=1 cos(2πbj + ϕ), with a reference frequency ω 0 /2π = 4.7 GHz (Fig. 1c).

Band structure spectroscopy
Band structures play an essential role in characterising topological phases of matter and discovering novel classes of intriguing topological materials [5]. Here, we directly measure the topological band structure of the integer quantum Hall system along the synthetic dimension using a dynamic spectroscopic technique applied in refs. [28,35,36]. This method detects quantised eigenenergies of the quantum system from the Fourier transformation (FT) of the subsequent response of the system given local perturbations. The experimental sequence of pulses of the band structure spectroscopy are shown in Fig. 1d. With 15 qubits initialised at their idle points, we place one target qubit Q j in the superposed state |+ j ⟩ = (|0 j ⟩ + |1 j ⟩)/ √ 2, using a Y π 2 pulse. Then, all qubits are detuned to their corresponding frequencies for the quench dynamics with a time t before the readout of the Q j at its idle point in theσ x andσ y bases. For each ϕ, time evolutions of ⟨σ x j ⟩ and ⟨σ x j ⟩ are recorded when choosing a target qubit, e.g., Q 8 in Fig. 2a. Figure 2b shows the squared FT magnitude |A j | 2 of the response function χ j (t) ≡ ⟨σ x j (t)⟩ + i⟨σ y j (t)⟩ for each qubit with ϕ = 2π 3 and ∆/2π = 12 MHz. With the summation of the squared FT magnitudes of all selected qubits I ϕ ≡ j |A j | 2 , the positions of its peaks clearly indicate the eigenenergies E/2π of the system for each ϕ (Fig. 2c).

Simulating 2D Chern insulators with a synthetic dimension
The topological nature of the 1D AAH model and the 2D integer QHE associated with Chern insulator can be identified from its band structure. When setting ∆/2π = 12 MHz, we plot in Fig. 2d the band structure of the 1D AAH model along ϕ with N = 15 and open boundary conditions along x-direction. As ϕ evolves, the gapless edge states (red curves) appear in two gaps between three "bulk band" regimes (dense blue curves). We experimentally map out this band structure by measuring I ϕ for ϕ ∈ [0, 2π] (Fig. 2e), which also agrees well with the numerical result by simulating the system dynamics ( Fig. 2f ). Two gapless Dirac bands within two band gaps are clearly observed in Fig. 2e, and each gap is related to a quantised Hall conductance σ = (e 2 /h)C, with the integer C being determined by the Chern number [4].
Moreover, the topological edge states predicted in the band structure can be verified in real space by observing localisation of an edge excitation during its QWs on the 15-qubit chain [37]. After the system initialisation, we excite one qubit with a X π pulse, tune all qubits to their corresponding frequencies and measure them at a time t after their free evolutions (Fig. 1e). For the topologically trivial case with ∆/2π = 0 MHz, the measured density distributions P j (t) of single-excitation QWs initialised at any position of the qubit chain show a lightcone-like propagation of information with boundary reflections [30].
When we set b = 1 3 and ∆/2π = 12 MHz, the P j (t) for the excitation initialised at a boundary qubit (Q 1 or Q 15 ) exhibits localisation for ϕ = 2π 3 , where the edge states appear in the middle of a band gap (Fig. 3a1, a3). This localised dynamical behaviour is due to the fact that the edge-excitation modes have a main overlap with the in-gap edge states which are topologically protected by the band gaps. This is also verified by the squared FT magnitudes |A 1 | 2 and |A 15 | 2 for edge qubits (Fig. 3b1,  b3), which mainly contain information of the in-gap edge states. The discontinuity of the FT signals for edge qubits results from the second and first band gaps, respectively, because the edge states localised at two boundaries have opposite propagation directions, like in the standard QHE. When a qubit away from either end, e.g. Q 8 , is excited, we observe the propagation of the excitation in Fig. 3a2 due to the absence of topological protection, and its FT signal |A 8 | 2 merely shows partial information of the bulk band (Fig. 3b2). Thus, by directly measuring the band structure and observing dynamical localisation of edge excitations, our experiments demonstrate the bulk-edge correspondence in the synthetic 2D Chern insulator.
In addition, a topological charge pump, entailing the charge transport in a 1D time-varying potential driven in adiabatic cycles [10,38], provides an alternative way to explore the 2D integer QHE. The charge transported in a pumping cycle is determined by the Chern number [39], which is defined over a 2D Brillouin zone with one spatial and one temporal dimensions. We experimentally simulate the charge pump by adiabatically varying ϕ in a 2π period starting from ϕ 0 = 5π 3 with ∆/2π = 36 MHz. Figure 3c1, c3 plot the evolutions of distributions P j (t) of an excitation initialised at the central qubit Q 8 for the forward and backward pumping schemes, respectively. In one pumping cycle, the excitation propagates through an integer number of unit cells, in this case three qubits, which is determined by the Chern number.
Following Laughlin's argument [40,41], the role of the threading magnetic flux is played by the adiabatic variation of ϕ, leading to the excitation's transport. The displacements of the centre of mass (CoM) δ x are shown in Fig. 3d, of which the slight deviations from the Chern numbers ±1 may result from the boundary reflection on the finite-size 1D qubit chain and the small fraction of excitations for the higher bands. A brief discussion of fast pumping [42] is given in the Supplementary Note 3. In comparison, there exists no excitation transport when ϕ is not pumped (Fig. 3c2, d).

Simulating bilayer Chern insulators
Next, using all thirty qubits, we simulate two different bilayer Chern insulators. The Hamiltonian of the 30qubit-ladder quantum processor is given in equation (1), where the on-site potentials on the qubits in two legs (s ∈ {↑, ↓}) are modulated as V j,s = ∆ s cos(2πbj + ϕ).
with b = 1 3 . For instances of two coupled AAH chains by varying ϕ with the same ∆ ↑(↓) /2π = 12 MHz and op-posite ∆ ↑ /2π = −∆ ↓ /2π = 12 MHz, we effectively simulate two different bilayer quantum systems [43], where the magnetic fluxes threading two layers are the same and have a difference of π, respectively. Using the above spectroscopic technique, we measure the band structures for these two cases (Figs. 4a, b), which agree well with the theoretical predictions (dashed curves) using the experimental parameters in the Supplementary Note 1.
For ∆ ↑(↓) /2π = 12 MHz with J ⊥ being comparable to J ∥ , gapless edge states are experimentally observed (Fig. 4a), even when the Chern number for half filling is predicted to be zero (see Methods). This novel topological phase results from the existence of a pair of counter-propagating chiral edge states in the second gap (Fig. 4c), which are protected by the inversion symmetry and lead to a zero Hall conductivity σ = 0 [44]. When we excite a corner qubit (e.g., Q 1,↑ ) and monitor the excitation's QWs for ϕ = 2π 3 , the time-evolved distribution P j,s shows oscillations between two corner states localised at the same boundary rung of the qubit ladder (Fig. 4f ). This dynamical behaviour can be further understood from the measured squared FT magnitudes |A 1,↑ | 2 and |A 1,↓ | 2 with respect to the corner qubits Q 1,↑ and Q 1,↓ , respectively, which share information of the same edge states (Fig. 4c).
By modulating opposite on-site potential fields with ∆ ↑ /2π = −∆ ↓ /2π = 12 MHz, we obtain a different band structure (Fig. 4b) from the one with two layers having the same magnetic flux. In this case, the magnetic fluxes threading two layers have a difference of π, and the gapless edge states are characterised by higher Chern numbers (see Methods). The QWs of the excitation initialised at a corner qubit (e.g., Q 1,↑ or Q 1,↓ ) for ϕ = 2π 3 show dynamical localisation at the initial qubit (Fig. 4g), and the measured squared FT magnitudes of two corner qubits show information of different edge states in different gaps (Fig. 4d). Thus, this synthetic bilayer quantum system with two layers having a π-flux difference presents a Chern insulator with chiral edge states, identified by higher Chern numbers. Since high-Chern-number insulators without the formation of Landau levels have attracted increasing attentions [45][46][47], our work provides a new perspective for exploring these emergent topological phases. In comparison, for ∆ ↑(↓) /2π = 0, the QWs of the corner excitation first show linear propagation and then indicate thermalisation on the qubit ladder as a nonintegrable system [29,48].

Discussion
In summary, we simulate 2D and bilayer Chern insulators with synthetic dimensions on a programmable 30qubit-ladder superconducting processor. By measuring the band structures and monitoring dynamical localisation of edge excitations, we implement the bulk-edge correspondence in the 2D Chern insulator. In addition, we synthesise two different bilayer Chern insulators with two layers having the same magnetic flux and a π-flux difference, respectively, and simulate distinct and novel topological phases. Our experiments, using a relatively large number of superconducting qubits with long coherence times and accurate readouts, show the future potential of using superconducting simulating platforms for investigating intriguing topological quantum phases and quantum many-body physics. For instance, by upgrading the ladder-type superconducting processor with tunable interactions between two qubit chains, we could observe the quantum phase transitions between topological phases with zero Chern number (Fig. 5b) and with higher Chern numbers (Fig. 5c). In addition, when tuning weak couplings between two chains and modulating b either positive or negative for two chains, respectively, the quantum spin Hall effect could be demonstrated by simulating bilayer quantum Hall systems with opposite Chern numbers [49].

System Hamiltonian
Our superconducting quantum processor can be described as a Bose-Hubbard ladder with a Hamiltonian (ℏ = 1) [28,29] whereâ † (â) is the bosonic creation (annihilation) operator, andn ≡â †â is the number operator. Here, J ∥ /2π ≃ 8 MHz and J ⊥ /2π ≃ 7 MHz denote the nearestneighbour (NN) hopping between nearby qubits on the same leg and on the same rung, respectively. Also, η is the on-site nonlinear interaction, and V j,s is the tunable on-site potential. Our device is designed to fulfil the hard-core limit |η/J| ≫ 1, and thus, the highly occupied states of transmon qubits are blockaded, which represents fermionisation of strongly interacting bosons [30]. The system Hamiltonian can then be simplified as whereĉ † (ĉ) is the hard-core bosonic creation (annihilation) operator with (ĉ † ) 2 =ĉ 2 = 0, and [ĉ † j,s ,ĉ i,r ] = δ ji δ sr .
Note that in addition to the hopping between nearestneighbour (NN) qubits, there also exist the hopping between next-nearest neighbour (NNN) qubits on different legs: and the hopping between third-nearest-neighbour (TNN) qubits on the same leg: Mapping the 2D quantum Hall model to various instances of Aubry-André-Harper chains An electron moving in a 2D lattice with a perpendicular magnetic field b is described by the integer quantum Hall model, where t x and t y are hopping strengths along the x and y-axes, respectively. We consider a periodic boundary condition in the y-direction and introduce the Fourier transformation (FT): Then, the Hamiltonian transforms to H IQH = ky H ky , with k y being the quasi-momentum in the y-direction, where +2t y N n cos(2πbj + k y )ĉ † j,kyĉ j,ky . (10) By replacing t x , 2t y , and q with J, ∆, and ϕ, we obtain the 1D Aubry-André-Harper (AAH) model where the second index has been omitted in the 15-qubit experiment.
Since we simulated a 1D tight-binding fermionic Hamiltonian (11) in our experiments, to investigate the behaviour of one excitation (a hard-core boson) effectively captures the topological property of the system. Furthermore, in the case when we only excite one qubit, the interaction term in equation (4) is zero, and the behaviour of the quantum system we simulated is not affected by the statistics of the particle (fermionic or bosonic). Since our device is designed to fulfil the hard-core limit |η/J| ≃ 25 ≫ 1, the highly occupied states of transmon qubits are blockaded, which represents fermionisation of strongly interacting bosons [30].
Thus, the dynamical behaviour of a multiple-excitation system in the 15-qubit-chain experiment is similar as a 1D fermionic system, as demonstrated in ref. [30].

Band structure spectroscopy and effect of decoherence
After initialising the selected qubits at their idle points, we prepare one target qubit Q j,s in the superposed state |+ j,s ⟩ = (|0 j,s ⟩ + |1 j,s ⟩)/ √ 2, using a Y π 2 pulse. Then, all qubits are tuned to their corresponding frequencies for the quench dynamics, and at a time t, we measure the Q j,s at its idle point in theσ x andσ y bases. For each ϕ, time evolutions of ⟨σ x j,s (t)⟩ and ⟨σ y j,s (t)⟩ are monitored. Then, we calculate the squared Fourier transformation (FT) magnitude |A j,s | 2 of the response function [28] With the summation of the squared FT magnitudes of all selected qubits I ϕ ≡ j,s |A j,s | 2 , the positions of its peaks clearly indicate the eigenenergies {E n } of the system for each ϕ.
Given |Ψ j,s (0)⟩ = n c (j,s) n |ψ n ⟩, with H|ψ n ⟩ = E n |ψ n ⟩ and c (j,s) n ≡ ⟨ψ n |Ψ j,s (0)⟩, and the evolution of the state can be obtained by solving the Schrödinger equation as We obtain the response function of the system after the initial perturbation as χ j,s (t) ∼ 2⟨Ψ j,s (0)|Ψ j,s (t)⟩ − 1 = 2 n |c (j,s) n | 2 e −itEn − 1, (14) and the FT of the response function is calculated as where the eigenenergies {E n } are indicated by the peaks of the FT signals. Furthermore, when considering decoherence with a decaying rate γ in a form (16) the FT of the response function is calculated as which indicates that the presence of decoherence increases the width of the peaks of the FT signals, and the locations of the peaks can still indicate the values of the eigenenergies. There also exists the unwanted zero frequency signal of the FT. In our experiments, to eliminate the unwanted zero frequency signal of the FT we calculate the FT of oscillations of the response function δχ(t) ≡ χ(t) − χ(t), i.e., the response function minus the its average over the time interval.
Avoiding rung-pair excitations in strongly interacting Bose-Hubbard ladders In the 30-qubit-ladder experiment, we excited either one corner qubit or a bulk qubit and monitored the quantum walks of the excitation, to study the topologically protected edge states of the bilayer quantum Hall systems. Localisation of the excitation initialised at a corner qubit and the propagation of the excitation initialised at the bulk qubit indicate the existence of a topological edge state protected by the topology of the bulk structure.
However, we do not consider to excite two corner qubits at the same edge simultaneously, because it has been experimentally and theoretically shown in refs. [29,50] that the dynamics of single-and double-excitation states have very distinct behaviours. Specifically, in the hard-core limit, there exists rung-pair localisation at the edges even for the topologically trivial case without the modulation of the on-site potentials, ∆ ↑ /2π = ∆ ↓ /2π = 12 MHz. In the centre-of-mass frame, the two-particle system can be mapped into an effective single-particle Hamiltonian, and there exists a zero-energy flat band in the hard-core limit, which is the origin of the localisation [50]. Therefore, in the 30-qubit-ladder experiment, we avoid exciting two qubits on the same rung simultaneously, when investigating the topologically-protected edge states.

Quantum charge pumping
In addition to the study of 2D topological systems, topological charge pumping provides an alternative way to explore the quantised transport with topological protection in a dynamical 1D system. The concept of topological charge pumping was first proposed by D. J. Thouless [10], and recognised a topological quantisation of charge transport in a 1D time-varying potential driven in adiabatic cycles. The charge transported in a pumping cycle is determined by the Chern number [4], which is defined over a 2D Brillouin zone with one spatial dimension and one temporal dimension.
The frequencies for all qubits are calibrated using the frequency calibration procedure as shown in the Supplementary Information. Then, by modulating the frequencies of all 15 qubits simultaneously, we slowly vary ϕ from ϕ 0 to ϕ 0 ± 2π in 1,100 ns with a relatively slow speed ∼ 1.8π/µs for the backward and forward pumping schemes, respectively. For the case of no pump, ϕ is fixed at ϕ 0 during the time evolution.
In our experiments, we measured the |1⟩-state occupation probability of each qubit at its idle point in theσ z basis for each evolution time t. To reduce the effect of the stochastic fluctuations, we maintained a fixed sample of 4,000 single-shot readouts and repeated the measurement procedure 10 times for estimating the mean values and standard deviations at each evolution time t. We obtained the average displacements of the centre of mass (CoM) as where P j (t) is the |1⟩-state occupation probability of the Q j,↑ . Note that the CoM is divided by 3, because we set b = 1 3 and there are 3 qubits in one unit cell. The experimental results can be found in Fig. 3 in the main text.

Characterisation of bilayer Chern insulators
The Hamiltonian of the 30-qubit ladder reads where s ∈ {↑, ↓}, and b = 1 3 determines the modulation periodicity. The typical hopping strengths for our sam- The AAH ladder can be mapped into two coupled Hofstadter lattices, subjected to the same effective magnetic fields for each layer.
We now proceed to discuss the topological properties in two cases, i.e., They correspond to the study of two types of bilayer quantum systems with two layers having the same magnetic flux and a π-flux difference, respectively. In our experiments, we set ∆/2π = 12 MHz.
Two qubit chains with the same periodically modulated on-site potentials. For ∆ ↑ = ∆ ↓ = ∆, two identical AAH chains are coupled to form an AAH ladder. In Fig. 5, we plot the band structures for different inter-chain hopping strengths J ⊥ . Figure 5 shows the energy spectra of different topological phase regimes, showing that topological phase transitions occur as J ⊥ varies.
We can identify the topologically nontrivial and trivial bands by using the Chern number of each band [4], which is defined as The nth band Berry connection A n is written as where |φ n (k, ϕ)⟩ is the nth band's wavefunction, and γ = k, ϕ. Then, the Hall conductivity reads σ = n C n , with n being summed over the occupied bands. When the inter-chain hopping strength J ⊥ is much smaller (Fig. 5a) or much larger (Fig. 5c) than the intrachain hoping strength J ∥ , the topological boundary states are well characterised by the Chern number. However, when J ⊥ is comparable to J ∥ , topological edge states with zero Hall conductivity appear ( Fig. 5b) for the half filling.
The zero Hall conductivity results from the contribution of a pair of counter-propagating chiral edge states. This novel topological phase was also studied in a dimerised Hofstadter model [44] and has never been experimentally observed before. Figure 6 shows the density distribution of a mid-gap state for the half filling, where the in-gap state occupies the end sites of both chains.
The Hamiltonian H(k) in equation (22) has the inversion symmetryP for ϕ = 2π 3 and ϕ = 5π 3 , where the inversion symmetry operatorP iŝ Then, we define an integer invariant N to characterize this novel topological phase, which is expressed as [44] where N 1 and N 2 are the number of negative parities (by applying the inversion symmetry operator to eigenstates) at the high symmetry points k = 0 and k = π, respectively. Figure 7 shows the Bloch band structures for ϕ = 2π 3 and ϕ = 5π 3 , respectively, and the topological invariant N at each band gap is indicated. For ϕ = 2π 3 , we find N = 1 for 1 4 -filling and half-filling, indicating a pair of edge states pinned at ϕ = 2π 3 . For ϕ = 5π 3 , we have N = 1 for 3 4 -filling and half-filling, indicating a pair of edge states pined at ϕ = 5π 3 . Therefore, two pairs of edge states appear at half-filling, which propagate in opposite directions at each edge.
Two qubit chains with opposite periodically modulated on-site potentials. For ∆ ↑ = −∆ ↓ = ∆, the AAH ladder shows two typical topological band structures, as shown in Fig. 8. This bilayer structure, with two layers having a π-flux difference, shows different topological features from the bilayer structure with two layers having the same magnetic flux. Furthermore, the in-gap states mainly occupy the end site of one of the chains, as shown in Fig. 9. Thus, this synthetic bilayer quantum system, with two layers having a π-flux difference, presents an integer quantum Hall effect with chiral edge states, identified by higher Chern numbers.

Data availability
The source data underlying all figures are available at https://doi.org/10.6084/m9.figshare. 23925009. Other data are available from the corresponding author upon request.

Code availability
The codes are available upon reasonable request from the corresponding author. Qubit number, j d e π Q 1, Q 2, Q 3, Q 15, Q 14, Q 13, Q 4, Q 5, Q 6, Q 7, Q 8, Q 9, Q 10, Q 11, Q 12, Energy, Energy, Energy, Edge state Figure 2. Band structure spectroscopy of the 2D Chern insulator with a synthetic dimension. a, Typical data of ⟨σ x ⟩ and ⟨σ y ⟩ versus time t when choosing Q8 as the target qubit for b = 1 3 , ∆/2π = 12 MHz and ϕ = 2π 3 . b, Squared FT magnitudes |Aj| 2 of the response functions χj(t) ≡ ⟨σ x j (t)⟩ + i⟨σ y j (t)⟩ for all fifteen qubits. Note that the edge states at the edge qubits: Q1 and Q15. c, Summation of the squared FT magnitudes I ϕ ≡ j |Aj| 2 . d, Band structure of the 2D Chern insulator for b = 1 3 and ∆ = 12 MHz with 15 sites along the x-direction and periodic along the y direction. Red curves show the edge states. e, f, Experimentally measured data for I ϕ (e) for different values of ϕ varied from 0 to 2π are compared with numerically calculated data of I ϕ (f ) obtained by numerically simulating the dynamics of the 15-qubit system without decoherence.  , ∆/2π = 12 MHz and ϕ = 2π/3 after initially exciting the leftmost qubit Q1 (a1), the central qubit Q8 (a2), and the rightmost qubit Q15 (a3). b1-b3, Experimental data of the squared FT magnitudes |A1| 2 , |A8| 2 , and |A15| 2 , when choosing Q1 (b1), Q8 (b2), and Q15 (b3) as the target qubits. c1-c3, Time evolution of an excitation initially prepared at the central qubit Q8 when it is forward pumped (c1), not pumped (c2) and backward pumped (c3), respectively, with ∆/2π = 36 MHz for an initial ϕ0 = 5π/3. d, Displacement of the centre of mass (CoM) δx versus time t in one pumping cycle with period T for the cases in (c1-c3). The error bars are 1 SD, calculated from all 10 groups of experimental results.     Fig. 8b are not clearly shown due to finite-size effects. Our experiments are performed on a superconducting circuit consisting of 30 transmon qubits (Q j,s , with j varied from 1 to 15 and pseudo-spin s = {↑, ↓}), which constitute a twolegged qubit ladder, as shown in Supplementary Fig. 1. Each qubit, coupled to an independent readout resonator (R), has an independent microwave line for XY and Z controls. This ladder-type 30-superconducting-qubit device is fabricated on a 430 µm thick sapphire chip with standard wafer cleaning.

Supplementary
First, a 100 nm thick layer of Al is deposed on a 10 × 10 mm 2 sapphire substrate, which is patterned with optical lithography using a 0.70 µm thick layer of positive SPR955 resist. Then, we apply the wet etching method to produce the large components of the superconducting chip, such as the microwave coplanar waveguide resonators, the transmission lines, the control lines, and the capacitors of the transmon qubits.
The second step is the preparation process for Josephson junctions. After patterning a bilayer of MMA and PMMA resists with electron beam lithography, the Josephson junctions are made using double-angle evaporations, which include a 65 nm thick Al layer at +60 • followed by an oxidation process in pure oxygen for several minutes and then a 100 nm thick second layer of Al at 0 • . Finally, in order to suppress the parasitic modes, several airbridges [1] are constructed on the chip.

B. Wiring
A schematic diagram of our experimental setup and wiring is shown in Supplementary Fig. 2. In our experiments, the superconducting quantum processor is cooled in a dilution refrigerator (BlueFors, XLD1000sl) with a base temperature of mixing chamber (MC) about 15 mK. This device has three readout lines, and the readout signals are amplified by the Josephson-based parametric amplifiers (JPAs) and the high- Supplementary Fig. 1. Optical micrograph of the thirty-qubit sample. Thirty superconducting qubits constitute a ladder. Each qubit has an independent readout resonator and an independent microwave line for XY and Z controls. Thirty readout resonators are separated into three groups, and ten resonators in each group share a microwave transmission line (highlighted in black, blue or red) for measurements.
electron-mobility transistors (HEMTs). The specific arrangement of the attenuators and filters are used for noise suppression. The up-conversion principle is used to generate the XY and read-in signals, while the down-conversion principle is for readout signal. The Z signals are directly generated by the arbitrary waveform generator (AWG). The XY and Z signals are coupled together at the MC plate of the dilution refrigerator by using a directional coupler.

C. Device information
For the 15-qubit experiment, fifteen superconducting qubits (Q j,↑ with j from 1 to 15) on one leg of the qubit-ladder are used, which constitutes a one-dimensional (1D) qubitchain. The idle frequencies of 15 qubits, denoted by ω idle j,↑ /2π, are carefully arranged to minimise the crosstalk error among qubits during the single-qubit operations, while the 15 unused qubits (Q j,↓ with j from 1 to 15) on the other leg are detuned far from 3.5 GHz to avoid any unwanted interaction. The qubits characteristics for the 15-qubit experiment are shown in Supplementary Table 1.
For the 30-qubit experiment, all thirty superconducting qubits (Q j,s with j from 1 to 15 and s ∈ {↑, ↓}) are used. The qubit's characteristics for the 30-qubit experiment are shown in Supplementary Table 2. As shown in Supplementary Table 2, the decoherence times T 1 and T 2 * are much longer than the system's evolution time within 1 µs. ) on the same leg. ω idle j,↑ /2π is the idle frequency of the qubit, at which the single-qubit gates are performed.T 1 j,↑ denotes the average energy relaxation time from the frequency about 4.0 GHz to the maximum frequency ω max j,↑ /2π of Qj. Also, ω R j,↑ /2π is the readout resonator's frequency of Q j,↑ with g j,↑ /2π being the coupling strength between the qubit and the readout resonator. The anharmonicity η j,↑ /2π of Q j,↑ is defined by η j,↑ ≡ ω 10 j,↑ − ω 21 j,↑ , with ω 10 j,↑ /2π (ω 21 j,↑ /2π) being the energy difference between the |1⟩ and |0⟩ states (|2⟩ and |1⟩ states). F 0 j,↑ (F 1 j,↑ ) is the measurement probability of |0⟩ (|1⟩) when the qubit is prepared at |0⟩ (|1⟩), which is used to mitigate the readout errors.

D. System Hamiltonian
Our superconducting quantum processor can be described as a Bose-Hubbard ladder with a Hamiltonian (ℏ = 1) [ whereâ † (â) is the bosonic creation (annihilation) operator, andn ≡â †â is the number operator. Here, J ∥ /2π ≃ 8 MHz and J ⊥ /2π ≃ 7 MHz denote the nearest-neighbour (NN) hopping between nearby qubits on the same leg and on the same rung, respectively. Also, η is the on-site nonlinear interaction, and V j,s is the tuneable on-site potential. Our device is designed to fulfil the hard-core limit |η/J| ≫ 1, and thus, the highly occupied states of transmon qubits are blockaded, which represents the fermionisation of strongly interacting bosons [4]. The system Hamiltonian can then be simplified as whereĉ † (ĉ) is the hard-core bosonic creation (annihilation) operator with (ĉ † ) 2 =ĉ 2 = 0, and [ĉ † j,s ,ĉ i,r ] = δ ji δ sr . Note that in addition to the hopping between nearestneighbour (NN) qubits, there also exist the hopping between next-nearest-neighbour (NNN) qubits on different legs: and the hopping between third-nearest-neighbour (TNN) qubits on the same leg:

E. Qubits coupling strengths
The coupling strengths between the nearest-neighbour (NN), next-nearest-neighbour (NNN), and third-nearestneighbour (TNN) qubits are shown in Supplementary Fig. 3a and 3b for the 15-qubit and 30-qubit experiments, respectively. The coupling strengths between the selected qubits are measured at 4.7 and 4.5 GHz for the 15-qubit and 30-qubit experiments, respectively. For the 15-qubit experiment, we only measured the nearest and next-nearest coupling strengths, which are considered in the numerical simulation, since the TNN coupling in the 15-qubit chain is too weak and can be neglected. For the 30-qubit experiment, we measured the NN, NNN and TNN coupling strengths, which are used in the numerical simulation.

F. Single-qubit gates
In our experiments, we initialise a qubit from the |0⟩ state to the |1⟩ and (|0⟩ + |1⟩)/ √ 2 states by applying a X π and Y π 2 gates on this qubit for quantum walks and topological . ω idle j,s /2π is the idle frequency of the qubit, where single-qubit gates are performed.T 1 j,s denotes the average energy relaxation time from the frequency about 4.0 GHz to the maximum frequency ω max j,s /2π of Qj. The dephasing time T * 2 is measured at the idle frequency ω idle j,s /2π. Also, ω R j,s /2π is the readout resonator's frequency of Qj,s with gj,s/2π being the coupling strength between the qubit and the readout resonator. The anharmonicity ηj,s/2π of Qj,s is defined by ηj,s ≡ ω 10 j,s − ω 21 j,s , with ω 10 j,s /2π (ω 21 j,s /2π) being the energy difference between the |1⟩ and |0⟩ states (|2⟩ and |1⟩ states). F 0 j,s (F 1 j,s ) is the measurement probability of |0⟩ (|1⟩) when the qubit is prepared at |0⟩ (|1⟩), which is used to mitigate the readout errors. band structure measurements, respectively. Single-qubit operations, including the Y π 2 and X π gates, have a duration of 120 ns and a full width at half maximum (FWHM) of 60 ns. By using the derivative removal by adiabatic gate (DRAG) theory, the quadrature correction terms with a DRAG coefficient are optimised to minimise the leakage to higher energy levels [5]. The readout pulses for all qubits have a duration of 2.0 µs, and the readout pulses powers and frequencies are optimised to realise high-visibility and low-error readouts. The readout fidelities, denoted by F 0 and F 1 , are shown in Supplementary Table 1 and Supplementary Table 2 for the 15-qubit and 30-qubit experiments, respectively.
Here, to identify the accuracy of the initial state preparation, we calculate the fidelity of the prepared states compared with the perfect |1⟩ and (|0⟩ + |1⟩)/ √ 2 states, respectively, by using the quantum state tomography (QST) technique. The QST measurements of the prepared states require individually measuring the qubit in bases formed by the eigenvectors of σ x ,σ y , andσ z , respectively. In our experiments, the mea-surement along the z-axis of the Bloch sphere can be directly performed, and those along the x-and y-axes are realised by inserting the Y − π 2 and X − π 2 gates on the qubit before the measurements, respectively. In total, we have three tomographic operations {I, X − π 2 , Y − π 2 } and obtain two occupation probabilities {P 0 , P 1 } for each operation, which allow us to reconstruct the density matrix ρ exp of the initially prepared state.
The fidelity between the experimentally measured density matrix ρ exp and the ideal one |ψ⟩ can be quantified as F (ρ exp , |ψ⟩) = ⟨ψ|ρ exp |ψ⟩. We maintain a fixed sample of 3,000 repetitions of readouts and repeat the measurement procedure ten times for estimating the average value and the standard deviation of the state fidelity for each qubit. Experimental results of the fidelity for the (|0⟩ + |1⟩)/ √ 2 and |1⟩ states are both above 0.990, which are shown in Supplementary Fig. 4a and 4b, respectively.     Fig. 4. Single-qubit quantum state tomography (QST) for the (|0⟩ + |1⟩)/ √ 2 state (a) and the |1⟩ state (b), which are prepared by applying Y π 2 and Xπ gates, respectively.

A. Pulse amplitudes for Z and XY controls
In our experiments, we periodically calibrate the amplitude of the Z pulse (about 10 µs) at the idle frequency and the pulse amplitude of the single-qubit gate on each qubit to avoid the variations of the performances of the Z and XY controls with time.

B. Readout probability
We monitor the time evolution of the readout fidelity during the experiment in order to perform the real-time correction of the readout errors [6], by applying the inverse of the tensor product of the matrices for the selected qubit Q j,s to the measured occupation probability vector (P 0 j,s , P 1 j,s ) T .

C. Z pulse distortion
In the waveform of the Z pulse to adjust the frequency of a qubit, the rising edge, the falling edge, and the flatness of the Z pulse are very significant for realising a high-fidelity quantum gate and the long-time quantum state evolution. The Z pulse distortion would cause an unwanted frequency drift before the readout process. This issue causes errors when we measure in theσ x andσ y bases, because the frequency drift would generate an unwanted phase shift that greatly affects the measured result [7]. Therefore, it is of fundamental importance to correct the unwanted Z pulse distortion. In our experiments, we develop two types of pulse sequences, which aim to correct the Z pulse distortions occurring in the short-time and longtime scales, respectively.
The existence of the DC component in the waveform of the Z pulse will shorten the dephasing time T 2 * of the qubit. In our experiments, a DC-block element (with a capacitance ∼ 1-2µF) is employed for blocking the DC component. However, the use of the DC-block element inevitably distorts the shape of the Z pulse with a long time (over 500 ns), which decreases the modulation precision of the qubits frequencies.
Thus, we design a pulse sequence for correcting the Z pulse distortion as shown in Supplementary Fig. 5a. First, a X π pulse around the target frequency and a Z pulse with a fixed duration and a fixed amplitude are applied to the qubit. Then, we vary the starting time of the X π pulse and measure the population of the |1⟩ state for different time delays. The qubit can be excited to |1⟩ when the frequency of the qubit matches that of the applied X π pulse. As shown in Supplementary  Fig. 5c, the measured population of |1⟩ before the Z correction shows the shape of the Z pulse that is sensed by the qubit. After obtaining the maximum value of the population of |1⟩ for each delay, we can fit the data of the Z pulse amplitude (Zpa) with an exponential decay function: with α and β being the fitting parameters and T 0 being the estimated input of the characteristic decay factor T d . As shown in Supplementary Fig. 5b, T d = 155.45 µs is fitted for the DC-block element with a capacitance 2 µF. With the obtained characteristic decay factor and the deconvolution method, we can send a corrected waveform to the arbitrary waveform generator (AWG). Finally, the qubit can sense a relatively flat Z pulse as shown in Supplementary Fig. 5d. The above procedure can correct the distorted Z pulse for the long-time scale. To calibrate the Z pulse distortion within a short-time scale (e.g., < 100 ns), the pulse sequence as shown in Supplementary Fig. 6a is designed to record the shape of the tail of the Z pulse. Note that the measurement is performed at the idle frequency of the qubit. Therefore, it is better to bias the qubit to a low frequency (1 GHz below the sweet point ω max j,s /2π) for obtaining a Zpa-sensitive region. The experimental result of the pulse-shape measurement is shown in Supplementary Fig. 6b, and an obvious distortion is observed before the correction. In our experiments, a finite impulse response (FIR) filter and several first-order infinite impulse response (IIR) filters are employed for the correction [8]. The FIR filter can be described by a polynomial function with 16-24 parameters, while the IIR filters are an integration of several exponential functions. The experimental result after the first round of the correction is shown in Supplementary  Fig. 6c. Generally, in order to obtain a better performance, we repeat the correction procedure twice and finally obtain a stationary step response behaviour as shown in Supplementary  Fig. 6d.

D. Z pulse crosstalk
The Z pulse crosstalks between different qubit pairs is calibrated by experimentally measuring the Z-crosstalk ma-trixM Z using a similar method as shown in Supplementary Ref. [9]. Each element ofM Z denotes the Z pulse amplitude (Zpa) that is sensed by Q j,s , when a unity Zpa (1 arb. units) is applied to Q i,r . Note that all elements (M Z ) i,r;j,s are less than 4.5%, as shown in Supplementary Fig. 7. The Zpa applied to each qubit is calibrated based on the mapping relationship where Z sensed and Z applied are for the Zpas that are sensed by the qubit and applied to the qubit, respectively.

E. Qubits frequencies
In the 15-qubit experiment, we simulate the twodimensional (2D) integer quantum Hall effect (QHE) using 15 qubits on one leg of the qubit ladder by tuning the qubits frequencies as ω j (ϕ) = ω 0 + ∆ 15 j=1 cos(2πbj + ϕ), with the reference frequency being fixed at ω 0 /2π = 4.7 GHz. Here, ϕ is varied from 0 to 2π in order to obtain various instances of the 1D Aubry-André-Harper (AAH) chains with different on-site potentials V j (ϕ) = ω j (ϕ) − ω 0 , in the rotating frame with respect to ω 0 . Then, we experimentally simulate the 1D topological charge pump, which provides an alternative way to explore the 2D integer QHE, by adiabatically varying ϕ initially from ϕ 0 = 5π 3 with ∆/2π = 36 MHz. In the 30-qubit experiment, we tune the qubits frequencies as ω ′ j,s (ϕ) = ω ′ 0 + ∆ s 15 j=1 cos(2πbj + ϕ), with the reference frequency being fixed at ω ′ 0 /2π = 4.5 GHz. Thus, it is very significant to precisely adjust the qubits frequencies in our experiments. For the single-qubit experiment, it is simple to tune the qubit to the target frequency with a  Supplementary Fig. 5. Calibration of the Z pulse distortion for the long-time scale caused by the DC block. (a) Schematic of the calibration pulse sequence. We apply a Z pulse with a fixed duration (∼ 10 µs) and a fixed amplitude (∼ 1 arb. units) and excite the qubit around the target frequency by a Xπ pulse. Then, we measure the population of the |1⟩ state for different time delays between the Xπ pulse and the Z pulse, given a difference of the frequency δf . (b) An exponential decay function is used to estimate the characteristic decay factor as Td = 155.45 µs. (c and d) Z pulse shapes obtained by measuring population of |1⟩ before (c) and after (d) the Z correction, respectively. determined relationship between the Zpa and the qubit frequency, as shown in Supplementary Fig. 8a. However, for the multi-qubit case, the target qubits frequencies will shift slightly when other qubits Z signals are applied simultaneously, even though the Z crosstalks have been calibrated by using the Z-crosstalk matrixM Z . This indicates the existence of residual Z crosstalks. Although we can obtain a larger Zcrosstalk matrixM ′ Z when other qubits Z pulses are applied, it is not very worthwhile, because to calibrate all the matrix elements are very time consuming.
Here, we perform direct Rabi oscillation measurements on each stand-alone qubit around the reference frequency [from (ω 0 /2π − 36 MHz) to (ω 0 /2π + 36 MHz)]. As shown in Supplementary Fig. 8b, to calibrate the Zpa for detuning Q 1,↑ to 4.7 GHz, we design two Zpa configurations for other qubits to cancel the effect of the residual Z crosstalks. Given the experimental on-site potentials and the coupling strengths between qubits as shown in Supplementary Fig. 3a, the nearest-neighbour (NN) qubit (Q 2,↑ ), the next-nearestneighbour (NNN) qubit (Q 3,↑ ) and other qubits are detuned ±80 MHz, ±30 MHz and ±20 MHz with respect to ω 0 /2π, respectively. Note that there are also other Zpa configurations in addition to those as shown in Supplementary Fig. 8b, which hardly affect the final result.
We fix the driving magnitude of the local transverse field with a resonant frequency at 4.7 GHz and scan the Zpa around the guessed one (marked with a blue star in Supplementary  Fig. 8a). For the two Zpa configurations, the 2D Rabi oscillation results are plotted in Supplementary Fig. 8c and 8d, respectively. The slowest Rabi oscillation characterises the optimised Zpa, corresponding to the smallest driving strength as shown in Supplementary Fig. 8e.
Finally, the Zpa used in our experiments (marked with a red star in Supplementary Fig. 8f) is fixed at the average value of the two optimised Zpa configurations. Note that the final Zpa indicates a frequency shift of about 2.5 MHz with respect to that shown in Supplementary Fig. 8a. For simplicity, we can select some typical frequencies and repeat the above optimisation procedure for obtaining a new relationship between the Zpa and qubits frequencies in the region studied (see Supplementary Fig. 8f).  Fig. 6. Calibration of the Z pulse distortion for short-time scales. (a) Schematic of the experimental pulse sequence. We apply a Z pulse with a fixed duration (∼ 1 µs) and a fixed amplitude (∼ 1 arb. units) to the qubit for obtaining the shape of the rising (falling) edge by exciting the qubit around the its idle frequency (by varying the Z offset parameter). (b) Pulse shape before the correction. (c and d) Z pulse shapes after the first (c) and second (d) correction procedures, respectively.

F. Phase calibration
As shown in Fig. 1d in the main text for the pulse sequences for the band structure spectroscopy, a Y π 2 rotation pulse is applied on the target qubit for measuring the time evolutions of ⟨σ x (t)⟩ and ⟨σ y (t)⟩. Before the readout of Q j,s at its idle point in the bases ofσ x andσ y , we detune all qubits to their corresponding frequencies for the quench evolution with a time t by applying a Z pulse to each qubit. This operation will accumulate a dynamical phase that needs to be compensated when applying the rotation pulse after the time evolution. For each qubit Q j,s , the dynamical phase accumulated during the evolution can be calculated by δω j,s t, where δω j,s ≡ ω j,s − ω 0 , and ω 0 /2π is the target frequency of Q j,s for the quench dynamics. Although the Z pulse distortion has been carefully calibrated before the experiment, there are still some imperfect factors, leading to an additional phase shift, as compared with the theoretical calculations.
Therefore, we design a pulse sequence for the phase calibration. First, a Y π 2 pulse is applied to Q j,s , and then, all qubits frequencies are arranged with two configurations, similar as the case in Supplementary Fig. 8b for the calibration of the qubits frequencies. Then, we apply another Y π 2 pulse to Q j,s and vary its microwave phase for obtaining the probability when Q j is in |1⟩ as a function of t and the phase difference (φ − δω j,s t). For these two Zpa configurations, we obtain the experimental results shown in Supplementary Fig. 9a and 9b, respectively. To calculate the phase compensation for the second rotation pulse, we perform a standard cosine fitting to the probability of a qubit Q j,s in the |1⟩ state, P j,s , as a function of (φ − δω j,s t) for each evolution time t. The results are shown in Supplementary Fig. 9c, and the final phase compensation can be obtained as (φ 1 + φ 2 )/2, where φ 1 and φ 2 are the phase compensations under two Zpa configurations, respectively. After the above phase calibration process and adding the phase compensation on each qubit, we obtain the experimental results shown in Supplementary Fig. 9d, which indicate that the phase shift is compensated for all evolution times.

Supplementary Note 3. SUPPLEMENTARY EXPERIMENTAL DATA
A. Single-excitation quantum walks After initialising the system with a state |0⟩ ⊗N , we excite one qubit Q j,s with a X π pulse, tune all qubits to their corresponding frequencies and measure them at a time t after their free evolutions. The total repeat count for the measurement at each time on each qubit is 3,000.
In the 15-qubit experiment, for b = 1 3 and different values of ∆ ↑ /2π and ϕ, the experimental results of the time evolution of the excitation probability P j,↑ are shown in Supplementary  Fig. 10, 11, 12 and 13 for: respectively, which are compared with the theoretical predictions using the parameters in Supplementary Table 1.
To evaluate the performance of the experimental results of the quantum walks, we calculate the fidelity of the evolved state using where p j,s (t) and q j,s (t) are the experimental and theoretical probability distributions of the measurements on Q j,s in the |1⟩ state. We consider the mean fidelity that is averaged for different initial states with the excitation on different qubits. As shown in Supplementary Fig. 19, the relatively high fidelity (above around 0.9 and 0.85 within 500 ns for the 15qubit and 30-qubit experiments, respectively) for all on-site potentials configurations imply that our experimental results are well consistent with the theoretical simulations.

B. Band structure spectroscopy
After initialising the selected qubits to their idle points, we place one target qubit Q j,s in the superposed state |+ j,s ⟩ = (|0 j,s ⟩ + |1 j,s ⟩)/ √ 2, using a Y π 2 pulse. Then, we detune all selected qubits to their corresponding frequencies for the quench dynamics with a time t. Then, we measure Q j,s at its idle point in theσ x andσ y bases, and the total repeat count for each measurement on each qubit is 2,400. For each ϕ, the time evolutions ofσ x j,s andσ y j,s are recorded. We then calculate the squared Fourier transformation (FT) magnitude |A j,s | 2 of the response function χ j,s (t) ≡ ⟨σ x j,s (t)⟩ + i⟨σ y j,s (t)⟩ for each qubit.
For the case ∆ ↑ /2π = −∆ ↓ /2π = 12 MHz, we choose 20 qubits for the band structure measurements, which include all four corner qubits, i.e., Q 1,↑ , Q 1,↓ , Q 15,↑ , and Q 15,↓ , and 16 bulk qubits, i.e., Q 2,↑ , Q 4,↑ , Q 5,↑ , Q 6,↑ , Q 8,↑ , Q 10,↑ , Q 12,↑ , Q 14,↓ , Q 12,↓ , Q 10,↓ , Q 9,↓ , Q 8,↓ , Q 6,↓ , Q 5,↓ , Q 4,↓ , and Q 3,↓ . In addition to the topological charge pumping data as shown in Fig. 3c1-3c3 and 3d in the main text, we have also performed experiments with a faster pumping rate, which clearly demonstrate fast pumping. In details, we have decreased the period T from 1,100 ns to 410 ns, corresponding to a faster ramping speed of ϕ. The experimental results for the displacement of the centre of mass (CoM) δx versus the time t/T are compared with numerical simulations in Supplementary Fig. 25a. It is obvious that the pumping process with a period T = 410 ns shows fast pumping with a large derivation from the Chern number +1. The experimental data of the evolutions of distributions P j (t) during the topological pumping and the displacement of the CoM with a period T = 410 ns are compared with the data with T = 1,100 ns in Supplementary Fig. 26. For fast pumping with T = 410 ns, we maintained a fixed sample of 3,000 single-shot readouts and repeated the measurement procedure 8 times (16 times for the no pump case) for estimating the mean values and standard deviations at each evolution time t. Furthermore, Supplementary Fig. 25b shows the numerical results of the displacement of the CoM in one pumping cycle with different periods, which indicates a fast pumping region in our system when T ≲ 610 ns.

D. Numerical details
Numerical simulations are performed using the QUTIP [10,11] (the quantum toolbox in PYTHON) and NUMPY. The time evolutions of the systems are numerically simulated using QUTIP's master equation solver mesolve, where the parameters in Supplementary Fig. 3 are used. Because the evolution time is much shorter than the qubits' decoherence times t ≪T 1 , T 2 * , we do not consider the effect of decoherence in simulations. Qubit number, j  Fig. 11. Time evolution of the excitation probability P j,↑ with ∆ ↑ /2π = 12 MHz, b = 1 3 , and ϕ = 2π 3 after initially exciting a qubit on the 15-qubit chain.  Qubit number, j  Fig. 13. Time evolution of the excitation probability P j,↑ with ∆ ↑ /2π = 12 MHz, b = 1 3 , and ϕ = π after initially exciting each qubit on the 15-qubit chain.