Observation of Bloch Oscillations and Wannier-Stark Localization on a Superconducting Processor

The Bloch oscillation (BO) and Wannier-Stark localization (WSL) are fundamental concepts about metal-insulator transitions in condensed matter physics. These phenomena have also been observed in semiconductor superlattices and simulated in platforms such as photonic waveguide arrays and cold atoms. Here, we report experimental investigation of BOs and WSL simulated with a 5-qubit programmable superconducting processor, of which the effective Hamiltonian is an isotropic $XY$ spin chain. When applying a linear potential to the system by properly tuning all individual qubits, we observe that the propagation of a single spin on the chain is suppressed. It tends to oscillate near the neighborhood of their initial positions, which demonstrates the characteristics of BOs and WSL. We verify that the WSL length is inversely correlated to the potential gradient. Benefiting from the precise single-shot simultaneous readout of all qubits in our experiments, we can also investigate the thermal transport, which requires the joint measurement of more than one qubits. The experimental results show that, as an essential characteristic for BOs and WSL, the thermal transport is also blocked under a linear potential. Our experiment would be scalable to more superconducting qubits for simulating various of out-of-equilibrium problems in quantum many-body systems.


Introduction
The transport phenomena in solids is one of the central topics in condensed matter physics. About 80 years ago, Bloch and Zener predicted that electrons cannot spread uniformly in a crystal lattice under a constant force, and instead, they would oscillate and localize [1][2][3]. This oscillation is called Bloch oscillations (BOs), and the corresponding localization is called Wannier-Stark localization (WSL). BOs and WSL are typical quantum effects which reveal the wave properties of electrons.
However, they can hardly be observed directly in normal bulk materials due to the requirement of long coherence times. It is not until the 1990s that these phenomena were observed experimentally in semiconductor superlattices [4]. Nevertheless, the relaxation time in this type of material is still a bottleneck for studying BOs and WSL. During the last two decades, the developments in quantum technology have made it possible to simulate these quantum phenomena in artificial quantum systems [5,6]. Compared with the semiconductor superlattice systems, these artificial quantum systems have much longer decoherence times making them suitable for the experimental study of BOs. BOs in bosonic systems have been observed in the cold atoms [2,[8][9][10][11][12][13][14][15][16] and photonic waveguide arrays [17], etc.
In this work, we experimentally investigate BOs and WSL of spin system on a 5-qubit superconducting processor. The effective Hamiltonian can be described by an isotropic XY chain. By manipulating the frequencies of superconducting qubits precisely, we can construct a linear potential. Under this type of potential, we observe that the spin can hardly propagate through the  [43][44][45][46]48]. The nearest-neighbor two qubits are coupled capacitively, and the negative anharmonicity U of each qubit gives the on-site attractive interaction. The frequency of each qubit is tunable by individual microwave driving through the flux-bias line (Z line). The qubit flipping can be realized by XY lines. Each qubit is coupled to a resonator for individual and simultaneous readout. More experimental details of the system are presented in the Supplemental Note 1. (b) The sketch of the corresponding Bose-Hubbard chain with a linear potential. The detuning of the nearest-neighbor two qubits is the potential gradient F . The red ball represents the photon (the excitation of the qubit), which can tunnel to the nearest-neighbor sites (black arrows). (c) Pulse sequences for studying the transport of spin. The qubits are ordered by their frequencies, and initialized at their idle frequencies and the state 0⟩, see Supplemental Note 1. A. Then, we excite Q1 to the state 1⟩ by a X gate and bias all qubits at the work point with square waves. After the system evolves for a specific time t, all qubits are biased back to their idle frequencies, and we finally read out each qubit. (d) Pulse sequences for studying the thermal transport. We use two X 2 gates at Q1 and Q2 two prepare the initial state X+X+000⟩, and other sequences is the same as (c).
lattice during the quench dynamics. It tends to oscillate at the vicinity of initial positions, which is a typical phenomenon of BOs and WSL. In addition, using the maximum probability of a photon propagating from one boundary to another to represent the WSL length, we can demonstrate that the localization length is inversely correlated to the potential gradient. By performing precise simultaneous readout of two superconducting qubits, we can also study the thermal transport of the system. It is shown that the energy transport is suppressed as well by the linear potential.

