Cooper-pair transistor as a minimal topological quantum circuit

The outlook of protected quantum computing spurred enormous progress in the search for topological materials, sustaining a continued race to find the most experimentally feasible platform. Here, we show that one of the simplest quantum circuits, the Cooper-pair transistor, exhibits a nontrivial Chern number which has not yet been discussed, in spite of the exhaustive existing literature. Surprisingly, the resulting quantized current response is robust with respect to a large number of external perturbations, most notably low-frequency charge noise and quasiparticle poisoning. Moreover, the fact that the higher bands experience crossings with higher topological charges leads to all the bands having the same Chern number, such that there is no restriction to stay close to the ground state. Remaining small perturbations are investigated based on a generic Master equation approach. Finally, we discuss a feasible protocol to measure the quantized current.


I. INTRODUCTION
Topological phases are an important research topic in condensed matter physics [1] most notably with the goal to realize inherently protected quantum computing [2]. The most common approach is to search for topological transitions in the band structure of the materials themselves [3][4][5], such as in topological insulators [6][7][8][9][10], Chern insulators [11][12][13], Weyl and Dirac semimetals exhibiting Fermi-arc surface states [14][15][16][17][18][19], or topological superconductors hosting Majorana fermions [20][21][22][23][24][25][26][27][28][29][30][31][32]. However, the realization of topological materials turns out to be challenging due to various reasons, such as a lack of tunability, detrimental effects of impurities [33][34][35][36], or quasiparticles in Majorana-based systems [37][38][39][40][41]. A further challenge concerns the direct observability of the topological invariant; often, the topological phase is only indirectly measured through the density of states, e.g., via ARPES [17] or STM [42]. This is why alternatives are actively researched, where the topological phase is encoded in other degrees of freedom. Circuit lattices [43][44][45][46][47][48][49][50][51][52][53] may form metamaterials where topological numbers are defined through the lattice degrees of freedom, which however requires the control of a large number of circuits. Topological materials are also very straightforwardly simulated when considering the space spanned by the control parameters of superconducting qubits [54][55][56], where it remains however unclear, how physics related to protected edge states may be observed. Such limitations may be circumvented by recently proposed topological transitions in multiterminal Josephson junctions [29][30][31][57][58][59][60][61][62][63][64][65] which give rise to topological phases even when using only trivial materials [66][67][68][69][70][71][72][73][74][75][76]. Here, Weyl points are found in the space of superconducting phase differences acting as quasimomenta, and a Chern number can be directly accessed through a quantized transconductance [66]. While topological transitions in transport degrees of freedom offer a * t.herrig@fz-juelich.de Figure 1. Circuit of the Cooper-pair transistor and quantized dc current responses. The two Josephson junctions are described by the energies EJL and EJR. We apply a voltage V between the left and right lead with superconducting phases φL and φR, which drives the phase difference adiabatically according to the second Josephson relation. The charge and phase of the superconducting island are described by the conjugate variables N and ϕ. The linearly time-dependent gate voltage Vg (t) induces a dc current I ind ∝Vg which flows a) entirely through the right junction if EJR > EJL and c) entirely through the left junction if EJR < EJL. b) + d) The dc parts of the expectation values of the currents coming from the right and left lead, IR and IL , respectively, are depicted under adiabatic conditions as a function of EJR/EJL. promising new approach, the proposal in Ref. [66] came with the experimental complication of needing small multiterminal junctions containing only a few channels, which are at the same time strongly tunnel-coupled. An important simplification was recently proposed by means of Weyl points in standard SIS junction circuits [77][78][79][80][81]. However, these proposals require a control of the offset charges on the order of a single electron charge e which may be experimentally challenging [82][83][84], unless offsetcharge feedback loops are employed [85].
Here, we consider the Cooper-pair transistor, consisting of two tunnel junctions with a superconducting island arXiv:2012.10655v4 [cond-mat.mes-hall] 13 May 2022 in between; see Fig. 1. Although this circuit has been studied to a great extent [86][87][88][89][90][91][92][93][94][95][96][97][98][99], the Weyl points it exhibits in its band structure have, to the best of our knowledge, not yet been discussed. Additionally to the phase difference across the two junctions, we use the island offset charge to define the Chern number which gives rise to a topological phase transition when the asymmetry of the Josephson energies changes its sign. This Chern number leads to a quantized dc current response into a particular lead when driving both the offset charge and the phase difference time dependently (see Fig. 1).
Remarkably, the quantization of the current response is insensitive to low-frequency offset-charge noise; in fact, it is actually beneficial for the convergence of the response. Furthermore, we find that the Chern number is insensitive to quasiparticle poisoning. Moreover, the crossings in the higher bands occur via Weyl points with higher topological charges. As a result, each band carries the same Chern number such that it is not required to be in the ground state to observe the effect. This is a surprising exception because usually in quantum systems, the ground and excited states exhibit different topological numbers, a fact which recently lead to the effort of generalizing topological phase transitions to systems out of equilibrium [100][101][102][103][104][105][106][107][108][109]. Motivated by the above remarkable protection, we analyze the influence from the environment more closely, by means of a generic master equation. Based on this, we expect that the leading-order deviation of the quantized current remains small. Finally, we discuss an experimentally feasible protocol to measure the quantized current. This protocol is to some extent reminiscent of Cooper-pair pumps [79][80][81][110][111][112][113], with the notable difference that previously studied pumps do suffer from quasiparticle poisoning [113].
Compared to the above vast existing literature on topological quantum systems, the effect we study here has the following advantages. First, we do not need any topological materials nor a network of coupled circuits to find topological phase transitions -a single circuit made of regular s-wave superconductors suffices. What is more, compared to similar recent proposals, this circuit has only two terminals and works already by means of standard SIS junctions. In fact, as we will show below, all independent system parameters entering the Hamiltonian play a crucial role for the observed topological effect, leaving no "spare" degrees of freedom. It is in this sense that we refer to our system as a minimal topological circuit. This paper is structured as follows. In Sec. II, we review the Hamiltonian of the Cooper-pair transistor and discuss its topological features. Afterwards, in Sec. III, we show how to access the Chern number by varying specific system parameters in time. In Sec. IV, we start with a short discussion of the robustness of the resulting quantized dc current response with respect to the most common perturbations, followed by the introduction of a generic open-system description to discuss the leadingorder perturbation to the otherwise protected quantization of the current response. The concrete experimental Displayed are the lowest four energy bands for φ = φR − φL = π as a function of the offset charge Ng for three cases where the junction asymmetry changes from EJR < EJL over EJR = EJL to EJR > EJL. In the symmetric case, b), one can see the Weyl points as band crossings, each associated with a topological charge C. In the asymmetric cases, a) and c), each band n can be assigned a Chern number CR,n, which is zero in the trivial phase and changes only when passing through a Weyl point. Due to the topological charges increasing by one with each higher Weyl point, the Chern numbers for the different bands all change by the same value.
realization will be the topic of the final Sec. V, where we propose a practical measurement scheme.

