Thouless pumping in Josephson junction arrays

Recent advancements in fabrication techniques have enabled unprecedented clean interfaces and gate tunability in semiconductor-superconductor heterostructures. Inspired by these developments, we propose protocols to realize Thouless quantum pumping in electrically tunable Josephson junction arrays. We analyze, in particular, the implementation of the Rice-Mele and the Harper-Hofstadter pumping schemes, whose realization would validate these systems as flexible platforms for quantum simulations. We investigate numerically the long-time behavior of chains of controllable superconducting islands in the Coulomb-blockaded regime. Our findings provide new insights into the dynamics of periodically driven interacting systems and highlight the robustness of Thouless pumping with respect to boundary effects typical of superconducting circuits.


Introduction
Josephson junction arrays (JJAs) and their intricate many-body physics have captivated the attention of researchers since ground-breaking experiments in the 1990s (see the review [1]).Due to the possibility of engineering complex networks and their long-range coherence, JJAs have also been one of the leading candidates for quantum simulators in solid-state devices.Their practical application as quantum simulators, however, has been hindered by technological limitations: on one side, difficulties in tuning their physical parameters entailed the necessity of fabricating multiple devices to explore their phases of matter (see, for instance, Refs.[2][3][4][5][6]), thus impeding detailed investigations; on the other, irregularities in the selfcapacitance, induced charge, and Josephson coupling of the superconducting elements, resulted in an uncontrolled disorder.
These limitations have been mitigated by recent advancements in epitaxial growth techniques [7] that have paved the way for the realization of clean superconductor-semiconductor (SC-SM) interfaces.These breakthroughs enable an unprecedented tunability of the Josephson couplings through electrostatic gates [8][9][10][11][12], as well as the fabrication of multiple quantum dots on the same hybrid device [13,14], thereby revolutionizing the potential of JJAs as quantum simulators.Moreover, SC-SM platforms allow for on-chip patterning of arbitrary geometries in one and two dimensions, a precise control of the magnetic fluxes in these systems [11,12], and a relatively easy scalability.All these elements motivate the theoretical design of novel phases of matter and quantum simulation protocols on controllable JJAs.
In this respect, topological phases of matter are an ideal target for quantum simulations in solid-state platforms due to their intrinsic robustness against disorder, noise, interaction, and possibly dissipation.When combined with a time-periodic driving, one can engineer novel out-of-equilibrium states with no static analogs, known as Floquet topological phases [15][16][17][18][19][20].Among these, one of the simplest yet profoundly intriguing examples is Thouless pumping [21][22][23], a phenomenon arising in one-dimensional (1D) insulators with suitable timeperiodic driving of the system parameters.The topology behind this phenomenon leads to quantization of the charge adiabatically pumped during each driving period; in 1D JJAs, this corresponds to a current I = 2eΩC, with Ω being the pumping frequency and C a suitably defined Chern number characterizing the filled energy bands.Although theoretically well understood, experimental implementations of Thouless pumping have so far been confined to systems of ultracold atoms [24][25][26][27][28], optical waveguides [29], magneto-mechanical systems [30], and superconducting quantum processors [31] which fall short in capturing genuine transport phenomena of charged particles.
In this paper, we propose an innovative approach that combines a JJA with the ability to finely tune the induced charge on each SC island and the corresponding Josephson couplings.Through numerical simulations of 1D arrays in the Coulomb-blockaded regime, we focus on the long-time behavior of such systems.To investigate this limit, we connect a Josephson junction chain with superconducting leads that act as Cooper pair (CP) reservoirs.We study the effect of the electrostatic repulsion arising from the cross-capacitance between neighboring islands and show that topological pumping is remarkably robust with respect to both the coupling to the leads and nearest-neighbor interactions.Our proposal is qualitatively different from former experiments based on geometrical pumping which rely on optimal control of the pumping protocol (e.g.Ref. [32] for electron junction systems and Ref. [33,34] for superconducting transistors).Topological pumping, on the other hand, is predicted to be robust against disorder [35][36][37] and imperfections in the modulations [22,38].Its resilience against disorder, in particular, has also been experimentally investigated in ultracold Yb gases [27].
Our findings shed new light on the role of interaction and dissipation in topological pumping schemes, which are currently at the core of intense debate [23,[39][40][41].Through this research, we expand our understanding of JJAs as versatile platforms for quantum simulations while unveiling new insights into topological phenomena and the dynamics of interacting systems.
The paper is organized as follows: In Sec. 2 we derive a bosonic model from the Hamiltonian of a 1D JJA.In particular, we specialize on two models characterized by nontrivial topological properties: the Rice-Mele (RM) model [42,43], which we introduce in Sec.2.1, and the Harper-Hofstadter (HH) model [44,45], presented in Sec.2.2.Our analysis of topological pumping in Josephson junction chains based on the RM model is presented in Sec. 3, where we investigate in detail the role of the coupling with external superconducting leads and the effects of nearest-neighbor interactions.In order to show that the findings are not modelspecific, but hold on a wide class of periodically driven systems, we report further numerical analysis on the HH model in Sec. 4. In Sec. 5 we discuss the energy scales and constraints relevant to realistic experimental implementations.Finally, we draw our conclusions and present future outlooks in Sec. 6. Appendix A provides a brief overview of the relation between the quantization of the pumped charge and the Chern number associated to the pumping scheme.Appendix B presents data about the RM protocol for weak charging energies.

Josephson junction arrays in hybrid superconductor -semiconductor platforms
Recent developments in the epitaxial growth techniques of hybrid SC-SM materials [7] allow for the fabrication of 1D and 2D arrays of superconducting islands lithographically patterned in arbitrary geometries [10][11][12] and contacted to a semiconducting substrate through atomically pristine interfaces [46].In such devices, the filling of the substrate can be controlled via a global electrostatic gate (a top gate in the experiments in Refs.[10][11][12]).Additionally, smaller gates can be used to locally change the potential and density of states of the SM [9].
In the following, we consider devices defined by a 1D chain of superconducting islands of sub-micrometer size.For Al islands, the typical critical temperature of such systems is ∼ 1.6K [10], and we regard the related superconducting gap as the largest energy scale in the description of these systems.As a consequence, when operating at temperatures of the order of a few tens of milliKelvins, as customary in experiments, we can neglect effects determined by the Bogoliubov quasiparticle excitations of the islands.We can thus describe the low-energy physics of these JJAs by considering solely the dynamics of their CPs.
For a neighboring pair of SC islands, the dynamics can be modeled through the interplay of two kinds of interaction; first, the electrostatic interaction determined by the capacitance Figure 1: Chain of superconducting islands (dark blue) on top of a semiconducting substrate (green).Next to every island is a T-shaped side gate with tunable voltage V g, j .At each end of the chain, there is a superconducting lead.matrix which describes not only the potential induced on one island by the charge of the other, but also the charges induced in both islands by the surrounding environment.Second, the tunneling of CPs between the two islands, which define an effective SC-SM-SC junction where the coherent hopping of CPs is mediated by Andreev states [47].These states are induced in the SM layer below the two islands, where superconductivity is induced by proximity [9], and in the in-between region.Both interactions can be affected by neighboring electrostatic gates.Let us consider, for instance, a 1D chain as the one depicted in Fig. 1.In this setup, each superconducting island is addressed by a side gate at potential V g, j , controlling the induced charge 2en g, j = V g, j C g j .Additionally, we assume that a cutter gate at potential V c, j can control the filling of the semiconducting region in proximity of each Josephson junction and thereby its transparency.In this way, the effective coherent hopping amplitude of CPs between the islands can be modulated by varying the carrier density in the SM [9] and potentially turned off by totally depleting the substrate.
At low temperature, the JJAs can be modeled through a standard quantum phase model (see, for instance, the review [1]).To describe a chain as the one depicted in Fig. 1, we assign a superconducting phase operator φj and a number operator Nj to each superconducting island.Nj defines the number of CPs in the island with respect to an arbitrary offset.The two operators obey the standard commutation relation [ Nj , e i φj ] = − e i φj , and e i φj annihilates a CP in the island j.The Hamiltonian for the chain is written as Ĥ = ĤC + ĤJ and describes the interplay between the electrostatic interactions and CP tunneling, respectively.Let us first consider the charging energies of the islands in the chain, defined by Here C −1 is the inverse capacitance matrix of the islands and M is the number of islands in the chain.However, by assuming that the semiconducting environment and the electrostatic gates in the hybrid device effectively screen the charge of the SCs, we reduce the sum to single-island and nearest-neighbor interactions: In this equation, we introduced two customary energy scales: E C = 4e 2 /C self sets the charging energy of a single island, with C self being the sum of all capacitances to the other islands and environment elements; E C C = 4e 2 (C −1 ) j, j+1 determines the electrostatic energy between neighboring islands.We will assume for simplicity that these quantities are translationally invariant, which is the expected behavior for a strong screening imposed by the environment.Weak variations along the chain, however, do not affect our analysis and can be easily accounted for.n g, j is the charge induced in each superconducting island and can be controlled by the surrounding electrostatic gates (see Fig. 1).In particular, we assume that each induced charge n g, j is primarily controlled by the voltage of the side gate addressing the island j, n g, j = V g, j C g j /2e.Here C g j defines the mutual capacitance between the island j and its neighboring side gate.More complex scenarios to account for the charge induced by all electrostatic gates can be easily investigated.The energy scale E C is determined by the geometry of the islands and their electrostatic environment.We consider, as an example, gated devices with Al islands patterned over an InAs 2D electron gas; for rectangular islands of size ∼ 750nm×80nm, the resulting charging energy is approximately E C ≈ 0.125meV≈ h30GHz [48].
The coherent tunneling of CPs is modeled by the Hamiltonian Here E J, j is the Josephson coupling between island j and j + 1 and the Peierls phase x is the line integral of the vector potential ⃗ A along a path between island j and j + 1, which accounts for the role of the magnetic field when embedding the chain in a closed superconducting loop.Eq. ( 3) corresponds to tunneling of single CPs between neighboring islands.Hence, we are neglecting terms characterized by higher harmonics of the phase difference φj+1 − φj .In SC-SM-SC junctions, such an approximation is justified for cases in which the transmissibilities of the Andreev channels connecting the islands are low [47].This, in turn, corresponds to a sufficiently depleted semiconducting substrate.Higher-harmonic terms can be considered as well, giving rise to non-quadratic interactions in terms of the annihilation and creation operators of CPs (for instance, coherent tunnelings of charge 4e objects).These are, however, considerably weaker than the energy scale E J, j of the single-CP tunneling, and we expect them not to play a crucial role in the implementation of Thouless pumping.
Concerning the Josephson amplitudes E J, j , these can be globally controlled by a top gate as in Ref. [10], and the maximal value vastly depends on the width and length of the junction.Additionally, for the implementation of the RM model, we require the ability to separately address them through the voltage V c, j of suitable cutter gates.In general, the function E J, j (V c, j ) may be complicated (see the experimental data related to gate-tunable devices in Refs.[49][50][51]).We consider, however, a regime in which E J, j = 0 when V c, j is below a certain threshold V * c, j , such that the substrate is totally depleted, and, for simplicity, we impose that E J, j is approximately linear in V c, j above this threshold [see Fig. 3(a)].Importantly, the topological character of the RM model makes the pumped current robust against the details of the function E J, j (V c, j ) as long as it is monotonic and sufficiently regular.Therefore, we consider a device in which the average value of the Josephson coupling E J, j is set by a global gate, whereas an additional periodic time modulation can be imposed by the cutter gate voltages.Experiments on hybrid Josephson junction showed that, for a width of about 0.3µm, the amplitude E J at zero cutter voltage is of the order of h50GHz [9] and it can be switched off by applying sufficiently strong negative potentials (see, for instance, Ref. [50]).
When the Josephson energy dominates over the electrostatic terms, the system displays global phase coherence and behaves as a SC, allowing for coherent transport of CPs.Instead, the regime in which the electrostatic interaction E C dominates over the Josephson energies E J, j , results in an insulating phase unless all the induced charges n g, j are fine-tuned close to 1/2.In this scenario, the transport of CPs across the chain is suppressed, as we can consider each island in a Coulomb-blockaded state.In order to devise charge pumping protocols, we consider the regime E C, j > E J, j , E C C, j , and, initially, we neglect the nearest-neighbor electrostatic interactions.We can rewrite the Hamiltonian Ĥ by introducing the operators Σj = e i φj Figure 2: Charging energy of a single island as a function of the induced charge n g, j for states with N − 1, N , etc. CPs.The thicker red and blue curves mark the two charge states we consider in the hardcore-boson approximation and the shaded area is the maximum region of validity of such an approximation.The effective chemical potential is the energy difference between the N + 1 (blue) and N (red) parabolas.and Σ † j = e −i φj which, respectively, lower and raise the number of CPs in the island j by 1, such that [ Nj , Σj ] = − Σj .We obtain: This expression shows explicitly that the quantum phase model corresponds to a tight-binding Hamiltonian for the CPs.The dispersion of the charging energy of a single island is depicted in Fig. 2 as a function of the induced charge n g, j for different number states.When n g, j is a half-integer, states that differ by one CP are degenerate; let us consider, for instance, the case n g ≈ 0.5.If E J ≪ E C , we may assume that only the two lowest-energy states, with charge N = 0 and N = 1, are significantly occupied.States corresponding to the other parabolas are separated in energy by a gap ∼ 2E C and their population in the many-body ground state is negligible.Therefore, under the assumption E C > E J, j , T , where T is the system temperature, we can further simplify our description of the JJA and map the Hamiltonian Ĥ into a hardcore boson model.We may then replace Σj by a new hardcore boson operator bj , such that b † j bj = {0, 1} and b2 j = ( b † j ) 2 = 0.The energy difference between the two lowest charge states is E C (1 − 2n g, j ) for n g ≈ 0.5 (Fig. 2).Hence, we can define an on-site potential, resemblant of an effective chemical potential, which vanishes for n g, j = 0.5.We rewrite the total Hamiltonian (with E C C = 0) as a tightbinding model of hardcore bosons: This Hamiltonian shows that, in the hardcore limit E C ≫ E J, j , the system can be mapped into a chain of free fermions via a Jordan-Wigner transformation.Any time modulation of the electrostatic gates would translate into a time modulation of the onsite potentials µ j and the tunneling amplitudes E J, j .Hence, by implementing a suitable periodic modulation of these parameters, we expect that the Josephson junction chain is able to reproduce the physics of periodically driven systems of non-interacting fermions.
In the following, we will focus on the two most known schemes for the realization of adiabatic Thouless pumping: the periodically driven RM model and HH model.We restrict our analysis to the hardcore boson description: its breakdown and the effect of onsite interactions in the RM model have been theoretically investigated in Ref. [40].We provide a further analysis on the hardcore boson approximation specialized to JJA platforms in Appendix B.
Both models rely on a periodic drive of the onsite potential.We consider a generic time modulation of the side gates with V g, j (t) = V 0, j + δV g, j cos ωt + χ j , (7) where ω = 2πΩ.Such modulation yields: If µ 0, j ≈ μ is approximately independent of the position along the chain, it plays the role of an overall chemical potential for the hardcore bosons.The static voltages V 0, j can therefore be used to set the average filling of the system.In case of position-dependent variations of the C g j parameters, instead, these voltages can be tuned to reduce the fluctuations of the onsite potentials µ 0, j which, essentially, play the role of onsite disorder in the Hamiltonian (6).Finally, the oscillation amplitudes δV g, j determine the modulation of the onsite potential, which is additionally characterized by a position-dependent phase χ j that we will suitably set to implement the RM and HH models, as described in the next subsections.Throughout this paper, we set ħ h = 1.

Rice-Mele Hamiltonian
The RM model [42] offers the most paradigmatic example of topological pumping in 1D systems.It is a model with a two-site unit cell, and its hardcore boson formulation is defined by a time-periodic Hamiltonian of the form Here we consider a chain with an even number of islands M and open boundary conditions.The instantaneous spectrum of ĤRM in the thermodynamic limit displays two bands with a linear level crossing at E J,1 = E J,2 and µ A = µ B .For µ A = −µ B the ground state is found at halffilling, and it is useful to define the two-component parameter vector The time-dependent single-particle gap is given by | ⃗ h(t)|, and a topological charge pumping is obtained when ⃗ h(t) winds around the gapless point ⃗ h = 0 during one period.
To implement Thouless pumping, we adopt a modulation of the kind in Eq. ( 8) for the onsite potentials and, in particular, we choose V 0, j such that all islands are tuned close to the charge degeneracy point µ 0, j = 0 between the lowest energy parabolas (Fig. 2).In order to enforce we set χ j = (−1) j+1 π/2 and choose modulation amplitudes δV g, j such that δµ is approximately constant along the chain.We observe that a residual μ may reduce the energy gap of the system, but in the limit of adiabatic pumping, it does not affect the pumped charge as long The blue curve displays the approximation considered to model the Josephson energies E J, j as a function of the cutter gate voltage V c, j (arbitrary units).The dashed vertical line indicates the threshold V * c below which E J, j is considered zero.Yellow and blue shaded areas schematically indicate the voltage range for the sinusoidal and clipped modulations, respectively.Panels (b) and (c) display the corresponding modulated E J, j (ωt); green and purple curves refer to odd and even junctions, respectively.as it remains sufficiently smaller than the single-particle gap at all times.The modulation of the Josephson energies E J, j (V c, j ) is achieved by time-periodic voltages in the cutter gates.In particular, we adopt the following signals: We assume that all junctions in the chain are characterized by the same parameters, such that this modulation approximately results in Josephson amplitudes of the kind where f l r is a linear rectifier function, with f l r (x) = x if x > 0 and f l r (x) = 0 otherwise, see Fig. 3(a).V c is used to control the offset J 0 and the modulation δJ is roughly proportional to δV c .If V c − δV c > V * c , the Josephson amplitudes E J, j are always positive and display a sinusoidal modulation [see the example in Fig. 3 c , E J, j (t) = 0 for a fraction of the time period and its modulation is clipped below zero [Fig.3(c)].
The time-periodic nature of the Hamiltonian H RM allows us to consider the time coordinate as a second dimension for the momentum, such that we can assign a Chern number C to each energy band [21].For the considered modulations, the Chern numbers are where the index n = 0, 1 labels the single-particle bands, resulting in an adiabatic Thouless pumping at half filling with average current I = 2eC 0 Ω.

Harper-Hofstadter Hamiltonian
The 2D Hofstadter model [44] provides the simplest description of charged particles moving on a square lattice and subject to a uniform out-of-plane magnetic field.The dynamics of the model is determined by the magnetic flux per plaquette and the corresponding Aharonov-Bohm phase Φ acquired by a particle that moves around it.Incommensurate values of the flux result in the celebrated Hofstadter butterfly fractal spectrum, while for commensurate values Φ = 2πp/q, the system can be modeled in the Landau gauge by introducing a q-site unit cell, and displays q energy bands with non-trivial Chern numbers.
The realization of this model in the Josephson junction chains requires exclusively the modulation (8) of the onsite potentials and is obtained by setting a position-dependent phase χ j = Φ j in Eq. (7).In particular, setting Φ = 2π/3 yields the simplest example of the gapped Hofstadter model.It corresponds to a 3-island periodicity of the phase χ j and, in the thermodynamic limit, it yields three bands of the instantaneous spectrum (see Fig. 4), characterized by Chern numbers C n = 1, −2, 1.The band gap is proportional to the minimum between E J /2 and δµ.
The realization of the corresponding pumping scheme does not require a modulation of the Josephson amplitudes.Furthermore, the uniform and constant potential μ can be changed by regulating the voltage offsets V 0, j and be used to set the overall filling of the CPs.The Thouless pumping is achieved when μ lies in either of the band gaps depicted in Fig. 4. In particular, the charge pumped adiabatically in the thermodynamic limit is proportional to the sum of the Chern numbers of the filled bands: this implies that by varying μ from the first to the second gap, we expect the pumped current to change sign from 2eΩ to −2eΩ.

Coupling with external superconducting leads
So far, we have discussed the modeling of isolated Josephson junction chains.To investigate their transport properties, however, it is necessary to embed these systems in a larger environment.Specifically, our aim is to model the driven transport of CPs across such JJAs.Therefore, the most convenient choice is to include two superconducting leads coupled to the first and last island of the chain (depicted in light blue in Fig. 1 and labeled by SC L and SC R ).To this purpose, we supplement the Hamiltonian Ĥ with a boundary term that describes a Josephson junction between the two extreme islands and the related superconducting leads [52] where E L describes the associated Josephson energy which we consider constant in time.By assuming that both leads are described by standard BCS coherent states, we can replace the related phase operators φL/R with classical phases ϕ L/R .When embedding the JJA in an external SQUID, the phase difference φ = ϕ R − ϕ L would correspond to the magnetic flux threading it.
We observe that such a phase difference can be absorbed in the Peierls phases θ j, j+1 in Eq. ( 3) by a suitable gauge transformation.
In the hardcore limit, the boundary Hamiltonian Ĥb becomes This boundary term violates the conservation of the total number of CPs in the chain, such that, in general, the many-body ground and Floquet states will be superpositions of different particle numbers.Hence, the adiabatic condition no longer involves the single-particle band gap e g but a many-body energy gap E g that depends on the interplay between the charging energies of the SC islands, the Josephson energies of the junctions, and the coupling with the leads.
We emphasize that the boundary Hamiltonian ( 15) provides a reliable description of the coupling with an external environment, despite the fact that it constitutes a coherent coupling.This is due to the fact that the external leads are treated as condensates of CPs, such that only the exchange of particles with zero energy is allowed with the environment.This avoids the non-coherent dissipation processes that would characterize, instead, a coupling with a fermionic or bosonic thermal bath (as in the case of normal metallic leads of phononic baths) and would require a thorough study of the related non-unitary dynamics.

Results of the Rice-Mele model
In the following, we will investigate Thouless pumping in short chains coupled to external SC leads to identify experimentally relevant regimes where charge quantization can be observed.In particular, we will focus on the role played by the Josephson coupling to the leads, which breaks particle number conservation, and on the effects of the nearest-neighbor interaction E C C .To this end, we adopt a Floquet description of the driven Josephson chain in the hardcore boson limit.We consider the Hamiltonian ĤRM,tot (t, φ) = ĤRM (t) + Ĥb (φ) , (16) which describes a RM chain, as defined in Eq. ( 9), embedded in a superconducting circuit through the contact with the leads in Eq. (15).ĤRM,tot (t, φ) is periodic in both t and φ, and we can analyze the dynamics of the related system by defining a many-body Floquet operator which, for a generic Hamiltonian Ĥ, reads Here T indicates time ordering and τ = 1/Ω is the time period such that Ĥ(t + τ) = Ĥ(t).We consider system sizes up to M = 8, compatible with realistic devices in which all islands can be independently addressed by gate voltages, and we numerically diagonalize U(τ) to obtain the many-body Floquet eigenstates {|Φ ν (τ)〉}: where E ν labels the many-body Floquet quasienergies.In the infinite-time limit, for a given phase difference φ, the time-averaged pumped charge per cycle corresponds to [22,53,54] The table represents the deviations from perfect quantization of the phase averages pumped charge Q as a function of the systems size M . with Here Ψ 0 is the state of the Josephson chain at the beginning of the pumping protocol, which we set to the ground state of ĤRM,tot (t = 0, φ), ρ(t) is the density matrix of the system at time t, and 2e∂ φ Ĥ = Ĵ corresponds to the current operator.In the adiabatic limit, the ground state |Ψ 0 〉 of ĤRM,tot (t = 0, φ) has a large overlap with the Floquet state with the lowest energy expectation value [45], which we denote by |Φ 0 (τ)〉, such that N ν → δ ν,0 .Hence, Q ∞ is carried by a single Floquet state.An equivalent expression for the pumped charge per cycle in the infinite-time limit is [54] where Ĵ (t) is the time-dependent current operator, and |Φ ν (t)〉 are the time-evolved Floquet states within one period.Given the difficulty of accurately differentiating numerically the quasienergies with respect to the external SC phase, due to the many (avoided) level crossings in the Floquet spectrum, we will use Eq. ( 21) to extract the pumped charge from the time evolution.Our numerical calculations are based on exact diagonalization and performed with the QuTiP Python framework [55,56].The phase difference between the two external superconducting leads, rescaled by the number of islands in the chain, φ/M , can be interpreted as a shift of the quasimomenta of the system CPs via a suitable gauge transformation.For small systems, this shift causes a nontrivial dependence of ∂ φ E ν on φ that yields, in turn, a dependence of Q ∞ on the external phase.This is expected to be relevant only for small systems and vanishes completely in the thermodynamic limit, provided the system is in an insulating phase [57].To evaluate these finite-size effects for the Hamiltonian ĤRM,tot , we depict in Fig. 5 the time-averaged pumped charge Q ∞ (φ) for E C C = 0, E L = 2J 0 and a ratio ∼ 1/16 between ω and the many-body gap E g .For M = 4 and M = 6 islands, the dependence on φ is significant, but the oscillation around the average quantized value reduces to ∼ 2% for M = 8.We find analogous results for the HH model, see Sec. 4. To recover an exact quantization in short chains, it is therefore necessary to average the pumped charge over φ.By following the Thouless construction [35], we define For systems with periodic boundary conditions conserving the particle number (such that φ is a phase twist in the boundary conditions), this quantity is quantized in the adiabatic limit at half-filling and corresponds, when N ν = δ ν,0 , to the Chern number of the lowest energy band (see Appendix A).More in general, Eq. ( 22) implies that only Floquet many-body states with a quasienergy that winds in the Floquet-Brillouin zone as a function of φ contribute to the pumped charge Q (see, for instance, Fig. 7).Our results show that Q displays a precise quantization towards the adiabatic limit within an error of 1%, even when including the boundary term Ĥb .The residual small deviation originates from non-adiabatic effects due to the finite frequency and the choice of the driving protocol [43,45].Adiabatic perturbation theory gives a leading-order estimate of this correction which scales as 1 In the following, we will consider small system sizes, for which the phase dependence is sizable, and we will investigate the dependence of the pumped charge on the modulation δJ of the Josephson coupling, the coupling E L , and the nearest-neighbor interaction E C C .

Role of the Josephson modulation
The amplitude δV c of the cutter gate modulation and, consequently, the amplitude δJ in Eq. (12), determine the waveform of E J, j (t) (see the examples in Fig. 3).In this respect, the role of δJ is twofold: on one hand, it determines the minimum single-particle bulk gap ε g = min | ⃗ h(t)| over one period, which is given by min δJ, To estimate these effects, we consider the dependence of Q ∞ on δJ for different values of φ (Fig. 6) in an array with M = 6 islands.In the sinusoidal regime, where δJ is small, finite-size effects induce a large dispersion of Q ∞ as a function of φ, and averaging over the phase is fundamental to obtain the charge quantization.The dispersion rapidly decreases by increasing δJ, until the variation of Q ∞ (φ) drops below 0.8% when δJ ≥ J 0 .The weak dependence of Q (red crosses) on δJ, instead, originates from non-adiabatic corrections.These are suppressed with δJ since the band gap ε g increases accordingly, bringing the system deeper into the adiabatic regime without changing the driving frequency.In particular, the residual nonadiabatic and finite-size corrections of Q are less than 0.6% for ω = 0.05J 0 in the clipped regime.
This dependence on the Josephson energy modulation can be understood by comparing the behavior of populated Floquet many-body states as a function of φ for δJ = 0.5J 0 (sinusoidal waveform) and δJ = 1.2J 0 (clipped waveform).They are shown in the left and right panels of Fig. 7, respectively, where we plot the Floquet quasienergies as a function of the phase φ, using the corresponding occupation number N ν to set the color of each data point.In both cases, a single low-energy Floquet state |Φ 0 (τ)〉 has an almost perfect overlap with the initial ground state, allowing us to follow the winding of E 0 around the Floquet-Brillouin zone and confirming the quantization of Q, which follows from Eq. ( 19).In the sinusoidal regime, E 0 is, in general, far from linear in φ, leading to a strong dependence of Q ∞ (φ) on this phase.In the clipped regime, instead, the derivative ∂ φ E 0 of the most populated state is practically constant, and, for the chosen value of ω, ) for all values of φ.We conclude that even with small systems (M = 6), the phase dependence is very weak in the clipped regime, and a good quantization can be observed without averaging over the SC phase difference.

Coupling with the leads
Besides the dependence of the pumped charge on the external phase φ, it is important to investigate the role of the coupling E L between the edge islands and the superconducting leads.When E L dominates over the other energy scales, the first and last islands are effectively merged with the two external leads and, in practice, the system behaves as a shorter chain.If instead E L is weak compared to J 0 , the transport of CPs across the system is hindered, as the pumped current is faster than the rate of charge transfer to or from the leads.Consequently, the system behaves as an open chain.To examine the interpolation between these two limits, it is instructive to consider the behavior of the system for a clipped modulation δJ > J 0 .In this situation, the onsite potentials µ A/B and E J,1 vanish in the middle of the pumping period, t = (n + 1/2)τ.Therefore, the ground state of the Hamiltonian ĤRM (τ/2) corresponds to the extreme topological dimerized state of the Su-Schrieffer-Heeger (SSH) model, and it displays the typical four-fold degeneracy associated with localized zero-energy boundary modes.The boundary Hamiltonian Ĥb , however, gaps the SSH edges; the ground states of the first and last islands in this dimerized limit result and are separated by an energy gap E L from the localized excited states in the same islands.Therefore, when E L < J 0 +δJ 2 , 2δµ, it sets the many-body energy gap E g .Notably, the limit E L → 0 corresponds to the isolated RM chain in which the zero-energy SSH edge modes cause transitions between the two energy bands, thereby disrupting the Thouless pumping.Only when E L is sufficiently larger than ω, this non-adiabatic effect vanishes and the quantized pumping can be restored.
Fig. 8(a) displays the pumped charge Q ∞ (φ = 0) and the many-body gap E g as a function of E L /J 0 for a modulation δJ = J 0 and a driving frequency ω = 0.05J 0 .The two quantities are clearly correlated; Q ∞ saturates to the quantized value as the gap induced by the coupling with the leads becomes sufficiently large.As expected, E g eventually saturates as well at the value E g = J 0 , when E L ≥ J 0 .In the case of sinusoidal modulation, instead, the dimerized SSH limit is never reached during the driving period.Therefore the many-body gap saturates at smaller values for large E L /J 0 .
The region 0 < E L < J 0 is the most susceptible to non-adiabatic errors since the ratio between the many-body gap and the driving frequency changes continuously between 0 and J 0 /ω.Hence, we expect the pumped charge to saturate faster to Q ∞ = 2e for smaller values of ω.In the adiabatic limit ω → 0, a finite gap opens for any value of E L > 0, leading to the quantization of Q ∞ .However, in a realistic scenario, the driving frequency is finite and sets the average magnitude of the current flowing through the array.Thus, E L needs to be large enough to prevent the CPs from accumulating on one side of the chain, similar to what happens  in an open system, suppressing the charge pumping.In the Floquet framework, this can be understood by observing the winding of the quasienergies associated with Floquet states with the highest occupation.When the charge is quantized and the driving is sufficiently slow, there is a clear correspondence between Floquet states and energy eigenstates, leading to an almost perfect fidelity N 0 ≃ 1, for any value of the phase difference φ.Instead, when E L < J 0 , there are multiple Floquet states that display a large overlap with the initial ground state, as shown in Fig. 8(b), which might interfere destructively.Moreover, large avoided crossings appear in the spectrum as φ changes, further suppressing quantized transport.We conclude that, in order to avoid additional non-adiabatic corrections to the pumped charge, E L must be larger than the single-particle band gap of the chain.

Nearest-neighbor interactions
In the most common trapped ultracold atom quantum simulators [23], nearest-neighbor interactions are negligible.Only very recently, motivated by experiments with dipolar atoms with long-range interactions [58], the role of nearest-neighbor repulsion in topological pumping schemes has been investigated in extended spinful RM models [59].
As in the case of dipolar gases, JJAs can be characterized by sizeable nearest-neighbor interactions.These interactions are bounded by E C C < 2E C , however, in the hardcore limit, they can be highly relevant because nothing prevents E C C from being of the same order of J 0 or the single-particle gap.For hardcore bosons, the nearest-neighbor interaction reads: Ĥint has two effects on the chain: a static repulsion, E C C b † j+1 b † j bj+1 bj , and a shift in the chemical potential, −E C C (n g, j−1 + n g, j+1 ), due to the charge of the neighboring islands.The staggered modulation (10) of n g, j in the RM scheme implies that the interaction (24) always favors, to a certain extent, the dimerization characterizing the half-filled topological state and the onset of charge-density states with occupations |0101 . ..〉 or |1010 . ..〉 that evolve one into the other.
In Fig. 9 we plot the pumped charge Q ∞ (φ = 0) at the onset of the clipped regime (δJ = J 0 ), as a function of both E L and E C C .Q ∞ is remarkably stable against the nearest-neighbor interactions and there is a wide region of the phase diagram where the pumped charge is quantized.The role of E C C seems indeed to be negligible for the Rice-Mele pumping: although the Hamiltonian H RM ,tot does not conserve the CP number, its many-body ground states are in good approximation the half-filled charge-density states where CPs are localized on every other SC island for most of the driving protocol.Therefore, the role of moderate interactions does not affect the pumping efficiency.In particular, for realistic ranges of E C C , our simulations indicate that the nearest-neighbor interactions do not cause any topological phase transition in the effective 2D Floquet topological insulating state that determines the onset of the quantization of the pumped charge.

Results of the Harper-Hofstadter model
The HH model contains two main differences from the RM model.First, its implementation does not require modulation of the tunneling amplitudes but only of the induced charges.This makes it more suitable for experimental implementations but, at the same time, more susceptible to finite-size effects since this protocol does not rely on any dimerization.Second, when the Aharonov-Bohm phase Φ is set to 2π/3, the model displays two insulating phases with opposite Chern numbers and filling 1/3 and 2/3. 1 This implies that by varying the parameter μ in Eq. ( 13), it is possible to invert the direction of the pumped charge.The two insulating states correspond to intervals in μ with bounds approximately given by the band gaps.These states are separated by a gapless phase in which pumping is not quantized.In addition, since the topological insulating Floquet states of the HH model appear at fillings 1/3 and 2/3, the role of nearest-neighbor interaction is less trivial than in the RM model: on one side, the HH pumping does not rely on the staggered charge configurations favored, for a broad range of parameters, by the interaction (24); therefore, for large values of E C C , we expect the HH pumping to be suppressed.On the other, with respect to the RM scheme, the different position dependence of the induced charges in the HH model reduces the interaction energy difference between the half-filled charge-density waves and density-wave states at filling 1/3 (for instance |100100〉) and 2/3 (for instance |011011〉), whose time evolution is beneficial for the HH pumping schemes and matches the E C C → 0 limit of this protocol.
In the hardcore limit, the two topological insulating phases of the HH model are related by a particle-hole-like symmetry which is fulfilled also in the presence of the nearest-neighbor interaction in Eq. ( 24).This symmetry is defined by br The second equation corresponds to the mapping µ 0, j → −µ 0, j and t → t +τ/2.This symmetry relates states with complementary fillings and implies that the pumped charge changes sign when reflecting n g, j around 1 2 .In particular, this holds for any value E C C , implying that the effect of the interaction is the same in both insulating phases.To control the filling factor and, therefore, the direction of the current, we tune the offset of the chemical potential μ through the average value of the induced charge during one driving period: Here we consider, for simplicity, uniform C g and V 0 across the chain.To investigate the stability of Thouless pumping in the presence of interactions, we consider the dynamics dictated by the Hamiltonian defined by the Hamiltonians in Eqs. ( 13), ( 15) and ( 24) with the drive obtained by a modulation of the induced charge with position-dependent phases χ j = 2π j/3.We calculate the charge using Eq. ( 21): the system is initialized in the ground state of the Hamiltonian in Eq. ( 28) at t = 0 and the time dependence of the Floquet states is computed by solving numerically the Schrödinger equation.We first examine the pumping for values of μ where the system lies deep in one of the Floquet topological insulator phases.Fig. 10 shows the behavior of Q as a function of the lead coupling E L for different values of the nearest-neighbor interactions, up to half the charging energy, E C C < E C /2 = 2E J .The behavior is qualitatively analogous to that observed for the RM pumping with clipped modulation (compare Fig. 8(a) and Fig. 10): a good quantization of Q is achieved for sufficiently strong E L ≳ 0.4E J .Importantly, even interactions comparable with the many-body gap (E C C = 1.9EJ in the data) do not significantly alter the behavior of the pumped charge, which appears remarkably stable.
We emphasize that for the HH model, averaging over the phase difference φ is necessary to obtain the predicted quantization of the pumped charge since we consider small system sizes (M = 6, 9).This is caused by the winding of the eigenenergies of the populated Floquet states, displayed in Fig. 11.For the most occupied state at filling 1/3 [panel (a)] the quasienergy winds three times in the positive direction and twice in the negative direction, resulting in a total winding number +1.Panel (b) in Fig. 11 shows instead the opposite winding for the gapped phase at filling 2/3, corresponding to pumping in the opposite direction.In both cases, the current ∂ φ E(φ) strongly depends on φ and even changes sign (see Fig. 12(a) for the state at filling 1/3).Therefore, a quantized pumping can be obtained only by averaging over φ.This is a characteristic of small-size systems.By comparing the results for M = 6 and M = 9 in Fig. 12(a) we see that the phase dependence considerably decreases with the chain length.However, in general, the HH model is considerably more influenced by finite-size effects than the RM model in the clipped regime (see Fig. 5

for comparison).
In an experimental setup, control over the voltage offset V 0 of the side gates provides a possibility of interpolating between the two gapped phases and exploring the phase diagram of the system as a function of μ [Fig.12(b)]: our results clearly show the appearance of the two topological Floquet phases with quantized charge Q = ±1 • 2e both without interactions, E C C = 0, and for strong interactions, E C C = E J /2. Moderate values of E C C seem indeed to be beneficial for the stability of the topological phases, as the interactions increase the width of the plateaus where the pumped charge is quantized.Above a certain threshold (E C C ≳ 2E J , not shown), however, the repulsions in Eq. ( 24) favor the charge-density waves states at filling 1/2, thus destroying the quantized HH pumping.
In the data of Fig. 12(b), the symmetry in Eqs. ( 25) and ( 26) is not exactly reflected.The discrepancy results from the fact that the initial time of the pumping protocol is the same for all values of μ.The missing translation t → t + τ/2 violates the requirement (26), but the effect is only seen when μ is close to the band extrema and the initial conditions strongly affect the pumping outcome.This strong dependence is shown explicitly in the inset of Fig. 12(b), where Q fluctuates with ng around the transition between the gapless metallic states (n g ≲ 0.35) and the Floquet topological insulator, (n g ≳ 0.35).Indeed, when μ lies close to a topological band edge, a small change in the SC phase difference φ can change the initial state between being either metallic or insulating.Consequently, Q ∞ (φ) has discontinuities in φ and its phase average Q strongly depends on the driving frequency, the precise sampling of φ, and the system size, which is reflected in the single-particle level spacing within the bands.
We show this irregular behavior in Fig. 12(a), where we plot the infinite-time-averaged pumped charge Q ∞ as a function of the SC phase φ.A comparison of the data obtained at the edge (n g = 0.35, orange triangles) and inside (n g = 0.4, teal circles) the topological phase clearly shows the different role of φ in the two cases.While the amplitude of the variation of Q ∞ is the same, in the topological phase the pumped charge has a smooth dependence on φ and it oscillates around the quantized value Q = 1 marked by the horizontal dashed line.On the edge between the metallic and the insulating phase, on the other hand, Q ∞ has many discontinuities that are reflected in the fluctuations of Q in the inset of panel (a).

Experimental implementation of the pumping schemes
Several aspects must be considered for designing an experimental realization of the proposed pumping protocols.In the following, we analyze first the tuning of the physical parameters to mitigate non-adiabatic and disorder effects.Then we discuss the main dissipative effects related to the coupling with the environment.

Physical parameters to achieve the topological pumping
The analysis presented in the previous sections and the observation of a quantized Thouless pumping rely on several constraints, important in devising an experimental realization of the presented models.
The most fundamental parameter in both driving protocols is the frequency Ω.On one side, Ω determines the ideal pumped current I = 2eΩC; therefore, in order to obtain clearly measurable currents (I ≳ 10pA), we require Ω ≳ 300MHz.On the other, the frequency determines the onset of nonadiabatic errors, typically scaling as (hΩ/E g ) 2 .Therefore Ω must be sufficiently small compared to the many-body gaps E g .For both the RM model with clipped modulation and the HH model at Φ = 2π/3, the gap is given by the minimum between 2δµ and a Josephson term proportional to E J (for sufficiently large E L ). δµ is proportional to E C , and to enter the regime of quantized pumping it is necessary that E C ≳ E J .E J can be tuned more easily than E C , for example by applying a suitable voltage to a global gate as in the devices analyzed in Ref. [10].Therefore, we consider the charging energy as the limiting factor in determining the gap which, for small superconducting islands, is approximately of the order of h30 GHz [48].This poses an upper limit to the frequency, Ω ≲ 10 GHz, to avoid excessive nonadiabatic errors.The frequency range Ω ∈ (300 MHz, 10 GHz) of the voltage drive of the electrostatic gates can be explored with standard waveform generators.
The RM model has the disadvantage of requiring modulation of both the induced charges and the Josephson couplings.This can be particularly challenging as the cutter gates controlling the Josephson energies may additionally induce a charge on the neighboring islands.Such an effect can, however, be compensated with suitable corrections of the voltages V g of the side gates.At the same time, the RM model displays, for M = 6 islands, smaller finite-size effects than the HH model, and, in particular, does not show a strong dependence of the pumped charge Q ∞ on the lead phase φ, especially in the clipped regime.For the same number of superconducting islands, the HH model requires control of only half of the electrostatic gates but displays instead a very strong dependence on the phase φ for M = 6 [Fig.12(a)].To avoid this limitation, two strategies can be envisioned: either implementing the pumping in longer chains or embedding the system in a device that allows for averaging in time over the phase φ, similar to the proposals of Refs.[52,60].
Concerning the scalability of the HH model at flux 2π/q, only q independent voltage signals are needed for sufficiently uniform SC chains, such that the capacitances C g j and charging energies present only minor variations along the system.These q signals need then to be distributed across suitable gate structures.
Regarding the averaging over the phase φ, several options can be envisioned.The first possibility is based on introducing a suitable voltage bias V b between the two external SC leads (see Ref. [52]).This yields a linear winding of the phase φ with period τ φ = h/2eV b .The number of pumped CPs in this period is given by Q(τ φ ) = CΩτ φ , and the average over the phase φ is suitably implemented if Ωτ φ ≫ 1.An alternative method would require instead embedding the HH Josephson junction chain in a superconducting ring, in order to impose a phase bias φ that can be varied in time through a driven magnetic flux.
In our analysis of short Josephson junction chains, we did not consider explicitly the role of disorder.For both the RM Hamiltonian ( 9) and the HH Hamiltonian (13), the onsite disorder corresponds to a position dependence of the time-independent part of µ j in Eq. ( 8).This can be caused by non-uniform capacitances C g j and a failure in balancing them with the voltages V 0, j .The effect of this kind of disorder on Thouless pumping has been extensively studied (see, for instance, Refs.[27,[35][36][37][38]); in general, the quantization of the pumped charge is robust as long as the random variations of the onsite potential are weaker than the energy gap.In Coulomb-blockaded JJAs, the disorder amplitude of the onsite energy depends on the variance of E C ng , whereas the energy gaps are determined by the Josephson energies.Therefore, to approach an accurate quantization of the pumped charge, we need to consider a balance between the following competing constraints.
On one side, the ratio E C /E J cannot be too large.Specifically, the standard deviation of the island-dependent E C ng (where ng represents the targeted n g average) must be considerably smaller than δJ and E J in the RM and HH pumping schemes, respectively.Our calculations rely on the CP hardcore assumption E C ≫ E J ; however, numerical investigations [40] of the RM model with interactions of the form (2) reveal that the RM Floquet topological insulator phase survives even when the ratio E C /δJ decreases to E C /δJ ∼ 3 if δµ is sufficiently strong.Therefore, keeping a moderate value of E C /E J ≳ 3 may be beneficial to reduce onsite disorder and nonadiabatic effects.On the other side, the Josephson energies cannot exceed the threshold corresponding to the insulator-SC transition in the static case.Thouless pumping can indeed be realized only when the instantaneous energy spectrum of the system and the related ground states at each time t correspond to insulating phases.In the appendix, we provide data acquired by relaxing the hardcore constraint.The short time dynamics of the RM model for δJ = J 0 confirms the results in Ref. [40] and suggests that quantized pumping is observed for charging energies above a threshold roughly given by E C > 2J 0 .By decreasing the ratio E C /E J , however, the nonadiabatic errors are expected to become more consistent, therefore lower pumping frequencies are required to observe a good quantization.
A different kind of disorder characterizing JJAs are random variations in the Josephson energies and, in the RM model, differences in the functions E J, j (V c, j ) associated with the junctions along the chain.Also in this case, the related disorder in the hopping terms of Eq. ( 6) becomes detrimental for an accurate quantization of Thouless pumping when the amplitude becomes comparable with the energy gaps of the systems.
Given the accuracy of the lithographic techniques adopted for the fabrication of hybrid JJAs, however, we expect the typical disorder amplitudes in the Josephson energies to be below 10% (see related experimental estimates in Ref. [61]).A disorder strength in this range is not harmful for the implementation of Thouless pumping as it is way below the necessary threshold to close the many-body gap and suppress quantized transport.

Coupling with the environment and dissipation
Our simulations considered exclusively the coherent and unitary time-evolution of the Josephson junction systems, such that they do not capture dissipative effects.The coupling with the environment and, in particular, to phononic baths can, in principle, affect the quantization of Thouless pumping.The most relevant dissipation effects will involve the population of states above the many-body gap E g .We observe, however, that, based on the above values for E C , we may estimate E g ≈ E C /3 ≲ K B 500mK.This gap must be compared with the typical temperature T ≲ 50mK of operation of dilution refrigerators.Therefore, in standard conditions, we expect the environment temperature to be considerably lower than the many-body gaps.In such a regime, it has been shown in the case of the non-interacting RM model that the quantization of the pumping is actually improved by the coupling with a cold environment with respect to the coherent time-evolution for a broad range of parameters [62].We do not expect a qualitatively different behavior in the presence of interactions and, even relaxing the assumption of hard-core bosons, the additional excited states of the system will be characterized by energies proportional to E C , thus much higher than the temperature.A similar conclusion holds also when considering the effect of Bogoliubov quasi-particles in the system, whose energy is higher than the SC gap, which corresponds to about 1.6K.
Besides the phononic environment studied in Ref. [62], another potential source of dissipation is provided by thermal charge fluctuations.This problem has been addressed in Ref. [53] for non-topological pumping protocols in driven SC devices.The charge noise affects the induced charges of the system n g, j and can be modelled, for instance, by a resistive thermal noise of the voltages V g, j parameterized by an effective resistance R.
The analysis of the impact of this kind of dissipation effects on the system dynamics can be performed via a generalized master equation approach.Also for this non-unitary kind of time evolution, the many-body Floquet representation can be efficiently adopted to calculate the expectation value of the observables of interest [63].For a weak coupling between the system and the thermal environment, the master equation approach models the dissipation through decay rates that depend on the environment via its spectral function, which can be taken Ohmic to model the noise of the gate voltages.
At low temperatures, an estimate of the relaxation rates associated to the noise of the induced charges results in Γ ≈ 2E g R(e 2 /h) 2C g /C self 2 [53] corresponding to Γ /E g ≲ 10 −3 for realistic parameters.
Within a Markovian approximation, this kind of dissipation dictates the relaxation of the system towards a Floquet steady state that assumes the form of a time-periodic density matrix diagonal in the Floquet eigenbasis.For temperatures considerably below the many-body gap, as the ones expected in experiments, this implies that the noise of the induced charges only provides weak corrections of the parameters N ν with respect to the unitary evolution.This applies in a straightforward way to the HH protocol, that relies exclusively on the driving of the n g, j parameters.A similar analysis, however, can be performed for the RM model, in which thermal noise can also affect the Josephson energies.In this case, the dissipation rates are proportional to E g (R/h) (δJ/δV c ) 2 , roughly corresponding to 10 −4 E g Re 2 /h based on the data in Ref. [51].Analogously to the charge noise, we expect that also these rates provide only small corrections to the populations N ν .
For a strong coupling with the environment, instead, the Markov approximation breaks down and alternative methods to study the dynamics are required.A first step in this direction has been considered in the study of the RM model coupled to a non-thermal phononic mode [64].
In general, however, the experimental results on the dynamics of SC circuits seem consistent with what expected in a Markovian description: we predict that the robustness against heating and dissipative effects of the presented pumping protocols in JJ chains will be analogous to what experimentally achieved for non-topological pumping schemes in different SC circuits [33,34], which allow for achieving pumping in the non-equilibrium Floquet steady state.This provides an important difference from ultracold atom setups, which are often more susceptible to heating given their considerably stronger isolation and where topological pumping is usually achieved only during a transient time window.Recent experiments with ultracold potassium atoms, however, allowed for the observation of Thouless pumping in a Rice-Mele-Hubbard protocol over more than a hundred cycles [65].
In conclusion, we expect that by operating the presented pumping protocols in standard conditions for superconducting circuits, such that the temperature of the environment is considerably lower than the relevant energy gaps, the coupling with a colder environment can actually assist against non-adiabatic effects [62].A beneficial effect of dissipation to improve the transport quantization of the RM pumping scheme has also been observed in experiments based on plasmonic waveguide arrays, in which it was shown that time-periodic and spaceperiodic dissipation can lead to the restoration of quantized transport for nonadiabatic driving conditions [66].

Conclusion
Motivated by recent advances in fabrication techniques, we numerically investigated possible driving protocols to implement topological Thouless pumping in short 1D arrays of tunable Josephson junctions.We considered the JJAs to be in a Coulomb-blockaded regime, where the charging energy of each SC island and the SC gap are the dominant energy scales, allowing us to approximate CPs as hardcore bosons.To study quantized transport in the long-time regime, we connected the array to two grounded SC leads, breaking the conservation of particle number and allowing for the observation of currents in the resulting non-equilibrium steady states.In particular, we used Floquet theory to extract the pumped charge at a small but finite driving frequency in the limit of an infinite number of driving cycles.
We focused on two prototypical models for topological quantum pumping, the periodically driven Rice-Mele and Harper-Hofstadter models.For both, we analyzed the role played by the coupling with the SC leads and their phase bias, as well as the effect of nearest-neighbor interactions originating from cross-capacitance between the SC islands.These ingredients are specific to the solid-state implementation of Thouless pumping with JJAs, and they extend the recent analyses of the role of interactions in Floquet topological insulators inspired by ultracold-atom experiments [23,28,40,59,65,67] which, in contrast, are closed systems with a fixed number of particles such that topological pumping is embodied in the motion of their center of mass.For the hybrid SC-SM devices we analyzed, instead, the understanding of the role of the external SC leads is essential for a successful experimental realization of topological quantized transport and the determination of the resulting currents.
Both the RM and the HH models display remarkable robustness with respect to nearestneighbor interaction, which does not affect their transport properties, even when it becomes larger than the Josephson tunneling.Moreover, the coupling to the leads helps to stabilize transport as it gaps out the zero-energy edge modes, which would otherwise appear in an open chain, and allows for an adiabatic evolution at sufficiently low driving frequency.
In the RM model, finite-size effects are most easily suppressed in a clipped driving regime where the modulation of the Josephson coupling is tuned such that the SM substrate below the Josephson junctions is depleted at half periods, effectively dimerizing the chain.However, the implementation requires simultaneous tuning of both the Josephson couplings and the charge induced on each island, resulting in a more complicated experimental protocol.
The HH model, on the other hand, requires fewer control gates as the Josephson couplings are constant in time.The drawback is that for small system sizes (M = 6, 9 islands) we find that it is necessary to dynamically average over the phase difference of the external leads to obtain a good quantization of the pumped charge.This behavior appears analogous to previous proposals to achieve quantized pumping in short Josephson junction chains [52].However, even at constant phase bias, the discrepancies from the ideal quantized case rapidly decrease with the system size.Moreover, the HH model displays a richer topological phase diagram since the number of bands with nontrivial Chern number depends on the effective magnetic flux Φ which in turn is determined by the position-dependence of the on-site energy modulation.In the simplest scenario, where Φ = 2π/3, the tuning of the average induced charge ng on the SC islands allows for control of the chemical potential.Tuning this to lie in different band gaps leads to a quantized current flowing in different directions.
Our results suggest that quantized Thouless pumping is experimentally accessible with the recently developed SC-SM JJAs, paving the way for a new generation of experiments investigating topological transport in Floquet systems as an alternative to ultracold-atom approaches.However, there are several further effects worth investigating to improve our analysis.First, one could relax the hardcore boson approximation and allow for a standard onsite repulsion, which has been studied for the RM model with periodic boundaries [39][40][41]68].We expect that an intermediate ratio E C /E J ≳ 3 is beneficial for experimental realizations because it would mitigate the effects of disorder in the induced charges without driving the system away from its insulating phases.In this situation, our estimates indicate that driving frequencies Ω ≳ 300MHz result in currents I ≳ 10pA without introducing excessive non-adiabatic effects.Another important element to consider is the possible presence of longer-range interactions.Indeed, the inverse of the capacitance matrix is approximately tri-diagonal only if the selfcapacitance of the SC islands is much larger than their cross-capacitance.While this can be a reasonable approximation, depending on the geometry of the JJA, it is important to examine the pumping robustness when it breaks down.Far from being necessarily a problem, longer-range electrostatic interactions can be exploited to study the transport of fractions of CPs pumped at each cycle [60].Finally, solid-state devices are subject to many decoherence channels, either due to quasiparticle poisoning, incoherent tunneling, charge noise, or scattering from crystal defects and imperfect interfaces.While precise modeling of these effects is prohibitive, an interesting future perspective is investigating the robustness of quantized transport with respect to simpler decoherence mechanisms, such as local dissipation or dephasing, similarly to Ref. [62], or incoherent coupling between the external (superconducting or normal) leads and the SC islands of the JJA.In particular, the analysis in Ref. [62] suggests that dissipations effect are beneficial to reduce non-adiabatic effects when operating in a cold environment.Based on our estimates, we expect indeed that Thouless pumping is observable in the Floquet steady state achieved under realistic conditions; this would provide a considerable advantage with respect to cold atom quantum simulators in which heating is often one of the main limiting factors.
Since the system is always in an insulating state, in the thermodynamics limit is possible to average over the phase twist φ at the boundary without affecting the expectation values of local observables [57].Thus we modify Eq.(A.The integrand is easily recognised as the Berry curvature B 0 (t, φ) = 2Im ∂ φ Φ0 (t) ∂ t Φ0 (t)〉 and its integral over the torus [0, t] × [0, 2π] gives, indeed, the first Chern number of the ground state manifold.Equation (A.7) is thus the equivalent of the TKNN formula [69] for particle transport with a periodic driving, thus showing the same topological origin of the IQHE and quantized Thouless pumping.Since, in the adiabatic limit, the time-evolved Floquet eigenstates converge to the instantaneous eigenstates, Eq. (A.7) shows the convergence of Q in Eq. ( 22) to a quantized value.

B Breakdown of the hardcore regime in the Rice-Mele protocol
Throughout the paper we employed a hardcore-boson approximation to describe the dynamics of Cooper pairs close to the charge degeneracy point n g = 1/2.This approach was motivated by the observation that, if E C is sufficiently strong, we expect only two charge states per SC island to influence the system dynamics when the induced charge is modulated around semi-integer values.This approximation is expected to provide reliable results only when the charging energy is sufficiently stronger than the Josephson energy scales.To test the robustness of the hardcore approximation, here we briefly investigate the dynamics of the RM model in a larger Hilbert space, which includes instead N cut charge states per SC islands.Because the computational cost increases exponentially as N M cut , M being the number of islands, we focus our attention on N cut = 4.The allowed local charge states are now labeled by N j = 0, . . ., 3 and we aim at driving the islands around the charge degeneracy point between states N j = 1 and N j = 2, therefore we set ng = 1.5.The hardcore-boson approximation is expected to fail when the dynamics leads the system out of the N j = 1, 2 subspace, such that also the N j = 0, 3 charge states get considerably populated.
To monitor the effectiveness of the hardcore approximation and the breakdown of the pumping for low values of the charging energy, we investigate the pumped charge over one period of the driving and the instantaneous many-body gap for M = 6 islands, as we change E C .All other simulation parameters are fixed: E C C = E L = δJ = J 0 , ω = 0.05J 0 , and the initial state is taken to be the ground state of the system at t = 0.
Figure 13(a) shows the pumped charge during the first driving period for different values of E C and compares the behavior for N cut = 4 (solid lines) and the hardcore approximation N cut = 2 (dashed lines).The data For E C = 3J 0 and 4J 0 (the value used in the main text) do not display any appreciable difference between the two cutoffs, suggesting that the hardcore-boson approximation is valid and well describes the system dynamics in this regime.As E C is lowered, larger discrepancies start to appear and the two behaviors becomes markedly different as the Josephson energy dominates over the charging energy (E C < J 0 ).Here, our simulations show nonphysical results as they indicate the emergence of a quantized pumping in the opposite direction with twice the Chern number.This is due to the spread of the charge over all the available states, as highlighted by Fig. 13(c), where the upper bound N j = 3 is saturated.A larger cutoff N cut is therefore needed to model accurately the dynamics, but it becomes unfeasible to study the time evolution with exact diagonalization.For comparison, Fig. 13(b) shows the density profile when the hardcore-boson approximation is valid and the average charge oscillates between 〈 Nj 〉 = 1, 2. Further insight on the approximation breakdown is provided by the instantaneous manybody gap E g as a function of time, shown in Fig. 14.As E C becomes smaller, the difference between N cut = 4 (blue solid lines) and N cut = 2 (orange dashed lines) increases; in particular, the system displays two gapless points (t = τ and t = τ/2) when E C = 0.5J 0 and N cut = 4, signaling the failure of the hardcore-boson approximation.This behavior is consistent with the fast spread of the charge over all available states in Fig. 13, whereas, for larger values of E C , the sectors away from N j = ng ± 0.5 are protected by the finite many-body gap.

Figure 4 :
Figure 4: Spectrum of the Harper-Hofstadter model as a function of momentum and time at Φ = 2π/3 for δµ = 0.4E J and μ = 0.

Figure 5 :
Figure5: Pumped charge at infinite time Q ∞ as a function of the phase difference φ between the external superconducting leads for M = 4, 6, 8 islands.The data correspond to the following parameters: δJ = 0.7J 0 , δµ = 1.5J 0 , E L = 2J 0 , ω = 0.05J 0 .The table represents the deviations from perfect quantization of the phase averages pumped charge Q as a function of the systems size M .

Figure 6 :
Figure 6: Pumped charge Q ∞ of the RM model for 20 discretized values of φ (blue points) as a function of the Josephson energy modulation δJ.The red crosses correspond to Q.The phase dependence and the charge quantization considerably improve approaching the clipped modulation.The driving frequency is ω = 0.05J 0 .

Figure 7 :
Figure 7: Comparison of the many-body Floquet eigenenergies of the RM model for (a) small sinusoidal modulations δJ = 0.5J 0 (b) and clipped modulation δJ = 1.2J 0 of the Josephson tunnelings.The dot colors represent the overlap N ν (φ) of the initial ground state with the Floquet states.In the clipped regime ∂ φ E 0 is almost constant, leading to a pumped charge almost independent on the phase.The data refer to a system with M = 6 islands, ω = 0.05J 0 , and E L = 2J 0 .

Figure 8 :
Figure 8: (a) Pumped charge Q ∞ (φ = 0) and many-body energy gap E g as a function of E L for δJ = J 0 (threshold between clipped and sinusoidal modulations).The data are obtained in a system with parameters δµ = 1.5J 0 , ω = 0.05J 0 , and M = 6.(b) The many-body Floquet eigenenergies of the Rice-Mele model for δJ = J 0 and E L = 0.5J 0 .In this situation, multiple Floquet states overlap significantly with the initial ground state resulting in a non-quantized current.Data are obtained for 6 islands and a driving frequency of ω = 0.05J 0 .

Figure 9 :
Figure 9: Pumped charge Q ∞ (φ = 0) as a function of the coupling with the leads E L and the nearest neighbor interactions E C C for the RM model with 6 islands.The data are obtained with δJ = J 0 , ω = 0.05J 0 and δµ = 1.5J 0 .

Figure 10 :
Figure 10: Pumped charge Q as a function of the coupling with the leads E L for different values of the nearest-neighbor interaction.The data are obtained with parameters M = 6, E C = 4E J , δµ = E J , and ω = 0.01E J , in the gapped topological phase at filling 1/3.

Figure 11 :
Figure 11: Winding of the quasienergies of the many-body populated state as a function of the lead phase difference φ for M = 6, E C = 4E J , E C C = 0.2E J , δ µ = E L = E J and ω = 0.01E J .(a) Winding in the gapped phase at filling 1/3 (n g = 0.4).We observe that the quasienergy winds three times in the positive direction and twice in the negative direction, corresponding to a Chern number C 1 = 1 of the lowest band.(b) Winding in the phase at filling 2/3 (n g = 0.6).The winding is opposite with respect to panel (a).

Figure 12 :
Figure 12: Pumped charge in the hardcore boson HH model.(a) Dependence of the pumped charge on the phase difference φ for chains of M = 6 and M = 9 islands around the transition ng = 0.35 (orange) and in the topological insulating phase ng = 0.4 (teal and maroon).(b) Averaged pumped charge Q for M = 6 islands.The blue curve corresponds to vanishing nearest-neighbor interactions.The orange curve is obtained for E C C = E J /2. Interactions extend the Floquet topological phases.Inset: fluctuations of Q close to the transition between the metallic and the topological insulating phases for E C C = 0.

Figure 13 :Figure 14 :
Figure 13: (a) Pumped charge during the first period of evolution for different values of E C (color-coded) for a RM setup with M = 6 islands.Solid lines correspond to N cut = 4 and dashed lines to N cut = 2. (b-c) Time-resolved local charge for E C = 4J 0 (b) and E C = 0.5J 0 (c).All data are taken for φ = 0.