Results
Experimental setup and model. In this experiment, our superconducting processor contains 5 qubits arranged into a 1D chain, with the capability of high-precision simultaneous readouts and full controls, see Fig. 1 (1) whereâ † j (â j ) is the photon creation (annihilation) operator,n j ≡â † jâ j is the number operator, g j,j+1 is the nearest-neighbor coupling strength, U j < 0 is the on-site attractive interaction resulted from the anharmonicity, and h j is the local potential which is tunable by DC biases through Z lines. To realize BOs, we let h j vary linearly along the lattice sites, i.e, h j = F j, where F is the potential gradient or the detuning of nearest-neighbor two qubits, see Fig. 1 In this superconducting circuits, since U j g i,j ≫ 1 and U j is staggered to suppress higher order tunneling (see Supplemental Note 1. A), the Fock space of the photons at each qubit can be truncated to two dimensions. Thus, the model is equivalent to a spin-1 2 system, and the nonlinear term can be neglected. Therefore, the effective Hamiltonian of Eq. (1) can be reduced to an isotropic XY The experimental results for the time evolution of density distribution Pj(t) up to 300 ns. The initial state is 10000⟩, and the potential gradient (a) F 2π = 0, (b) F 2π = 5, (c) F 2π = 10,(d) F 2π = 15 MHz. The photon or the spin can exhibit a light-cone-like spreading in the system. As the increase of F , the photon becomes more difficult to propagate to the right boundary, instead, it tends to oscillate in the vicinity of the initial position. Each point shows the average of 6 × 100 single-shot measurements. (e-h) The corresponding theoretical results of (a-d). The numerical and experimental results are consistent with each other. The details of numerical method are presented in Methods.
model [30,32] whereσ ± = (σ x ± iσ y ) 2, andσ x,y,z are Pauli matrices. According to Eqs. (1-2), we know that the system has an U (1) symmetry, so that the total spins ∑ 5 j=1σ + jσ − j forĤ eff (or the total photon number ∑ 5 j=1nj forĤ) are conserved. In the following discussion, we do not distinguish the photons and spins. In addition, this system is time-independent, thus the energy is also conserved, where we can explore both spin and thermal transport in this system. Note that we only consider the evolution time t ≤ 300 ns in the experiment, which is much smaller than decoherence time (more than 17 µs, see Supplemental Note 1. A). Therefore, the above conservation laws are nearly unbroken under the impact of decoherence.
Spin transport. Firstly, we study the spin transport after a quantum quench. The explicit experiment sequences are shown in Fig. 1(c). We initially excite the leftmost qubit Q 1 from the state 0⟩ to 1⟩ by a X gate, i.e., the initial state is ψ(0)⟩ = 10000⟩. Then, each qubit is biased to the working frequency with the fast Z pulse, and the system will evolve under the Hamiltonian (2). Finally, we measure, for each qubit, the probability distribution of state 1⟩, i.e., the density distribution of the photon or spin, defined as where ψ(t)⟩ = e −iĤt ψ(0)⟩ is the wave function of the system at time t. As shown in Fig. 2(a), when F = 0, the spin displays a light-cone-like propagation without any restrictions and can exhibit a reflection when approaching the boundaries [30,32]   is the peak value of Gaussian curve shown in (a). The corresponding error bars are the fitting errors. The numerical data points are peak values of the first wavefronts of P5(t) without Gaussian fitting [30]. Here, ln P max 5 is nearly linearly related to the potential gradients F . The solid and dashed lines are the corresponding linear fittings of the experimental and theoretic results, respectively. Now we extract the oscillation amplitude or WSL length ξ W S . Generally, in the presence of linear potential, the single-particle wave function is localized and has the form ψ(x) ∼ Ae −x ξ W S , where A is a normalized factor. Hence, the probability that a particle can propagate the distance r, i.e., P (r), should satisfy P (r) ∝ e −r ξ W S . In this experiment, the existence of boundaries makes it challenging to obtain the localization length. To overcome this difficulty, we propose another method to extract the WSL length. With the maximum photon occupancy probabilities at Q 5 , defined as P max 5 ∶= max t>0 P 5 (t), we can obtain the WSL length using ξ W S ∝ 1 ln P max 5 . In Supplemental Note 2, we present a phenomenological proving of this relation. To extract more reliable P max 5 , we use Gaussian function to fit P 5 (t) and take the corresponding peak value as P max 5 , see Fig. 3(a). Now we study the relation between the potential gradient F and P max 5 . For a WSL system, the localization length ξ W S is inversely proportional to F , i.e., ξ W S ∝ 1 F . Hence, we expect that ln P max 5 ∝ F . According to Fig. 3(b), we can find that both the numerical simulation and experimental There is no crossing between the curves ⟨ρ K 1 (t)⟩ and ⟨ρ K 4 (t)⟩. The existence of this energy gradient between two edges shows that the energy transport is compressed. That is, the energies at one side can hardly spread to the other side.
results are consistence with this relation.
Thermal transport. Now we focus on the thermal transport in this system. For a 1D chain, the energy density at the j-th bond is defined asρ E j =Ĥ j,j+1 ≡ ρ K j +ρ P j , whereρ K j andρ P j denote kinetic energy and potential energy densities, respectively. From Eq. (2), these two quantities can be expressed aŝ In general, the thermal transport is closely related to the electronic charge transport in a classical metal system, which is known as the Wiedemann-Franz [49,50] law, i.e., λ σ = LT , where λ is thermal conductance, σ is electronic conductance, T is temperature, and L is Lorenz number. Eq. (5) shows that the potential energy only depends on the spin distribution, which displays BOs and WSL as discussed in the previous section. Here, we consider the time evolution of the kinetic energy density. In Fig. 1(d), the pulse sequences of this experiment are presented . To study the transport ofρ K j , the kinetic energy densities should exist a gradient between two edges at the initial state. Here, we choose the initial state as X + X + 000⟩, where X + ⟩ = 1 √ 2 ( 0⟩ + 1⟩) is the eigenstate ofσ x with eigenvalue 1 and can be prepared by X 2 gate. We can verify that, with this initial state, the kinetic energy at left edge is larger than one at right edge, so this initial state can be used to study the thermal transport. Then, ⟨ρ K 1 (t)⟩ and ⟨ρ K 4 (t)⟩, i.e., the kinetic energy densities of two edges, are measured, where the simultaneous readout of the nearest-neighbor two qubits is necessary. As shown in Fig. 4(a), when F = 0, the kinetic energies of two edges can exchange almost freely. Nevertheless, from Fig. 4(b), we can find that the difference between ⟨ρ K 1 (t)⟩ and ⟨ρ K 4 (t)⟩ always exist, when F 2π = 15 MHz. Therefore, similar to the spins, the thermal transport is also suppressed under the linear potential.
Due to the U (1) symmetry, the quench dynamics can be decomposed into different particle-number subspace, and different subspaces are decoupled with each other. For the initial state X + X + 000⟩, the photons only bunch at Q 1 or Q 2 . Despite existence of two-excitation populating for this initial state, Hamiltonian (2) can still effectively describe the dynamics of this system, since two excitations can hardly bunch at a same site due to large and staggered U j . Thus, we can use Slater determinant to calculate the dynamics of two-excitation sector. We can verify that the spins are localized among all of these subspaces with F ≠ 0. The spins can hardly propagate to the other side, so Q 4 and Q 5 almost remain at the initial state 00⟩. Therefore, the change of ⟨ρ K 4 ⟩ is small in this case, i.e., the kinetic energy can hardly transport from the left edge to the right edge. In this picture, we can know that the restriction of energy transport origins from the localization of the spins, which is identified with the classical Wiedemann-Franz law.