II. THE CIRCUIT AND ITS TOPOLOGY
We here investigate the topological properties of the Cooper-pair transistor, consisting of two Josephson junctions connected in series with energies E JL and E JR (see Fig. 1a) forming a charge island, which is capacitively coupled to a gate voltage V g . The Hamiltonian is given by The first term describes the charging energy of the superconducting island with the Cooper-pair number operator N and E C = (2e) 2 /(C L + C R + C g ), where C L , C R , and C g are the capacitances of the two junctions and the gate capacitor, respectively. The gate voltage induces the charge offset The phase operator ϕ is canonically conjugate to the number of Cooper pairs, such that ϕ, N = i. Due to charge quantization, we can write the Hamiltonian in the discrete charge basis with N = N N |N N | and e i ϕ = N |N N − 1|. We denote the eigenenergies and eigenvectors of H as n and |n , respectively, which correspond to the standard Mathieu functions [86,114]. The energy spectrum is shown in Fig. 2.
While this system has already been extensively studied [86][87][88][89][90][91][92][93][94][95][96][97][98][99], it has to the best of our knowledge not yet been explicitly remarked that it harbors nontrivial Chern numbers in the base space given by the charge offset N g and the phase difference φ = φ R − φ L , with the Berry curvature B α,n = −2 Im ∂ φα n ∂ Ng n (α = L, R). In fact, the transistor thus simulates a Chern insulator, where the parameter pair (N g , φ) acts as the Brillouin zone on a 2D torus [115]. Let us now explain the origin of the nontrivial Chern number in the remainder of this section. First, we note that the Chern number defined in Eq. (2), as well as the corresponding Berry curvature, carry the index α = L, R, indicating that there are actually two distinct Chern numbers (Berry curvatures) for one and the same band. The reason for this is the following. The eigenenergies only depend on the total phase difference, n (φ), such that ∂ φ L n = −∂ φ R n . This stems from the fact that φ α can be "eliminated" by the unitary transformation In other words, the eigenenergies are insensitive to a unitary change of basis. The eigenvectors on the other hand are not, and therefore the two Berry curvatures B L,n and B R,n are nontrivially related via ( This is a consequence of the current conservation law, α ∂ φα H = i H, N , which we will come back to in Sec. III. Note that while B α,n depends thus on α, it depends explicitly on only the phase difference, B α,n (N g , φ R − φ L ), which is why it is sufficient to integrate over φ in Eq. (2).
Based on Eq. (3), we can also relate the Chern numbers for different α, On a formal level, this difference between the Chern numbers for different α is a simple consequence of the definition in Eq. (2). On a physical level, this difference will become meaningful in Sec. III, as the difference of measuring the current from the left or from the right contact. In fact, the physics will be engrained in the very structure of the Berry curvature: a current measurement into contact α (∂ φα ) as a response to a driving of N g (∂ Ng ).
The Chern numbers, as defined above, must be quantized, in accordance with standard literature [116]. One way to explicitly calculate their actual values, would be by inserting the Mathieu functions into the Berry curvature and performing a numerical calculation of the values of the Chern numbers. However, we here resort to a simpler way which is inspired by Ref. [66].
Namely, the nonzero Chern numbers are a consequence of Weyl points (that is, topologically protected band crossing points) appearing in the 3D space given by (N g , φ, E JR /E JL ). These points appear for symmetric junctions, E JR = E JL , at φ = π + 2πm, where m ∈ Z. Here, the Josephson energies for the left and right junctions cancel in the Hamiltonian, such that only the charging energy remains, H = E C N +N g 2 /2. Therefore, the Weyl points simply represent the crossings of the shifted parabolas for different charge states on the island. Indexing the ground state as n = 0 and the excited states with n > 0 in ascending order, one can state that the crossings between band n and n + 1, for n odd (even) occurs at N g being (half) integer; see Fig. 2b. Hence, when we tune from E JR < E JL to E JR > E JL (see Figs. 2a and c), the Chern numbers for the different bands [Eq. (2)] change according to the topological charge (or the chirality) C of the corresponding Weyl points (see Fig. 2b). Importantly, while the Weyl points connecting the bands n = 0 and n = 1 are regular Weyl points with topological charge C = +1, the Weyl points connecting arbitrary higher bands n and n + 1 have in fact a higher topological charge, C = + (n + 1), as indicated in Fig. 2b. We will call such a point a multiple Weyl point, which can be considered as a merger of n + 1 regular Weyl points, each with charge +1, as we explain in a moment. Since each band n experiences a change in its Chern number by subtracting the topological charge of the Weyl point connecting to band n − 1 from the topological charge of the Weyl point connecting to n + 1, the Chern numbers of all the bands are the same. For C R,n , it follows that from a completely trivial spectrum for E JR < E JL , where all C R,n = 0 (Fig. 2a), we go to a spectrum where all bands have the same nonzero Chern number C R,n = +1, for E JR > E JL (Fig. 2c). Vice versa, for C L,n , we find C L,n = +1 for E JR < E JL , and C L,n = 0 for E JR > E JL . Thus, it is only meaningful to consider a band as "topologically trivial" if we refer to a specific Chern number, C L,n or C R,n , because one is zero and the other one is nonzero for every junction asymmetry, E JR /E JL ≶ 1, independent of the band n. The Chern number being the same for all bands is a remarkable feature, and renders the observation of the Chern number insensitive to whether or not the system is in the ground state (e.g., when including finite temperature); see also the discussion in Sec. IV.
Let us now discuss the physical origin of the multiple Weyl points. We here provide an explicit discussion for the lowest two Weyl points, which amount to the topological charges of +1 and +2. First, regarding the regular Weyl point connecting the ground and first excited state, the band crossing involves two charge states which differ by only one Cooper pair, |N − 1 and |N . Tuning the parameters close to the crossing, N g = −N + 1/2 + δN g , E JR = E JL + δE J , and φ R = π + δφ while at the same time φ L = 0 [117], we find the approximate Hamiltonian as derived in Appendix A 1 with x = δE J /2, y = E JL δφ/2, and z = E C δN g /2, and the Pauli matrices acting on the charge subspace This is the standard form of the Weyl Hamiltonian with the topological charge C = +1.
As for the double Weyl point, with charge +2, we have to tune to N g = −N + δN g , while the other small parameters (δE J and δφ) are defined as above. Here, the relevant subspace close to the crossing involves the charge states |N − 1 and |N + 1 . Importantly, here it is impossible to gap the two bands with the lowest-order Cooperpair tunneling process, because N − 1| e ±i ϕ |N + 1 = 0. Therefore, we need to go to higher-order processes involving the tunneling of two Cooper-pairs via virtual charge states, which can be done by means of a Schrieffer-Wolff transformation; see Appendix A 2. We find the Hamiltonian of the following form, The fact that the Weyl point, here, has a topological charge of +2 can be shown in different ways. One could in principle define a Berry curvature in the 3D space (x, y, z) and compute explicitly a closed surface integral enclosing the Weyl point. A more elegant and instructive way is, however, to add a cotunneling term E (2) J cos 2 ϕ to the full Hamiltonian, Eq. (1), which may originate from higher-order tunneling processes in the SIS junction [118][119][120]. Note that this cotunneling term is by no means relevant for our considerations. It however serves as a neat mathematical trick to visualize the topological charge. Namely, this term introduces a shift c = E As a consequence, the Weyl point at (0, 0, 0) for c = 0 splits into two regular Weyl points with topological charge +1 at (± √ c, 0, 0) for finite c. Consequently, without the splitting, it must have the topological charge C = +2 and, thus, is a double Weyl point. This proof (not explicitly shown here) can be extended to higher bands, where three or more regular Weyl points are merged, giving rise to a triple or higher multiple Weyl point, due to the gapping of the bands being third or higher order in Cooper-pair tunneling processes, respectively.
To conclude, let us note one subtle difference between our notion of a Chern number compared to the more common definition within solid state theory, where the Chern number is usually connected to the single-particle band structure. Here, our bands are already in a manybody formulation, such that there is no need to associate occupation numbers to the individual bands.

III. QUANTIZED CURRENT RESPONSE
We now show that the above discussed nonzero Chern numbers lead to a directly measurable effect, which is a quantized, directed current response, that is, a dc current flowing either precisely to the left, or precisely to the right (depending on the junction asymmetry), as depicted in Fig. 1. This effect emerges when driving N g and φ time dependently. The driving of the superconducting phase difference is accomplished by means of applying a voltage, whereas the gate-induced offset charge is linearly ramped up (or down) with a constant ramping speedṄ g . Let us emphasize that there is an arbitrary number of possibilities to satisfy Eq. (7). Strictly speaking, since these different choices are related through a time-dependent unitary transformation U (via H → U H U † ), they do not give rise to equivalent Schrödinger equations, due to the additional term −i U ∂ t U † . However, this leads merely to a shift in the initial condition on N g (t), which is irrelevant for the here considered dc current response (due to the time-averaging). Our subsequent results are formulated independent of this choice. Before proceeding, let us comment on one important point. Of course, the ramping up of N g with a constant ramping speed can in reality not be exerted for unlimited time, as the transistor would eventually break. However, as we show in Sec. V, this is no actual limitation, as this problem can be easily circumvented by choosing an appropriate driving protocol.
As we show now, the quantized current response resulting from the time-dependent driving is in close analogy to the proposal in Ref. [66]. In the four-terminal setup of Ref. [66], voltages were applied to two different contacts, resulting in a pure φ-driving and a resulting dc transconductance. Similar proposals have very recently emerged in pure SIS junction circuits [77,78]. Here, on the other hand, we have a two-terminal device with only one independent phase difference φ = φ R − φ L and driving in the "mixed" parameter space (N g , φ).
To proceed, we consider the dynamics of the system for slow, adiabatic driving. In this limit, the time-dependent Schrödinger equation i∂ t |ψ n (t) = H(t) |ψ n (t) has the solutions [121] |ψ under the assumption that, at the initial time t 0 , the system was in the eigenstate |n(t 0 ) . Here, H(t) |n(t) = n (t) |n(t) denote the instantaneous eigenbasis at time t. The time evolution gives rise to the (here irrelevant) dynamical phase α n (t) = −i t t0 dt n (t ) [122]. Note that for simplicity, we refrain from explicitly adding the time argument (t) from now on. This adiabatic approximation is valid for [123] | m| ∂ t |n | | m − n |, (n = m), which in our case requires Ṅ g , |V | inf | m − n |. Also note that adiabaticity is a standard requirement (see also Refs. [66,69,77,78]) and that nonadiabatic effects like Landau-Zener transitions [124] are exponentially suppressed away from the degeneracy point.
As indicated above, we are interested in the expectation values of the currents into the system from the left and right contacts in the presence of the drive. The current operators are formally defined as The reason that we have to consider both the left and the right current separately, is that the time-dependent driving of the gate charge N g produces a finite dc displacement current, such that the current expectation values do not simply add up to zero, as it would be the case for the stationary system. Instead, we have to carefully keep in mind the current conservation law where the right-hand side is nonzero, due to H not commuting with the island charge N . The expression˙ N is, of course, to be understood in the Heisenberg picture.
We can now compute the expectation values of the currents, by inserting the wave function given in Eq. (8). We find α,n is a correction term of first order in the driving parameters.
With the above result, we can now discuss the very important issue of current conservation, to understand how the expectation values of the left and right currents are related. The sum of the zero-order contributions α I (0) α,n has to vanish due to the eigenenergies depending only on the total phase difference (as discussed in Sec. II). Importantly, however, the first-order term I (1) α,n gives rise to a Berry curvature B α,n , as it appears in Eq. (2), due to the driving of the offset charge N g ; and since the Berry curvatures for different α being nontrivially related via Eq. (3), the sum of the first-order contributions is nonzero. This physical consequence is reflected in the current conservation I R + I L ψn = 2e ∂ t N n , in accordance with Eq. (10).
In the same spirit as in Ref. [66], we now find that when averaging the currents over long times I α,n = lim τ →∞ τ 0 dt τ I α,n (dc limit), the zeroth-order contributions I (0) α,n cancel and the first-order ones I (1) α,n average out to give the Chern numbers C α,n from Eq. (2), providing the result I α,n = −2eṄ g C α,n . This is due to the currents being periodic in N g and φ, such that the long-time integral will eventually converge into an area integral over the Brillouin zone of the (N g , φ)-plane. Importantly, the presence of the Weyl points and the resulting nontrivial Chern numbers (as discussed above and shown in Fig. 2) lead to the quantized dc currents into the system This central result is also visualized in Fig. 1. Namely, we find that the current injected into the system through the ramping of the gate voltage, 2eṄ g , is completely redirected to the right (left) contact, if E JR is larger (smaller) than E JL -crucially, irrespective of the precise ratio between E JR and E JL . As we see, the difference between C α,n for different α corresponds to current measurements at different contacts α, which are not equal due to the displacement current.
We note that the plots in Fig. 1 b) and d) are idealized in the sense that it is purely hypothetical to stay adiabatic in the vicinity of the symmetric point, E JR = E JL , where we necessarily come close to a degeneracy point. At such a point, Landau-Zener transitions [124] cannot be avoided. These transitions, and their importance for the observation of the Chern number have already been discussed elsewhere, for SNS type junctions (see Ref. [69]), which is why we do not repeat a similar effort here. Overall, it can be expected that the step in Fig. 1 b and d will be "rounded" due to Landau-Zener transitions, where the broadening decreases when reducing the driving frequency of N g and φ. We also note the importance of applying a voltage V across the two contacts. If we were to modulate N g only, the integral of the zeroth-order term I (0) α,n ∼ ∂ φα n would not in general vanish, nor would the integral of the Berry curvature extend over the entire Brillouin zone of (N g , φ). In the absence of any bias voltage, the total injected current through the gate drive will be partitioned to the left and right contact with proportions depending on many system parameters. The perfect directing of the injected current to exclusively either the left or the right, independent of the parameter details, is a pure topological effect, requiring the combined modulation of the two parameters N g and φ.
Finally, let us get back to the titular notion of a minimal topological circuit. In the Hamiltonian, there are overall the independent parameters E C , E JR − E JL , φ, and N g [125]. The last three out of these four parameters provide the base space in which the Weyl points are defined. Finally, E C guarantees that the energy spectrum is gapped, such that it is meaningful to consider discrete bands with individual Chern numbers. Thus, all involved parameters play an indispensable role for the observed topological effect.

