# Superadiabatic theory for Cooper pair pumping under decoherence

###### Abstract

We introduce a method where successive coordinate transformations are applied to decrease the error in the adiabatic master equation resulting from truncation in the local adiabatic parameter. Our method reduces the non-physical behaviour stemming from the lack of complete positivity. The strong environment-induced relaxation at high Cooper pair pumping frequencies leads to adiabatic ground-state pumping only in the lowest-order approximation. We illustrate the robustness of the frequency where the adiabaticity breaks down using the high-order theory and show the emergence of an optimal environmental coupling strength, for which ideal pumping is preserved for the highest frequency. Finally, we study the effect of quantum interference on the pumped current and give an estimate for the relaxation rate of an experimentally measured system.

###### pacs:

03.65.Yz, 03.65.Vf, 85.25.Cp, 85.25.Dq## I Introduction

Detection and manipulation of geometric phases prsla392/45 (); prl51/2167 (); prl52/2111 (); prl58/1593 (); prl60/2339 () in superconducting quantum devices has been an area of active research in recent years with one of the ultimate goals being the ability to realize holonomic quantum gates pla264/94 (). Even though alternative methods to experimentally generate and detect the geometric phases in such devices have been proposed and realized nature407/355 (); science318/1889 (), the link they have to Cooper pair pumping has attracted major interest prb60/R9931 (); prb63/132508 (); prb68/020502(R) (); prl91/177003 (); prl95/256801 (); apl90/082102 (); prl98/127001 (); prl100/027002 (); prl100/117001 (); prb77/144522 () as it provides means of detecting the phases by a measurement of the dynamically and geometrically transferred charges. A measurement scheme based on one superconducting island with two tunable Josephson junctions, the Cooper pair sluice in a superconducting loop, has been introduced prb73/214523 () and experimentally realized prl100/177201 (). Recently, a system based on a similar device structure has been proposed to execute fully geometric quantum computing using non-Abelian phases prb81/174506 (); pra82/052304 ().

Even though the principles of operation for Cooper pair pumps are well-known, accounting for system–environment interactions has been a work in progress. Cooper pair pumping being essentially a coherent process, a proper description of the operation must include the effect of the external environment. The usual methods for describing the dynamics of open steered quantum systems tToOQS (); API (); prl94/070407 () have been shown to violate the conservation of physical observables such as the electric charge and a variety of means have been employed in an attempt to properly describe the dynamics pra73/052304 (); pra73/022327 (); prl90/160402 (); prb77/115322 (); pra82/052107 (); pra83/012112 (); pra81/022117 (); pra83/032122 (). Recently, it was discovered that for a consistent description using the master equation for the reduced density matrix of the system, all the non-secular terms must be included to enable relaxation to a proper basis and to ensure conservation of the pumped charge in Cooper pair pumping prl105/030401 (); prb82/134517 (). The same master equation has also been derived using superadiabatic bases pra82/062112 ().

In this paper, we employ the recent methodology of superadiabatic bases pra82/062112 () in deriving the master equation for steered systems and apply it after multiple coordinate transformations of the time-local basis. This results in a master equation where the truncation error related to the local adiabatic parameter is potentially decreased as a function of the number of coordinate transformations. We apply our method to the problem of Cooper pair pumping. We show that for the zero-temperature environment and fast pumping, the ground-state adiabatic evolution revives in the relaxation dominated region only using the adiabatic basis, that is, in the lowest order of our description. Furthermore, we show that the overestimation of the pumped charge caused by the non-positivity of the density matrix, is alleviated dramatically by the usage of high-order bases. We simulate the breakdown of adiabaticity with increasing pumping frequency and show the emergence of an optimal coupling strength preserving ideal pumping up to the highest frequency. We present a condition for the highest transition probability caused by constructive interference between driving-induced excitations generated at different times and show that it corresponds to the downward resonance peaks in the pumped current. Finally, we obtain an estimate for the relaxation rate of the device employed in the experiments of Ref. prl100/177201, to pump Cooper pairs.

The structure of this paper is as follows. In the next section, we introduce the model describing a driven quantum system and demonstrate our method of defining successive effective Hamiltonians by coordinate transformations. In Sec. III, we write a master equation for the matrix elements of the reduced density matrix of the system taken in an -times transformed time-dependent basis. In Sec. IV, we use the master equation to model Cooper pair pumping. Furthermore, we simulate previous experiments on the pumped current and derive an estimate for the relaxation rate of a measured superconducting system. We conclude the paper in Sec. V.

## Ii Model

We study a quantum system with a Hamiltonian which depends on a set of real control parameters that vary in time. The system is assumed to be interacting with its environment such that the total Hamiltonian is