Discussion
In summary, we have reported the experimental observation of BOs and WSL on a 5-qubit superconducting processor. We provide another representation of the WSL length for a finite size system, i.e., the probability that a photon can propagate from one edge to anther edge. Using this representation, we verify that the WSL length is inversely proportional to the potential gradients. Furthermore, benefiting from the precise simultaneous readout of two qubits, the thermal transport in this system is also studied. The evolution of the energy densities shows that the thermal transport, akin to the spins, is not free under the linear potential, neither.
Comparing to the other artificial quantum manybody systems, one of the most significant advantages of the superconducting quantum circuits is that the states of superconducting qubits can be measured in an arbitrary basis. Thus, it enables us to study the thermal transport associated with BOs, which is generally a challenge for other platforms. Our results reveal that the superconducting quantum circuits can be considered as alternative synthetic quantum systems for experimentally exploring BOs and other quantum physics. Our platform may be useful for the further study of BOs, such as studying the BO frequency and spin current (see Supplemental Note 3), and imaging the Bloch band through BOs [16]. Our platform can also be extended to studying the transport phenomenon in other specific systems, for instance, in the presence of disorder potentials or engineered noises. In addition, it is meaningful to extend this system to the interacting case, and the Stark many-body localization may be realized in this system [51,52]. To explore these problems, our system could be scaled to include more qubits with longer decoherence time.