IV. STABILITY WITH RESPECT TO EXTERNAL PERTURBATIONS
In analogy to Ref. [66], the convergence of the dc current to the values in Eqs. (12) and (13) requires the driving frequenciesṄ g and 2eV either to be incommensurable or to come with a sufficient amount of low-frequency noise to make sure that the entire Brillouin zone is covered for sufficiently long times τ . Surprisingly, this means that, here, low-frequency offset-charge noise actually helps for the observation of the Chern number, instead of perturbing it. This is a considerable advantage with respect to other recent proposals [77,78] which rely on a control of the offset charge on the order of an elementary charge e during the entire integration time τ , which seems challenging given the experimental evidence for offset-charge noise [82][83][84].
A further important point concerns quasiparticle poisoning. Quasiparticles appear to be much more numerous than what should be expected in equilibrium [126,127] and induce stochastic switches between states of different parity. They are harmful for a large number of quantum devices such as Cooper-pair boxes [128,129], transmon and fluxonium qubits [130], Flux-qubits [131], Majoranabased qubits [37][38][39][40][41], or Cooper-pair sluices [113]. Also, the observation of transport Chern numbers defined purely in φ-space [66,77,78] is hampered by quasiparticle-induced parity flips, since the topological numbers differ in different parity sectors. The Chern number we consider here, however, is the same, independent of the parity. Namely, in our formalism, parity flips can be accounted for by simple shifts of N g by half an integer, N g → N g ± 1/2. Due to the periodicity of the Berry curvature in N g -space, it follows that 1 0 dN g B α,n (N g ± 1/2, φ) = 1 0 dN g B α,n (N g , φ), demonstrating the insensitivity of the Chern number, Eq. (2), on fermion parity. Moreover, as already indicated above, the Chern numbers for all the bands n are the same, due to the higher Weyl points having higher topological charges, such that the dc current does not depend on n, see Eqs. (12) and (13). Therefore, the here predicted effect is likewise not sensitive to finite temperature occupations of higher energy bands, which is a further advantage over the proposals in Refs. [66,77,78].
Encouraged by these striking facts, we now consider the effects of the environment in more detail. As we will show, stochastic transitions induced by the environment will introduce what we expect to be a small leading-order correction in the current response. Such stochastic transitions may be of various origins. In particular, the aforementioned beneficial effects of fluctuations in N g or φ α are restricted to low-frequency fluctuations, that is, in a frequency regime where the noise can be considered adiabatic, such that it does not give rise to stochastic transitions between different energy bands. Even though the finite-frequency power spectrum of the noise can be expected to be low, such transitions may still occur in reality and have to be taken into account. Likewise, quasiparticles may induce stochastic transitions within the time required to average the current signal.
Therefore, a detailed open quantum system description would have to encompass a large number of models to account for all perturbations. For instance, phase fluctuations are standardly described by an external impedance, which can be modeled by an ensemble of LCresonators [132]. Charge fluctuations are usually modeled via so-called two-level fluctuators [133][134][135]. However, such models have recently been put into question for 2D transmons, where a deviation from the typical 1/f -noise spectrum has been observed [84]. In order to avoid any dependence on such details, we here take into account the open-system dynamics as generally as possible, by means of a quantum master equation for the density matrix of the system ρ, The effect of the environment is described by the timelocal kernel W, which may in general be time-dependent (since the system is driven time dependently). In order to guarantee positivity of ρ for all times t, we assume that W can be cast into a Lindblad form (whose specific form however is, in a first approach, irrelevant). Apart from that, we assume only that in the absence of the time-dependent driving, the system will, up to small corrections, end up in the ground state as the stationary solution, ρ st ≈ ρ 0 = |0 0|. This is equivalent to assuming that the environment has a small temperature with respect to the band gaps, k B T inf | n − 0 |. Thus, we can understand W as a generic cooling mechanism. Note importantly that we do not assume the usually standard rotating-wave approximation (RWA) and instead keep processes in W which couple diagonal and off-diagonal contributions to the density matrix (that is, diagonal and off-diagonal with respect to the instantaneous eigenbasis of H). Strikingly, we find in the absence of the RWA a leading-order contribution which would otherwise be neglected and which is, at least nominally, of the same order as the Thouless result. Nonetheless, this correction can be shown to be small, due to other small factors; see below.
Let us now again focus on slow driving. For this purpose it is useful to cast Eq. (14) into the aforementioned instantaneous eigenbasis of H. We thus get Importantly, all the objects appearing in Eq. (15) differ from the ones in Eq. (14) by the time-dependent unitary transformation U (t), which changes the basis to the instantaneous eigenbasis. For example, for the density matrix, this corresponds to a mapping ρ → U ρ U † . Therefore, one should strictly speaking use different symbols for Eqs. (14) and (15) from which we will refrain for notational simplicity. As a consequence of the unitary transformation, the closed system dynamics receives an extra term to the ordinary L 0 • = H, • , denoted as For a consistent slow-driving approximation we have to consider the relationship between the three time scales W −1 , δL −1 , and L 0 −1 (where A is a suitably chosen norm to capture the magnitude of A). The regular closed-system dynamics scales with the instantaneous energy gaps L 0 ∼ | m − n | and the additional contribution due to the driving scales as δL ∼ m| ∂ t |n (both with m = n). Now, the regime of interest is L 0 −1 < W −1 < δL −1 . Thus, we explicitly expand the density matrix in orders of W and δL , where ρ (ν,µ) scales as δL ν W µ , and plug the result into the expectation value I α = tr I α ρ . The derivation of these results can be found in Appendix B. Note that in general, there may be specific types of system-reservoir interactions providing an additional contribution to the current, which can generically be taken into account by a "current kernel" W I [137], such that I α = tr I α ρ + tr W I ρ . For simplicity, we assume that the systemreservoir interaction is current conserving such that there is no such W I term. Since W > δL , the corrections of ρ due to the open-system dynamics are in principle dominant. However, in the absence of an explicit (timedependent) driving and the environment at equilibrium, the system cannot generate a finite dc current, such that the dc contributions of all I (0,µ) α = tr I α ρ (0,µ) must be zero [138]. Therefore, we need to consider at least first order in the driving parameter δL . As we will show in the following, the leading-order terms of the dc current will be the computation of which is detailed in Appendix B. Nominally, both of these contributions, I α,d , are of the same order (1, 0). However, they are of fundamentally different origin. As we will discuss now, the first term, I α,d , will provide the leading-order correction to the above Chern number term, and is due to a combination of driving and opensystem dynamics (even though it is nominally of zero order in W ), as we will see in a moment. Moreover, its existence is based on the kernel creating transitions between diagonal and off-diagonal matrix elements, which is not the case when applying the RWA.
As already foreshadowed in Eq. (17), in order to discuss the current it is useful to decompose each order, I (ν,µ) α = tr I α ρ (ν,µ) , as with Here, the projection superoperator P d is defined to project onto the diagonal sector (in the eigenbasis of H), giving n|(P d ρ )|m = δ nm n| ρ |m . Vice versa, P o projects onto the off-diagonal subspace, such that P d + P o = I, where I is the identity superoperator, leaving any input matrix unchanged. Thus, the first term, is arising from the density-matrix contributions that are diagonal (analogously to the zero-order term in the closed system, I α,0 = ∂ φα n ), while the second term is associated with the off-diagonal part of ρ, originating in the driving and the interaction with the environment. It can be brought into the form Here, we defined N α ≡ −i n,m n| ∂ φα |m |n m|. Evidently, N α can be formally related to the Cooper-pair number operator of contact α. However, note that care has to be taken with this interpretation. The contacts in the here considered model are macroscopically large and their actual charge operators do not have a well-defined expectation value, whereas N α is always well-behaved.
To avoid such unnecessary complications, we simply refer to it in the way it is defined: as the operator −i∂ φα expressed in the eigenbasis of H.
We now have to find the contributions to the density matrix by expanding it according to Eq. (16) and solving the Lindblad equation, Eq. (15), for leading orders in (ν, µ). The details of this calculation are shown in Appendix B. Let us first consider the contribution I α,o as defined in Eq. (21). When plugging in the solution ρ (1,0) , we find the term giving rise to exactly the same Berry curvature term as in the closed system; see Eq. (11). For I α,o one could now in principle compute higher-order terms, such as ρ (1,1) , to obtain corrections due to the open-system dynamics. However, as already indicated in Eq. (17), when including the impact of the environment, it turns out that there will actually be a correction nominally of zero order in W in the current contribution I α,d . Let us here discuss the origin of this surprising new term. When treating the closed system in Sec. III, we have given the solution of the diagonal part of the density matrix as P d ρ = |0 0|. However, strictly speaking, this is not the most general solution for the closed system. As a matter of fact, in the absence of any relaxation mechanisms, the most general solution for the diagonal part can be any mixed state P d ρ = n ρ nn |n n| with arbitrary ρ nn provided that ρ nn be constant in time. If all ρ nn are indeed constant in time, then Eq. (20) will average to zero when computing the dc component. In the open system, the kernel W accomplishes two things. On the one hand, it fixes the solution for the diagonal part of the density matrix, that is the ρ nn are now unique. On the other hand, these ρ nn now depend in general on time. Let us emphasize that this nontrivial behavior of the diagonal part is present only when including processes that couple the diagonal and the off-diagonal sectors in Win other words, it is a consequence of going beyond the RWA. Therefore, we find with Since ∂ t ρ (1,0) nn = 0 in general, the dc contribution of this current will be nonzero, marking it as the leading-order correction term. Note that the inverse 1/L 0 is defined only in the off-diagonal sector. However, the objects W ρ 0 and δL ρ 0 , on which this inverse acts, are both purely off-diagonal, which is also explained in the Appendix B. Similarly, W dd ≡ P d WP d has a zero eigenvalue, corresponding to a degenerate eigenspace including the ground state, ρ 0 = |0 0|, and any purely offdiagonal matrix. Therefore, its inverse is ill-defined. However, we can construct the inverse W −1 dd such that . For more details see Appendix B.
While Eqs. (23) and (24) are general under the assumptions made, for the sake of concreteness, we now compute them numerically considering background charge fluctuations [133][134][135], N g → N g + δ N g . Here, δ N g is an operator, whose dynamics is governed by an environment Hamiltonian, H env . The circuit Hamiltonian can now be approximated as follows, H(N g ) → H(N g + δ N g ) ≈ H(N g ) + V , where we find the interaction with the environment V = γE C N ⊗ δ N g . The total system plus environment is thus described by the Hamiltonian H(N g ) + V + H env . Let us now focus on the regime where E J (φ) E C with the Josephson energy E J (φ) = E 2 JR + E 2 JL + 2E JR E JL cos(φ) (which we refer to as transmon limit). This regime is of particular interest for the here considered topological effect, because the energy bands are flatter and therefore transitions between bands (via Landau-Zener) are suppressed. Note that we aim to remain in this transmon regime for all φ, which implies in addition, that the two junctions should be sufficiently asymmetric w.r.t. E C , |E JR − E JL | E C . Additionally, we assume ξ 1, with ξ ≡ inf{E JR /E JL , E JL /E JR }, to simplify the numerical calculations. In this transmon limit, we can describe our system as a damped harmonic oscillator (HO) with an energy spacing of ω 0 = E C E J (φ). The respective kernel is a well-known result from standard literature of open quantum systems [139], which can be computed by tracing out the environment degrees of freedom, resulting in with The four (real) correlation functions are where with λ = γ √ E C ω 0 . We observe that these correlation functions depend only on the transmon frequency ω 0 (and E C ), which in turn does not depend on N g and only very weakly on φ. Therefore, this already indicates that the last term of Eq. (24), i.e., the term with the time derivative ∂ t , is not the dominant term.
Due to microreversibility, the ratio of excitation and relaxation rates of two neighboring states is a Boltzmann factor Γ n→n+1 /Γ n+1→n = (∆ − ζ)/(∆ + ζ) = e −βω0 , with Γ n→n ≡ n | (W|n n|) |n . And since we focus on low temperature (compared to ω 0 ), we find that ∆ ≈ ζ. Thus, we have three independent correlation functions determining the kernel, which, for realistic predictions, would have to be identified experimentally along similar lines as Ref. [84]. Here, for demonstration purposes (to show that the above discussed correction does not vanish in general), we assume the reservoir to be an ohmic bath with a cutoff frequency ω c ω 0 . Herewith, one finds the following scalings of the correlation functions [139]  α,d , both relative to the current I ind = −2eṄg, induced by the ramping of the offset charge. Those numerical results are obtained by assuming a large junction asymmetry, EJR/EJL 1 or EJR/EJL 1. In the intermediate regime, the perturbation cannot be expected to remain small (or even finite) due to the energy gap decreasing to zero when approaching the point degeneracy. For a large asymmetry, the deviation is quadratically suppressed. We chose the parameters EJR+EJL = 50EC , ωc = 100EC , and eitherφ = 0.1Ṅg oṙ φ = 0.5Ṅg to demonstrate that the current correction scales withφ = 2eV , and not withṄg as for the closed-system current.
Assuming a junction asymmetry described by ξ 1 allows us now to expand the current correction from Eq. (23) w.r.t. ξ.
Since all φ-dependencies enter via a dependence on E J (φ) or ∂ φ δ(φ) = E JR (E JR + E JL cos φ) /E 2 J (φ), we find that the zeroth order (in ξ) is always a constant term (in φ) while the first order is associated with a harmonic dependence. Thus, odd orders of the current vanish when averaging over φ (corresponding to a time integral), and since ∂ φα n has a vanishing zeroth order, the correction has to be at least of quadratic order in the asymmetry ξ. This can be confirmed numerically: see Fig. 3, where we explicitly expand Eq. (23) up to second order in ξ, and subsequently perform an integration over φ to obtain the dc component of I We conclude that while the dc part of the open system correction to the current does not vanish, and gives rise to a deviation from the otherwise perfect current quantization, our analysis offers very concrete indications how to minimize its influence. Namely, by means of our general discussion above, we identify an important difference in the scaling behavior between the topological part of the current response, and the open system correction: while the former appears with a prefactor proportional toṄ g [see Eqs. (12) and (13)], the correction in Eq. (23) turns out to scale withφ = 2eV , such that for this deviation to be small, V should not be chosen too large with respect to the ramping of N g . To conclude, we find that the correction is small as long as ξ 2 V /Ṅ g < L 0 / W . Figure 4. Measurement protocol. a) The right Josephson junctions is replaced by a SQUID whose energy EJR (φext) can be tuned by the magnetic flux φext. The voltage V drives the phase difference continuously, while the linearly timedependent gate voltage Vg (t) induces a dc current I ind ∝Vg. b) A suggestion on how to tune Ng via Vg and EJR via φext as a function of time in order to be able to measure a quantized dc current. Whenever we ramp up Ng, the induced current flows into the right lead and when we ramp it back down, the current flows from the left lead into the system. In the intermediate steps, Ng is held constant to adjust φext such that EJR becomes larger or smaller than EJL.
Therefore, while the open system correction is not exponentially suppressed (which seems to be a generic feature of open systems; see, e.g., a recent discussion for topological insulators [140]), our above calculation provides very clear and stringent strategies to mitigate it, by choosing the driving and other parameters accordingly. In particular, as we explicitly show in Fig. 3, the overall opensystem correction can be kept very small, due to ξ 1 andφ <Ṅ g . In addition, given the level of generality of Eqs. (23) and (24), we expect that a similar perturbation occurs in the models considered in Refs. [77,78].