(1) |

where is the coupling between the system and its environment and is the Hamiltonian of the environment. We assume that the coupling is of the generic form , where is the system part of the coupling operator and acts in the Hilbert space of the environment. Let be the instantaneous eigenstate of and the corresponding eigenenergy defined by . In the context of adiabatic evolution, is referred to as the adiabatic basis. We assume that the adiabatic states are normalized and non-degenerate.

Let the Hamiltonian be diagonalized in a fixed time-independent basis using the eigendecomposition as , implying that . We define a similar transformation for the total density operator in the Schrödinger picture as . It follows from the Schrödinger equation that the evolution of is governed by the effective Hamiltonian for the adiabatic basis

(2) |

where and . The eigenbasis of is usually referred to as the superadiabatic basis.

We can further define a unitary transformation making diagonal in the fixed basis pra82/062112 (). Thus the evolution of the density matrix is governed by the effective Hamiltonian for the first superadiabatic basis

(3) |

where , , and . This method of successive coordinate transformations can be continued to yield for the th superadiabatic basis an effective Hamiltonian of

(4) |

where , and , where we omitted explicitly marking the temporal dependence of the operators for clarity. The operator product is defined as . If we define for and for , the density operators governed by the Hamiltonians in Eqs. (2), (3) and (4) obtain a more universal form . Defining successive diagonalizations in this manner proves useful in Sec. III as the recently derived master equation prl105/030401 (); prb82/134517 (); pra82/062112 () can be applied to solve the system dynamics using these high-order bases.

The iterative method described here is an adaptation of Berry’s concept prsla414/31 () he later referred to as adiabatic renormalization GPIPberry (). It is based on the idea that each transformation rotates the basis we use to describe the system dynamics ever closer to the exact evolving closed system state, that is, the time-dependence of the transformed system Hamiltonian is suppressed after each rotation. After transformations, we define the time-dependent basis as . This approach generally works only in the restricted sense, that is, after a number of iterations, the following rotations will not allow us to describe the dynamics of the system more accurately prsla414/31 ().

Finally, we introduce the local adiabatic parameter as , where we compare the Hilbert-Schmidt norm of the operator arising from the adiabatic evolution to an instantaneous minimum energy gap in the spectrum . Here denotes the trace over the system degrees of freedom and in the following we will use to denote the trace over the environment degrees of freedom. The parameter should give a good estimate for the degree of adiabaticity of the evolution prb82/134517 (); pra82/062112 (). In cyclic evolution with the period , the parameter scales as and, thus, in adiabatic evolution we should have .

## Iii Master equation

We consider an adiabatically steered two-level quantum system weakly coupled to its environment. We denote the ground and excited states of in the Schrödinger picture as and , respectively, with corresponding eigenenergies and . Using the interaction picture approach, a master equation was derived to describe the dynamics of such a system in Refs. [prl105/030401, ; prb82/134517, ; pra82/062112, ] up to the linear order in and the quadratic order in the system–environment coupling. The method of derivation employed in Ref. pra82/062112, is our starting point for developing a numerical scheme for obtaining a more accurate description of the dissipative system dynamics. In Ref. pra82/062112, , a master equation for nonsteered systems was used in conjunction with the effective Hamiltonian in Eq. (3) to derive the leading-order master equation under steering. However, a similar derivation can be carried out using for any to obtain a master equation for the matrix elements of . Notice that even though the method of defining successive coordinate transformations can be applied to a system with arbitrary number of energy levels, we constrain ourselves to the two-level case. This is practical since our main goal is to explore the implications of applying our scheme compared to previous results prl105/030401 (); prb82/134517 ().

We define the reduced density operator of the system as so that its diagonal element becomes and the off-diagonal element , where is the relevant fixed basis. These are simply the matrix elements of the usual density operator of the system in the Schrödinger picture taken in a time-dependent basis , where and . This is the rotated basis obtained through the iterative procedure we described in Sec. II. We emphasize that the basis states are not obtained using a perturbative expansion in the local adiabatic parameter and, thus, each iteration generally alters them by terms of all orders of prsla414/31 (); pra78/052508 (). However, if the method described in Ref. pra82/062112, is applied in the -times transformed basis, the error in the resulting master equation is defined by a perturbative expansion in . The norm of decreases, in the restricted sense, with increasing as the time-dependence of the transformed system Hamiltonian is suppressed. In regard to depicting the actual dynamics of the system, the only issue relevant to the selection of the basis is that the time-evolution of the density matrix elements can be accurately described using it. Thus, we can exploit the th iterative basis and define a master equation up to the quadratic order in the system–environment coupling and to the first order in , where such that and , as