Methods
Setup. This 5-qubit device is made in the following processes: i) Depositing Aluminum.
A 100-nmthick Al layer is deposited on a 10 × 10 mm cplane sapphire substrate by means of electron-beam evaporation with a base pressure lower than 10 −9 Torr. ii) Etching the wires, resonators, and capacitor. We use a direct laser writer (DWL66+) and wet etching to produce microwave coplanar waveguide resonators, transmission lines, control lines, and capacitors of the Xmon qubit. The resist used here is S1813, and wetetching process is carried out with Aluminum Etchant Type A. iii) Fabricating Josephson junctions. The Josephson junctions of qubits are fabricated by the double-angle evaporation process. In this step, the undercut structure is made by a PMMA-MMA double layer EBL resist following the process similar to one reported in Ref. [45]. During the evaporation, the bottom electrode is about 30 nm thick, while the top electrode is about 100 nm thick with intermediate oxidation.
We package the device in an aluminum alloy sample box and fix the box on the mixing chamber stage of a dilution refrigerator. The temperature of the mixing chamber is below 15 mK during measurements. In order to reduce the external electromagnetic interference, an aluminum can and a µ-metal can are placed outside the sample box.
For each qubit, microwave pulses are applied through XY lines to rotate the qubit state between 0⟩ and 1⟩. Such XY pulses are formed by modulating continuous microwave signals sent from arbitrary waveform generators (AWGs: Zurich instruments HDAWG) via IQ mixers. To control all 5 qubits, the signal from a microwave source is divided into 5 channels through a power splitter, and each channel is amplified by a 11 dBm level. Current pulses are applied through Z control lines to tune the qubit frequencies. We use a DC current source (Yokogawa GS220) to apply static direct current to bias a qubit to its idle frequency and use an AWG to apply a fast current pulse to tune the qubit frequencies dynamically. Such static direct current and fast current pulses are combined by a bias-Tee, of which the capacitor is removed.
Readout pulses are composed of five tones at 40-MHz intervals. Each pulse corresponds to one qubit and is applied through the readout line. The output signals are amplified by a broad band Josephson parametric amplifier (JPA) [1] and a low temperature HEMT amplifier before further enhancement by a room temperature amplifier. The amplified signal is demodulated by a IQ Mixer and acquired by an analogdigital converter (ADC: Alazar ATS9360).
Attenuators, filters and isolators are used to reduce and isolate the noise from the electronic instruments, active electronic components (such as JPA and HEMT) and passive components outside the mixing chamber.
Error estimation. In our experiments, for the singlequbit readout, e.g., Fig. 3(a), each point shows the average of 6 × 100 single-shot measurements. To estimate the errors, we equally divide these single-shot readout data into 6 groups (each group contains 100 readouts). Thus, we can obtain 6 expectation values for each point, and the error bar is the standard deviation of these 6 expectation values. For the two-qubit readout, e.g., Fig. 4, each point shows the average of 10 × 200 single-shot measurements. We use the same method to estimate the errors, where the readout data are equally divided into 10 groups.
Numerical methods. The numerical results are obtained by numerically solving the Lindblad master equation, which reads whereρ(t) is the density matrix at time t, and Lindblad operatorsΓ n = 1 T 1ân andÂ n = 1 T * 2â † nân represent the excitation leakage and dephasing, respectively. The corresponding parameters applied here have been calibrated experimentally, and the details are shown in Supplemental Note 1. A.