V. DC CURRENT MEASUREMENT
As we have indicated in Sec. III, a remaining experimental obstacle concerns the fact that N g cannot be ramped up indefinitely since at some point the transistor will break. This limitation can easily be circumvented with a simple procedure using the fact that the direction of the quantized dc-current is sensitive to the junction asymmetry. Namely, we replace the right junction with a superconducting quantum interference device (SQUID) consisting of two parallel junctions, each with an energy E JS > E JL /2 (see Fig. 4a). This introduces a tunable Josephson energy E JR → E JR (φ ext ) = 2E JS cos (φ ext /2), controlled by an external magnetic flux going through the SQUID [141].
In a first step, the protocol now simply consists of ramping up N g to a maximal value, related to the maximum gate voltage V max , in one configuration, e.g., where E JR (φ ext ) > E JL (such that the induced dc current flows into the right contact). Afterwards φ ext is changed to a value where E JR (φ ext ) < E JL , and N g is subsequently ramped back down to a minimal value related to the gate voltage V min (pumping the dc current into the system from the left contact), while keeping the bias voltage V on for all times (see Fig. 4b). The time required to switch the junction asymmetry is referred to as t switch , while a single ramping goes on for t ramp . For long times we will measure the averaged dc current which thus depends only on the single driving parame-terṄ g and the two relevant times of the cycle, which are completely controlled by the experimenter. In the limit of t switch t ramp , we find I dc = eṄ g . Note that in a single ramping process the current 2eṄ g flows, as per Eq. (12) and Eq. (13). However, we need two individual ramping processes (ramping the offset charge in both directions) to complete the cycle, which takes twice the time. Furthermore, we stress that as long as a transition from E JR (φ ext ) < E JL to E JR (φ ext ) > E JL can be achieved, the flux control does not even need to be very precise, nor is it susceptible to flux noise, apart from the above discussed finite-frequency perturbations. Note that the only important restriction is to make sure that V min and V max are chosen such that the system does not go through the degeneracy point (Weyl point) when ramping from E JR < E JL to E JR > E JL (and vice versa). However, we expect even such a dissipationinduced deviation will likely not be dramatic due to this event not being very probable and giving only a small contribution compared to the topological part as long as the difference V max − V min is chosen sufficiently large.
In fact, this protocol has a lot of similarity with Cooper-pair sluices [110][111][112][113]. However, one advantage of our approach is that we do not require a precise control of the tunnel couplings to the contacts: the quantization of the current requires merely the averaging in the (N g , φ)-space, which is guaranteed in the presence of the bias voltage V . Moreover, contrary to regular Cooperpair pumps [113], our proposal is insensitive to fermion parity, as we argued above.