(5) |

and

(6) |

Furthermore, we denote , , , and . The reduced spectral density of the noise source is defined as . Similarly to Refs. [prl105/030401, ; prb82/134517, ; pra82/062112, ], we assume that the system is in the Markov regime, the system time scales are longer than the environment autocorrelation time leading to neglecting the Lamb shift, and the approximation of adiabatic rates applies. These assumptions and the time scale separation they lead to are described in detail in Ref. pra82/062112, .

As described above, the benefit of defining the successive coordinate transformations of the Hamiltonian is that the corresponding master equation is up to the first order in , thus describing the evolution of the system more accurately, in the restricted sense, as increases. Defining a master equation of arbitrary order in using the original methods prl105/030401 (); prb82/134517 (); pra82/062112 () is possible, but the effort required renders such derivations highly unpractical.

It has been shown prl105/030401 (); prb82/134517 () that assuming a zero-temperature environment and taking the quasi-stationary limit, the master equation in the lowest order leads to and . This translates to and showing that in the first order in , the density matrix describes the evolution of a pure state. This is a remarkable result validating that the master equation in Refs. [prl105/030401, ; prb82/134517, ; pra82/062112, ] ensures relaxation to up to the first order in . Similarly, using our master equation for ensures that the relaxation takes the system to up to the first order in . Especially, in the limit , the rotational terms , , in the master equation vanish and the basis fully describes the steering assuming that the process of basis rotations converges. The requirements for the convergence or the number of transformations up to which the iterative procedure suppresses the time-dependence of the Hamiltonian of the system prsla392/45 () are not studied in this paper, as it turns out in Sec. IV that a small number of transformations allows one to capture the main effect of this scheme.

## Iv Cooper pair sluice

### iv.1 Definitions

We introduce the Cooper pair sluice prl91/177003 () as a physical realization of a steered two-level system. The charge pumped through the sluice establishes a connection to geometric phases prb60/R9931 (); prb68/020502(R) (); prb73/214523 (); prl100/117001 (); prl100/177201 () acquired during the adiabatic evolution and provides a physical observable. We aim to study and improve on the recent theoretical pumping results prl105/030401 (); prb82/134517 () using the high-order effective theory.

The Cooper pair sluice shown in Fig. 1(a) is comprised of a superconducting island separated by two SQUIDs ItS (), each involving two Josephson junctions.

If we assume that the self-inductances of the SQUID loops are negligible, the SQUIDs operate as tunable Josephson junctions whose Josephson energies , where , are determined by the external fluxes threading the loops. The sluice Hamiltonian is

(7) |

where the Coulomb energy for one excess Cooper pair is and the gate charge is defined in units of as . Here is the gate capacitance and stands for the total capacitance of the island. The operator describing the phase on the island is and its canonical conjugate, the operator describing the number of excess Cooper pairs on the island, is . The gauge-invariant phase difference over the device is determined by and kept constant during the evolution.

We denote () as the Cooper pair number operator of the th SQUID and write the average value of the current through the th SQUID as prb82/134517 ()

(8) |

where denotes the respective current operator. If the environment does not directly induce a current, that is, , Eq. (8) defines the usual current operator prb73/214523 ()

(9) |

If we set and , the dynamics are accurately described by the two lowest charge states allowing us to apply the preceding two-state theory. We denote and as the states with no and with one excess Cooper pair on the island defining our fixed basis. We study the capacitive coupling of the environment to the system by introducing voltage fluctuations at the gate of the sluice prl105/030401 (); prb82/134517 (). The coupling operator then becomes where and denotes the strength of the coupling. The coupling operator has been selected traceless in the two-state basis by adding an operator comparable to the identity operator to adopt a convention used in the derivation of the master equation pra82/062112 (). Such a selection can be applied to any coupling operator and it does not reduce the generality of the master equation. Since , Eqs. (7) and (9) imply that .

Assume that the noise source is a resistor in thermal equilibrium; a situation which can be engineered in the physical realization of the sluice. If we consider the voltage noise of a resistor grounded at one end and connected to the gate by a low impedance circuit at the other end, the reduced spectral density of the noise source at the gate becomes rmp82/1155 () , where is the effective resistance of the noise source and is the resistor temperature. We note that the detailed balance condition applies. Furthermore, we introduce dephasing to the system by assuming that , where is the effective dephasing temperature.

We denote the matrix elements of the current operator of the th SQUID by . The expectation value of the current using the adiabatic basis is

(10) |