Data availability
All data not included in the paper are available upon reasonable request from the corresponding authors.

SUPPLEMENTAL MATERIAL OBSERVATION OF BLOCH OSCILLATIONS AND WANNIER-STARK LOCALIZATION ON A SUPERCONDUCTING PROCESSOR
This Supplementary Information contains details of the experiment including: the experimental setup, amplification performance of the Josephson parametric amplifier (JPA), characterization of frequency multiplexed readout, decoherence time of each qubit, coupling strength between nearest neighbor qubits, Z control line crosstalk calibrations, delay time calibration for all control channels, square pulse distortion corrections, preparation of initial states like X + X + 000⟩, and calibration of dynamical phase induced dune to frequency tuning. Finally, we present a phenomenological analysis about the Wannier-Stark localization length in a finite-size system, and give a numerical result of the dynamics of spin currents.

D. Z control line crosstalk
Each qubit is tuned by current pulses applied to its Z control line. However, there are crosstalks between Z control lines, which means that the current pulse on one Z control line can also affect other qubits. To evaluate the influence of such crosstalk, we bias target qubit Q i at f i which is more than 1 GHz below its maximum frequency, so that the qubit frequency is sensitive to external magnetic flux as well as the Z pulse crosstalk. We apply a long (2µs) Gaussian shaped resonant X (π) pulse, and then measure the excitation probability of Q i . Applying a Z pulse with amplitude Z j to the Z control line of another qubit Q j will generate a crosstalk pulse on Q i , as shown in Supplementary Figure S6. Such a crosstalk pulse would lead to a frequency change for Q i , and hence the decrease of the excitation probability. We apply a compensation pulse with amplitude Z comp 1 to the Z control line of Q i to cancel this crosstalk. In practice, we monitor the excitation probability as a function of compensation pulse amplitude. When a complete compensation is achieved, the excitation probability is maximized. The linear relation between Z j and Z comp 1 is shown in Supplementary Figure S7.
After completing such crosstalk-compensation measurement for all qubits, a full transformation matrix for Z control line crosstalk calibration is constructed and shown in Eq.S1, where Z i represents the actual amplitude applied on the Z control line of Q i , and Z * i represents the effective amplitude after taking crosstalk compensation into account.
FIG. S6. Schematic diagram for Z control line crosstalk compensation. Qi is bias to a magnetic flux sensitive point with a frequency fi. When a square pulse with amplitude of Zj is applied to another qubit Qj, it could generate a crosstalk pulse with amplitude Z ji crosstalk . To cancel out the influence caused by this crosstalk and to maintain the qubit frequency to fi, we apply Z comp i = −Z ji crosstalk amplitude to the Z control line for Qi.