VI. CONCLUSION
We have found that the Cooper-pair transistor hosts topologically nontrivial Chern numbers, giving rise to a quantization of the dc current response, which is precisely steered either to the left or right contact. This circuit has various advantages to alternative systems, not least the simplicity and straightforward realizability of the circuit. Surprisingly, low-frequency charge noise is actually beneficial for the observation of the quantization effect. Moreover, the Chern number is insensitive to quasiparticle poisoning and to whether or not the system is in its ground state. The latter is due to the emergence of Weyl points with higher topological charges connecting higher bands. Remaining environment-introduced perturbations are found to be small. Finally, we presented an experimentally feasible protocol to carry out the dc current measurement. We conclude that the Cooperpair transistor presents a promising platform to realize a topological circuit, with a topological number defined in a "mixed" basis consisting of the phase difference φ and the offset charge N g .
Finally, when considering topological systems, the question of the existence of protected edge states inevitably occurs. Usually, in materials with a topological band structure in k-space, edge states naturally occur at the termination of the material (in position space). Since we here consider topological numbers in an alternative base space, the physical presence of edge states becomes somewhat elusive. Nonetheless, they are present in the following sense. A sharp boundary in x-space corresponds to a highly nonlocal feature in (canonically conjugate) k-space as per the Heisenberg uncertainty principle, such that the edge is able to "probe" the bulk topological number directly [142]. In our proposal on the other hand, the Chern number is probed via time-dependent driving, and subsequent time-averaging of the electric current response. Since the current corresponds to the total number of transported particles (per time), and this particle number being conjugate to φ, one can interpret the quantization of the former as an indirect probe of an edge state. The explicit creation of edge states in charge space (e.g., via an active shaping of the charge space itself) will be a crucial future research endeavor, in order to use the above topological features for protected quantum information processing. When tuning the Hamiltonian [see Eq. (1)] near the degeneracy, that is N g = −N + 1/2 + δN g , E JR = E JL + δE J , and φ R = π + δφ while at the same time φ L = 0, we find up to first order in δN g , δE J , and δφ Note that we can ignore the first term, because, within the subspace spanned by the states |N and |N − 1 , it can be regarded as a constant energy contribution.