since is Hermitian. The first two terms are the dynamic supercurrent through the junction and the third term describes the geometric part of the current prb73/214523 (); prl100/177201 (); prl105/030401 (); prb82/134517 (). The pumped charge corresponding to the geometric contribution becomes

(11) |

where is the length of the closed cycle.

The definition for the different terms in the average current only applies directly for the adiabatic basis. This implies that if we pursue to obtain the geometric current using higher order bases, should be written using the density matrix elements in the basis where the evolution takes place in our calculations. The adiabatic density matrix element can be written as

(12) |

Using this, we can rewrite the integrand in the pumped charge

(13) |

Equation (13) enables us to calculate the pumped charge when the time-evolution of the th basis density matrix is known. Note that for the adiabatic basis, Eq. (13) reduces to the form corresponding to Eq. (11).

### iv.2 Effect of the environment on the pumped charge

We use the parameter cycle shown in Fig. 1(b) for the pumping and ensure the smoothness of the parameter functions in time using trigonometric interpolation dividing the total cycle time into 201 equidistant points. This is a necessary step for the simulations as exploiting the high-order bases requires that the high-order temporal derivatives of the parameter functions stemming from are non-divergent. The dynamics of the quantum system are solved numerically from Eqs. (5) and (6) utilizing the effective Hamiltonians introduced in Sec. II. The density matrix and any physical observables are recorded in the steady state, that is, after sufficiently many cycles such that the system evolution of consecutive cycles is identical.

We begin by studying the effect of using the higher-order bases when describing the dynamics of the sluice in a zero-temperature environment. Using the lowest-order approximation , the analytical result of the environment inducing ground-state evolution in the adiabatic limit has been demonstrated using numerical simulations of the pumped charge prl105/030401 (); prb82/134517 (). Additionally, increasing the pumping frequency has been shown to induce regions where either the nonadiabatic transitions or relaxation dominates depending on the ratio . However, it turns out that a region in the -space emerges where the pumped charge is unphysically overestimated. This is a direct result of the master equation not strictly ensuring the positivity of the density matrix in any finite order, that is, it does not reduce to the standard Lindblad form.

We present the pumped charge using the basis with as a function of the coupling strength for different pumping frequencies in Fig. 2(a) assuming a zero-temperature environment. We explore the regime where the strength of the environmental coupling is small to ensure that we remain close to the weak coupling limit.

In the adiabatic region ( MHz), all orders of approximation indicate ground-state pumping for all environmental coupling strengths as predicted in Sec. III. By increasing the pumping frequency, we observe the emergence of the two pumping regimes mentioned above. Additionally, Figure 2(a) illustrates how increasing the coupling strength does not lead to ground-state pumping beyond the adiabatic region for . The reason for this phenomenon stems from the structure of the superadiabatic bases. For nonadiabatic evolution, increasing the coupling strength leads to relaxation to up to the first order in which translates to ideal pumping only for as the system is forced to the solution of the adiabatic limit, from which the asymptotic solutions for generally deviate in all orders of the local adiabatic parameter. The adiabatic solution for the pumped charge assuming similar SQUIDS has been derived previously prb73/214523 () as . To make the discrepancy between the high-order bases and the case more visible, we present the pumped charge up to high coupling strengths in Fig. 2(b). The regime of increased coupling strength is beyond the range of validity of our approach but shows how the master equation properly displays the distinctions between the different bases as they approach the asymptotic relaxation dominated solutions.

From Fig. 2, it is clear that the charge overestimation observed for is alleviated by utilizing the high-order bases. As the high-order bases follow the exact evolving closed system state more closely, the non-adiabatic transitions disturb the mixed state less. The major contribution in alleviating the lack of complete positivity is already given by and the subsequent third basis rotation has little effect on the pumped charge in comparison. The effect of using the higher-order bases is evident from studying not only a set of observables but from the density matrix itself. We present the lower eigenvalue of the steady state density operator in Fig. 3.

The non-positivity is greatly reduced and for some parameter values, completely removed by doing one or more further rotations beyond the adiabatic basis.

### iv.3 Pumped current and the breakdown of adiabaticity

So far, experimental results for Cooper pair pumping using the sluice have been scarce prb71/012513 (); apl90/082102 (); prl100/177201 (). However, the breakdown of adiabaticity has been observed by studying the pumped current as a function of the pumping amplitude , that is, at high pumping amplitudes the pumped current has been noticed to deviate from the analytical result in the adiabatic limit , where is the number of Cooper pairs transported ideally through the sluice per cycle prl100/177201 (). The number of transported Cooper pairs can be experimentally dictated by adjusting the gate voltage, more spesifically, by altering so that the ideally transported average current corresponds to . Unfortunately, our master equation cannot be directly used to simulate the effect of altering the pumping amplitude since it is defined in a two-state basis which requires that remains approximately half-integer during the evolution. Any selection of the two charge states as the fixed basis enables the maximum pumping amplitude of one Cooper pair per cycle.