E. Calibration of time delay for all control channels
To guarantee the pulses applied on different control channels reaching to chip at the expected time, the output delay time of AWG ports need to be adjusted to compensate the transmission time difference in different channels. In our experiment, we set Z control channel of Q 1 as the reference, and make all other control channels align to it. There are two steps: In the first step, we align Z control channels by aligning Q 2 to Q 1 , followed by Q 3 to Q 2 and so forth. In the second step, we align XY control channel of each qubit to its Z control channel.
Scheme for determining Z control channel delay length is shown in Supplementary Figure S8. The blue lines represent transmission time in control channels. The red lines are the output delays that we insert before AWG output. The black dots are the AWG output start times. By adjusting the length of output delay, we can make the Z pulses align to the reference channel. In experiment, we set Q i (or Q i+1 ) to 1⟩, then apply two square pulses to make Q i and Q i+1 resonantly coupled for a time duration 2π 4g. This means that if two pulses are aligned, the resonant coupling will induce a complete photon swap. Here g is the coupling strength between Q i and Q i+1 . During the measurement, we vary the length of output delay time of Q i+1 and measure Q i (or Q i+1 ) after the Z pulses. The final state of Q i (or Q i+1 ) depends on the overlap of two Z pulses. When two pulses are aligned, the overlap will be equal to 2π 4g, and the swap probability will be the highest. The experimental results and fittings are shown in Supplementary Figure S10. In our experimental setup, the cables connecting qubits and electronics are designed to be of the same length, so the final adjustment of the output delay is in the range of ±500ps.
After the alignment of all Z control channels, we align XY control channel of each qubit to its Z control channel in time. As shown in Supplementary Figure S9, the effective excitation probability of Q i is the integral of Gaussian shaped π pulse over time before the Z pulse. When we adjust the output delay of XY channel, and measure the final state of Q i , we obtain the data in Supplementary Figure S11. By fitting the data, we can get the proper output delay time of XY channel. With cables of the same length, the measured delay times are also small and in the range of ±500ps.
FIG. S9. Schematic diagram for aligning XY channel to Z channel. The same as Fig. S8, we use blue lines to represent the transmission time through control channels, and use red lines to represent the individually adjustable delay before AWG output ports. By adjusting the length of delay, we can align π pulse in XY channel to square pulse in Z channel.

F. Square pulse distortion correction
In our experiment, qubits are tuned to specified frequencies by square pulses. However, an ideal square pulse is usually distorted when reaching to the chip. To correct this distortion, we need to determine the deformed shape of the step edge of Z square pulse. The corresponding qubit is biased to a relative low frequency f z (more than 1GHz below its maximum frequency) to improve its sensitivity to the pulse shape deformation. As shown in Supplementary Figure S12, an amplitude fixed (Z amp ) step signal is applied to a qubit, and followed a π pulse with frequency f z and length 20ns. If there is no distortion, the qubit will be excited to 1⟩ by the π pulse. If there is small distortion after the step signal, qubit will not be completely excited to 1⟩. By varying the compensation amplitude ∆Z, we can find the value of full compensation when the maximum excitation probability is achieved. By changing the delay time t of π pulse, we can get a view of the distorted step signal response, as shown in the upper row of Supplementary Figure S13. With this response data, we can calculate the needed adjustments for the square pulses. Then, we repeat the measurement in Supplementary Figure S12 with corrected pulses , and find the step signal response is flattened, see the lower row of Supplementary Figure S13.
FIG. S12. Schematic diagram for measuring the step signal response. Qubit is biased to a low frequency fz (usually 1 GHz lower than its maximum frequency) to gain better sensitivity to flux variation. Zamp is fixed. The length of Gaussian shaped π pulse is set to 20ns for acceptable time resolution. The frequency of π pulse is set to fz and fixed. In experiment, for different delay time t, ∆Z are varied to find the full compensation to the distortion. When qubit excitation probability reaches the maximum, full compensation is achieved.

G. Calibration of initial state phase
In order to measure the energy density transport, we need to prepare 5 qubits in an initial state such as X + X + 000⟩, and measure the expectation value of correlated operator σ i x σ i+1 x (σ i y σ i+1 y ). In the experiment, we firstly rotate Q 1 and Q 2 to X-Y plane with a π 2 pulse and set the initial state of Q 1 as X + ⟩, then we need to adjust the state of Q 2 by changing the phase φ of π 2 pulse to make Q 2 at X + ⟩, see Supplementary Figure S14. Here, we use the Bloch sphere representation for describing the state of a qubit. Then these two qubits are resonantly coupled to each other at the frequency f c . (In our experiment, we need to resonantly couple all 5 qubits at this frequency, and set a linear frequency gradient around f c .) After time τ = 2π 4 ⋅ 2g, we measure the final state probability as a function of the phase φ. The corresponding result is shown in the leftmost figure in Supplementary Figure S15. The curves of P 01⟩ and P 10⟩ are consistent with theoretical and numerical analysis. The cross point, indicated by the black arrow, corresponds to the φ when the state of Q 2 is also X + ⟩. In a similar way, we can calibrate the phase between Q 3 and Q 2 , and so forth.
FIG. S14. Schematic diagram for calibrating initial state phases of qubits. Rotate both qubits to X-Y plane, and set Qi initial direction as the reference X+. Adjust the phase (direction) of Qi+1, and make them resonantly couple for a time duration τ = 2π 4 ⋅ 2g, where g is the coupling strength between two qubits. At the end, both qubits are measured.