Double Weyl point
Now, we derive the Hamiltonian near the crossing point between bands n = 1 and n = 2, describing a Weyl point with topological charge C = +2, as in Eq. (6).
As pointed out in the main text, we have to tune to N g = −N +δN g , E JR = E JL +δE J , and φ R −φ L = π+δφ to get close to the double Weyl point. Here, a gapping can occur only by changing between charge states |N − 1 and |N + 1 , which is achieved by a higherorder process involving the tunneling of two Cooper-pairs via virtual charge states. We tackle this problem by means of a Schrieffer-Wolff transformation. First, we are shifting the reference point of the energy by the average value of the charge subspace {|N − 1 , |N + 1 }, H → H− 1 2 tr H P , with the projector onto this subspace P = |N − 1 N − 1| + |N + 1 N + 1|. Afterwards, we can write the effective Hamiltonian in the low-energy regime approximately as when decomposing the Hamiltonian according to H = H 0 + V , with and keeping only the lowest order in δN g , δφ, and δE J . Inserting those two into Eq. (A3), we find and where we defined N ≡ E C /2 (N − N ) 2 − 1 and v ≡ (δE J + iE JL δφ)/2. We ignore the first two terms of Eq. (A7) since they give us only what can be regarded as constant energy contributions. Inserting Eqs. (A6) and (A7) back into Eq. (A3), we find arriving at Eq. (6).