Even though we cannot modify the pumping amplitude, we can still simulate the breakdown of the adiabaticity by altering the pumping frequency. If we assume that the deviation from the adiabatic behaviour in the experiments is due to increase in the pumping speed caused by the amplitude growth using a constant cycle time, the effect of decreasing the total cycle time should be similar. We choose a frequency range beyond the strict adiabatic limit to model experiments carried out with finite cycle times. The pumped current is shown in Fig. 4(a) using and in Fig. 4(b) using .

The behaviour of the pumped current with high frequencies is only suggestive but clearly indicates that the effect of the number of coordinate transformations increases with the pumping frequency. The physically most relevant features are found near the point where the adiabaticity breaks down. As we utilize a more accurate description of the dynamics, the point where the adiabaticity of the system is broken becomes more robust against changes in the environment. Furthermore, we observe the emergence of an optimal coupling strength, with which the ideal ground-state pumping is conserved up to the highest frequency. This feature could have also been anticipated from Fig. 2 where we observe the emergence of a maximum pumped charge as a function of the coupling srength for any given frequency far from the adiabatic limit if . This corresponds to the coupling strength, below which the non-adiabatic transitions reduce the pumped charge and above which the relaxation takes the system away from the solution in the adiabatic limit. The optimal coupling strength can be probed and exploited utilizing the environment engineering scheme presented in Ref. prb82/134517, .

We turn our attention to modeling the breakdown characteristics of the actual experimentally pumped current in Ref. prl100/177201, . The simulation is motivated by the above discussion on the pumping speed and we establish an equality between the experimental pumping speed , where is maintained at a constant value, and the pumping speed in our simulation , where is a constant to ensure that the two-state approximation holds. A comparison between the experimental results and our simulations is presented in Fig. 5.

The system parameters used in the simulations are estimated from the experiments and can contain up to 20% error. To simulate the experiments, the aforementioned equality defines the needed frequency by MHz so that the point where the adiabaticity breaks down according to the experimental data would require a simulation frequency of approximately 1.1 GHz. This is significantly higher than the breakdown frequency in our previous simulations and tests the limits of validity of our approach.

Assuming both vanishing dephasing in Fig. 5(a) and non-vanishing dephasing in Fig. 5(b), the simulated pumped current exhibits strong oscillatory behavior allowing for weak predictability of the exact breakdown characteristics. The observed behavior is due to quantum interference between driving-induced excitations generated at different times. The instantaneous transition probability is not only dependent on the energy gap but also on the phases accumulated during the quantum evolution. Especially, if the phase difference between two successive excitations in time is a multiple of , the transition probability reaches its maximum due to constructive interference and the geometric current obtains a downward resonance peak. Studying the phase accumulation during the quantum evolution allows us to estimate the resonance peak positions and compare them with the simulated pumped current.

The time-evolution of the energy gap of the system is symmetric with respect to the mid-point of each pumping cycle corresponding to times , . In addition, the energy gap decreases during the two gate manipulations and as both manipulations in Figs. 5(a) and 5(b) are symmetric with respect to the degeneracy point , the energy gap reaches its minima at points we denote as . Hence, it is enough to study the phase accumulation between different to depict the resonance behavior.

In adiabatic evolution, the number of coordinate transformations we perform has a profound effect on the observed accumulated phase, that is, each transformation takes us to a new basis in which the phase accumulates differently. More exactly, the closed system state follows the Schrödinger equation at all times. After transformations, the state follows the transformed equation . Assuming adiabaticity in this basis, that is, the exact evolving state remains in the th eigenstate of , the -times transformed evolving state can be written, similarly to Ref. ajp66/431, , as

(14) |

where we assumed that , , and describes the geometric phase contribution. Using the familiar notation, this implies and . The accumulated quantum phases can be obtained from Eq. (14) following the derivation in Ref. ajp66/431, . However, in our dissipative calculations, after transformations, we take as the perturbation. This means that for the dissipative simulations, the relevant reference frame after transformations is given by the eigenbasis of for and for , that is, and .

Concentrating on the two-state model, we present the accumulated phases after transformations. We assume that the transitions taking place at are instantaneous and the system evolves adiabatically between them. The dynamically accumulated phase obtained during adiabatic evolution between and for the eigenstates in the relevant reference frame after transformations is given by , , and it is equal for any successive two points due to symmetry. The geometrically accumulated phase