H. Calibration of dynamical phase
To measure the expectation value of correlated operator σ i x σ i+1 x (σ i y σ i+1 y ), we need to determine the phase of final state for each qubit. In our experiment, qubits are coupled at frequency f c , but are rotated and measured at their individual idle frequencies f idle . The frequency difference between f c and f idle can lead to a dynamical phase accumulation to the final state. Thus, in our experiment, we need to calibrate this accumulated dynamical phase.
The accumulation of the dynamical phase can be divided into three parts: rising edge part, middle flat part, and falling edge part. The summation, represented by the shadowed part in Supplementary Figure S16, can be expressed as Due to the presence of the rising and falling edges, the actual dynamical phase accumulation, comparing with ∆ω ⋅ τ , has an approximately fixed offset φ d = φ(τ ) − ∆ω ⋅ τ . To determine φ(τ ), we use the Ramsey-like method as shown in Supplementary Figure S16. The phase of the first π 2 pulse is φ 0 and fixed, and the phase of the second π 2 pulse is φ 0 + ∆ω ⋅ τ + ϕ. Thus, when ϕ = φ d , the phase of the second π 2 pulse is equivalent to φ 0 + φ(τ ), and the qubit will be in 1⟩.
The calibrated results of Q 1 , Q 2 , Q 4 , and Q 5 are shown in Supplementary Figure S17.

II. Supplementary Note 2. Wannier-Stark Localization Length
Here, by Jordan-Wigger transformation, the Hamiltonian (2) can be mapped to a free-fermion lattice model . Here we consider gij = g = 1, and the initial state 10000⟩.
whereĉ † j (ĉ j ) is the fermion creation (annihilation) operator. Here, without loss of generality, we consider g ij is siteindependent, i.e., g ij = g. Under a linear potential, i.e., h j = F j, the system can exhibit Wannier-Stark localization with localization length ξ W S = 2g F [S2]. This localization length can be understood classically as following [S2]: For a single particle, it can have the maximum kinetic energy E K max = max 2g cos(k) = 2g. In addition, the particle need to consume the kinetic energy h when hopping to the next site due to the linear potential. Thus, the maximum distance that a particle can propagate is about 2g h, i.e., Wannier-Stark localization length.
For a finite-size system, similarly, due to the linear potential, the particle cannot propagate completely from one boundary to another. In this picture, we can find that the maximum probability of a particle propagating from one boundary to another (i.e., P max 5 in the main text.) can indeed reflect the localization length. We assume that the localized wavefunction has the form where A is a normalized factor, n represents that the wave function is localized at site n, vac⟩ is the vacuum state ofĉ j . Therefore, intuitively, we expect that the maximum probability of a particle propagating from one boundary to another is proportional to e −α ξ W S , where α is the corresponding factor. In Supplementary Figure S18, we present the numerical result of the relation between g F and P max

5
, and we can find that P max 5 ∝ g F ∝ 1 ξ W S with F g > 0.6.

III. Supplementary Note 3. Spin Current
In this section, we consider the spin current under the linear potential. The spin current of Hamiltonian (2) can be defined asĴ where n represents the n-th bond of the chain. It can also be measured in our platform by the joint readout of two nearest-neighbor qubits. Here, we present the numerical results [see Supplementary Figure S19]. When F = 0, the spin current has nearly no decay as the increase of the propagation distance, while it decays quickly for F ≠ 0. Thus, this behavior of spin currents is an another signature of Wannier-Stark localization.