Appendix B: Open-system correction terms
Here, we apply perturbation theory to the master equation, Eq. (15), to derive the corrections to the current expectation value as shown in Eq. (17), arising through a coupling to the environment and a time-dependent driving of the system. We vectorize ρ and decompose it into diagonal and off-diagonal sectors (in the eigenbasis of H) |ρ) = (|ρ d ) , |ρ o )), each of which are represented by vectors in the Fock-Liouville space In the same manner, we write the superoperators in the corresponding matrix representation with four subblocks. The two diagonal subblocks of these matrices correspond to transitions from diagonal to diagonal, respectively offdiagonal to off-diagonal sectors. The off-diagonal matrix subblocks describe the coupling between the diagonal and off-diagonal sectors, such that the master equation has the following shape where L oo denotes the nonzero part of the Liouvillian superoperator L 0 . We demand that W dd has one (nondegenerate) zero eigenvalue (equivalent to it having a stationary state) W dd ρ (0,0) d = 0. Let us additionally suppose that W dd relaxes the system to the ground state (up to exponentially suppressed contributions, equivalent to a low temperature assumption) such that ρ (0,0) d = |0 d ) (the vectorized version of ρ (0,0) = |0 0| appearing in the main text as ρ 0 ). The subscript d simply expresses the fact that here, the size of the vector is that of the diagonal sector, such that Be aware that this state is not the exact stationary state of the total kernel W because we are explicitly allowing transitions between diagonal and off-diagonal states in the kernel, represented by the off-diagonal subblocks.
Corrections to ρ (0,0) are deviations from the commonly assumed rotating-wave approximation (RWA). Therefore, in the absence of the time-dependent driving, we will always end up in a stationary state of the shape ρ st = ρ (0,0) + ρ (0,1) , where ρ (0,1) denotes the first-order correction beyond RWA. The notation (0, 1) refers to this contribution being of zero order in the driving δL and of linear order in the kernel W ; see also Eq. (16). However, the stationary state, ρ st , does not contribute to the dc part of the current. As we already argued in the main text, this is due to the fact that dc currents can be induced only by driving the system out of equilibrium and therefore cannot be connected to terms of zero order in the driving with the environment at equilibrium [138]. Turning on the slow driving of the system now gives us a small correction to the density matrix, driving the system away from the stationary state to what we call a quasistationary state, ρ. Assuming that δL < W < L 0 , we can expand this state in both W and δL , analogously to Eq. (16), as |ρ) = ν,µ ρ (ν,µ) .
We now define the quasistationary state ρ such that the time derivative of each order ∂ t ρ (ν,µ) is of higher order in the driving (ν +1, µ). Here, the zeroth order, ρ (0,0) = |0), is constant, meaning it only implicitly depends on time due to the now time-dependent basis (which is the instantaneous eigenbasis of H).
Besides the previously mentioned ρ (0,1) , we now find the first-order correction ρ (1,0) which has a nonzero dc contribution and therefore gives rise to the leading-order terms of the dc current. Even though ρ (0,1) does not contribute to the dc current, we nonetheless have to compute its off-diagonal part ρ Reexpressing Eq. (B2) for the diagonal and offdiagonal sectors separately to find the quasistationary state, we find We begin by considering the stationary open system (i.e., in absence of the driving) equivalent to keeping only terms of zero order in the driving. In leading order of W , we find providing the stationary state with Here, W −1 dd is defined such that with I dd being the identity matrix in the diagonal subspace and (0 d | being the left eigenvector of W dd with eigenvalue zero, which is the trace (0| • = tr[•]. We use this notation to denote a map from an operator to a scalar, such that the scalar product (0|0) = tr ρ (0,0) = 1. We can derive the corrections for the quasistationary state in the presence of the drive by comparing the terms of Eqs. (B5) and (B6) that are of linear order in δL , where we find in leading order of W ∂ t ρ (B13) This second equation is fulfilled for We stress that this additional correction arises solely because of the nonzero off-diagonal subblock in the kernel, W do , which would vanish when making the RWA. Strikingly, it is thus an open-system correction of zero-order in W . This contribution, Eq. (B15), gives rise to the leading-order open-system correction to the current; see Eq. (23).
For a system-reservoir interaction that conserves the current (see discussion in the main text), we can write the current expectation value as I α = tr I α ρ = (I α |ρ), introducing the notation (A| • = tr A • . I α can be computed with the help of the operator N α = −i n,m n| ∂ φα |m |n m|, as already introduced in Eq. (22). It can easily be shown that from which directly follows I α = 2e n ρ nn ∂ φα n + 2e (N α | iL 0 |ρ) , which we have written in the form I α = I α,d + I α,o , analogously to the current in the closed system in Eq. (11).
Here, I α,d (I α,o ) represents the current contribution due to the diagonal (off-diagonal) elements of the density matrix. The diagonal term, I α,d = n ρ nn I (0) α,n , represents a generalized version of the zero-order contribution corresponding to the ordinary Josephson effect of the closed system, I (0) α,n = 2e ∂ φα n [see Eq. (11)]. However, when inserting the order-by-order expansion of the density matrix, Eq. (B4), we find that in contrast to the closed-system result, this expression contributes to the dc current, in leading-order via I α,0 = −2e (N α | iδL |0); see Eq. (11). The linear open-system correction I (0,1) α,o = 2e (N α | W |0) is again purely contributing to the ac current, like any higher-order term that does not depend on the driving parameters. We thus find the dc current response I α = I