Pumping current in a non-Markovian $N$-state model

A periodically modulated N-state model whose dynamics is governed by a time-convoluted generalized master equation is theoretically analyzed. It is shown that this non-Markovian master equation can be converted to a Markovian master equation having a larger transition matrix, which affords easier analysis. The behavior of this model is investigated by focusing on the cycle-averaged pumping current. In the adiabatic limit, the geometrical current is calculated analytically, and compared to numerical results which are available for a wide range of modulation frequencies.


I. INTRODUCTION
Master equations (MEs) are widely used in nonequilibrium statistical mechanics to model the time evolution of a range of classical and quantum mechanical systems. The mathematical foundation of MEs is the differential Chapman-Kolmogorov equation of stochastic analysis [1], and in its basic form, it only describes systems that carry no memory of their past: this property is referred to as Markovianity. Is is known, however, that due to a number of physical reasons, real systems do usually possess some degree of memory of their past evolution, and thus obey a non-Markovian ME (nMME). An archetypal example of this is the fluctuation-dissipation relation [2], which shows that any time-nonlocal correlations in the environment necessarily lead to memory effects. Theoretically, non-Markovianity has been attracting attention in classical mechanics as a link between continuous-time random walks and time convoluted MEs [3,4]. In quantum physics, the connection between the flow of information and non-Markovian processes, as well as quantum measurements of Markovianity have received considerable research effort [5][6][7]. Furthermore, advancements in experimental techniques have made it possible to directly measure non-Markovianity in the context of classical [8] and quantum [9,10] systems.
Modulating control parameters such as rate constants, bath temperatures or gate voltages of a physical system out of equilibrium can lead to net flow of a physical current, e.g. flow of product, heat flow or electron current, even in the absence of a net bias of the control parameters [11]. The origin of the pumping current is the Berry-Sinitsyn-Nemenman (BSN) phase, which was originally discovered in the context of quantum systems [12][13][14]. Geometrical current is generated when the modulation speed of the parameters is sufficiently slow, i.e., in the adiabatic limit. The BSN phase has a significant impact on quantum transport [15,16], and it also leads to a path-dependent geometrical entropy [17,18]. The effects of the BSN phase have been extensively studied in the spin-boson model [19]. Over the recent years, effects of * E-mail: ville.paasonen@yukawa.kyoto-u.ac.jp finite speed modulation, i.e., non-adiabatic effects, have also been investigated for this model [20]. Furthermore, it has been shown that the presence of the BSN phase engenders non-Gaussianity of the system fluctuations, leading to a modified form of the fluctuation theorem for geometrical pumping [21,22]. Moreover, a heat engine using Thouless pumping [11] associated with the BSN curvature has been theoretically designed [23].
The essential features of this geometrical effect have also been shown to exist in classical systems, such as the Sinitsyn-Nemenman (SN) model of reaction kinetics, in which some parameters are cyclically modulated by an external agent [13,14]. Recently, the adiabatic result has been extended to the non-adiabatic regime also for this model [24]. It was found that the pumping current reaches a peak after the adiabatic regime in which the current linearly increases with the modulation frequency, and eventually decays as the inverse of the modulation frequency in the asymptotic limit. Finite modulation speed means that the pumping current can no longer be expressed using strictly geometrical quantities, but a formally geometrical expression in terms of a line integral in parameter space is still possible. Furthermore, the effect of non-adiabaticity on the fluctuation theorem has been investigated in the context of the SN model [25]. Motivated by the attention these recent results have garnered, in the present paper a non-Markovian generalization of the SN model is presented as a natural extension of the aforementioned research, and its adiabatic and non-adiabatic behavior is investigated analytically and numerically.
The structure of this paper is as follows: In Sec. II, we first introduce our Markovian N + 4 model and then perform a tracing out procedure to obtain an N -state non-Markovian system. The results obtained for the pumping current are shown in Sec. III. We obtain the geometrical form of the adiabatic current in the low frequency limit and compare this to the full non-adiabatic current obtained using numerical methods. In Sec. IV, the results obtained and their physical implications are discussed, and a conclusion is presented in Sec. V. Supplementary information and details of the various calculations in the main text will be presented in the Appendices.

II. NON-MARKOVIAN SINITSYN-NEMENMAN MODEL
In this section, we introduce our non-Markovian model by deriving it from a Markovian model having a larger transition matrix. This procedure is analogous to the Nakajima-Zwanzig method in which we obtain a system obeying non-Markovian dynamics by tracing out degrees of freedom [26][27][28]. However, as will be discussed below, we can solve the Markovian model with a larger matrix in some situations, because the larger Markovian system is phenomenological and governed by simple dynamics. Therefore, in contrast to the relation between the full von-Neumann equation and the Redfield equation or Lindblad equation for a subsystem in the Nakajima-Zwanzig method, in our case the relation between the Markovian model and non-Markovian model is simple, affording exact conversion between them.

A. Markovian Formulation
Let us consider a Markovian master equation for a system with N + 4 states, where ρ i is the probability of state i and the element L ij of the transition matrix is the hopping rate between states i and j. For convenience, we have introduced the inverse time scale (characteristic hopping rate) k 0 so that all other quantities except for time t become dimensionless. The time dependence of the hopping rates arises from modulation by an external agent. This can be realized for instance by controlling the temperatures or chemical potentials of the environments, though the temperature control is not easy in realistic situations [29]. In the following, the rates will be assumed to undergo cyclic modulation in time. We also assume that the environments attached to the system are large enough such that the back-reaction from the system is negligible and the probability distribution in each environment is always described by an equilibrium distribution ρ ν,eqm i , where ν = L or R is the index to specify the left (L) or right (R) environment attached to the system, even when some parameters such as the temperature and chemical potential are controlled by an external agent.
We first briefly review the properties of the ME required for conservation and non-negativity of the probabilities as well as detailed balance to hold. For probability to be conserved, the matrix elements must satisfy [1] The condition for complete positivity is more restrictive: we are dealing with a continuous process in time, which means that the only way for a probability, say ρ i , to become negative is if it first passes through zero. We thus require [1] d dt ρi=0 for arbitrary non-negative ρ i . This in turn means that all off-diagonal elements of L ought to be non-negative. Finally, we impose detailed balance. The principle of detailed balance states that a process and its reverse happen at the same rate at equilibrium: The equilibrium probability distribution would usually have the Boltzmann form characterized by the inverse temperature β and energy E i of the state i, i.e., ρ eqm i ∝ e −βEi so that the above equation can be written as Let us explore what this means for our ME. Summing over j, we obtain where we used Eq. (2). Thus, we arrive at This is not only a statement of the intuitive notion that an equilibrium state should exist, but also an important condition to determine the hopping rule between the system and the environment. It can be guaranteed to hold by requiring det L = 0. (8) Note that this condition is already guaranteed by Eq. (2). Thus, requiring detailed balance to hold does not impose additional constraints. We note that the Markovian SN model treated in e.g. Ref. [13] corresponds to a Markovian two-state model. Now let us rewrite the Markovian ME in matrix notation, where we have introduced the phase variable θ := Ωk 0 (t− t 0 ) with the dimensionless angular frequency Ω of cyclic modulation and the time t 0 required to reach a cyclic state which is independent of a specific choice of initial conditions. We assume that the system can be partitioned into one N -state site and an adjacent two-state site on each side: i.e. the N -state reduced system p and the left and right adjacent states q L and q R , respectively. We note here that this partitioning can be considered analogous to the projection operators P and Q of the Nakajima-Zwanzig approach, though the sizes of adjacent states q L and q R are small in our model. Provided that the adjacent states do not couple to each other directly, the ME becomes Here, S(θ) is a N ×N square matrix, W ν in (θ) are N ×2 and W ν out (t) are 2 × N rectangular matrices, and K ν (θ) with ν = L, R are 2 × 2 square matrices. This small matrix K ν (θ) enables us to understand the relationship between non-Markovian dynamics in a subsystem and Markovian dynamics in the full system.
Next, we specify the structure of the submatrices introduced above. Firstly, K ν (θ) are assumed to be given by where, to make the model concrete, we choose [24] k where the internal rate constant w > 0, amplitude r, phase difference φ and dimensionless modulation frequency Ω of Eq. (9) are used as control parameters. This choice can be motivated as follows. Firstly, in the space of the rate coefficients, the protocol trajectory must enclose a finite area for finite pumping current to be observed in the adiabatic limit [24]. The phase difference φ in the last equation of Eq. (13) ranges from 0 to π/2. Below we mainly focus on the choice φ = π/2 which gives a circle of radius r centered at k 1 = 1, k 2 = 1, but for illustrative purposes we will also consider a straight line which does not enclose any area, obtained by setting φ to zero. Secondly, as can be seen from Appendix A, to obtain a finite pumping current, there needs to be a phase difference between the left and right rates as well as the rates of the two states on each side.
We note that the concrete calculations below are performed for N = 2, where S, W ν in and W ν out were assumed to be proportional to the 2 × 2 identity matrix 1: S = −2w1, W ν in = W ν out = w1, though we do not specify N in the general framework. The full transition matrix for N = 2 is shown explicitly in Eq. (A2).
As can be checked by direct calculation, the Markovian equation (11) together with the above structure of the submatrices guarantees the conservation and nonnegativity of ρ, thus also satisfying detailed balance. Physically, an important feature of this model is that the two states at each site, e.g. q R 1 and q R 2 for the right adjacent site, can only mix in the left and right environments, so there is no coupling between them within the system. Furthermore, we assume that only k ν 1 and k ν 2 have time dependence. Physically, this corresponds to externally modulating only the coupling between the adjacent states and the environment.
We note that K ν in Eq. (12) can be split as where, K 0 and K L 1 (θ) are, respectively, and K R 1 (θ) is obtained by exchanging cos θ with cos(θ − φ). Furthermore, with the quantity of interest, the pumping current, in mind, we introduce the counting field χ: We may now consider the modified ME, where the partial derivative is used to emphasize the dependence of ρ on χ. The modified transition matrix L(θ, χ) is defined as the original matrix with K R (θ) replaced by K R (θ, χ); its explicit form is shown in Eq. (A2). Thus, the modified probability vector ρ(θ, χ) is simply defined as the solution of the modified ME, Eq. (17). Solving for ρ(θ, χ), computing the cycle averaged characteristic function z(θ), where and defining the cumulant generating function, will now give us access to all moments of the pumping current. It can be readily shown that the pumping current, defined as the first moment, can equivalently be obtained by taking the cycle average of the instantaneous current: where Physically, this corresponds to the current flowing from the right adjacent states to the right environment. In the adiabatic limit, this pumping current has a geometrical formulation and can be obtained analytically, as shown in the next section.

B. Non-Markovian Form
It is well-known that memory effects arise when we trace out environmental degrees of freedom. Physically, when the environment is allowed to back react with the system, the environment will retain information about the past states of the system, which can affect its dynamics.
We now demonstrate how this arises for our system. Specifically, we show that the N +4 -state Markovian ME is equivalent to an N -state non-Markovian ME. First, we rewrite the Markovian ME, Eq. (11), in the split form Next, we formally solve for q ν , where transient terms are neglected, and the time evolution operator U ν (θ) obeys the differential equation While U ν (θ) can be written with the help of the time ordered exponential operator, this expression is not useful for concrete analysis. Nevertheless, if we restrict our interest to the case of small r, Eq. (27) is reduced to an exponential function. Accordingly, we expand U ν (θ) as which, upon substituting into Eq. (26) and solving at each order, leads to Substituting this into the differential equation for p, we indeed obtain a non-Markovian ME of the form where the memory kernel is given by On the other hand, we may obtain coupled equations for q ν (t) instead, as We emphasize that we have not used any approximations in deriving these time-nonlocal equations, so the important properties of the original Markovian model, namely that ρ remains non-negative and conserved throughout the time evolution, are retained. However, the effectiveness of the above procedure hinges on there being only a few adjacent states q ν connecting the system to each environment. In the general case, as treated by the Nakajima-Zwanzig formalism, tracing out a large number of irrelevant degrees of freedom does not lead to an analytically workable equation. In many situations, the memory kernel M (θ, θ ′ ) is described by a multiple exponential function. In fact, it is known in the literature on Markovian embedding of non-Markovian processes that an exponential memory kernel is the easiest case to treat analytically [30,31]. However, our model does not lead to a multiple exponential function, but a more complicated form due to the time dependence introduced by the external modulation and small adjacent vectors q L and q R , and thus, treating this system anaytically is non-trivial.
The advantage of our approach is that there is no need to deal with the non-Markovian dynamics: we can simply solve the equivalent Markovian system because the eliminated adjacent systems are small. Working with the non-Markovian form of the equations does not present any mathematical or physical merit. Indeed, we will not employ the non-Markovian form of the dynamics in the treatment that follows. However, deriving the timenonlocal form serves as an interesting example of how non-Markovianity arises naturally in the context of externally modulated systems.
There are a number of alternative ways to deal with nMMEs that have been explored in the literature. The most straightforward approach is to transform the nMME into Laplace space and solve the resulting algebraic equation there. However, it often happens that M (θ, θ ′ ) does not depend only on the phase difference (as is the case here), so that the convolution theorem cannot be utilized. This issue can be circumvented by performing a Taylor expansion of M (θ, θ ′ ) around one of the time variables so as to create terms which only depend on the difference between θ and θ ′ [32]. While this allows transformation into Laplace space, it is generally difficult to perform the inverse transformation explicitly. Furthermore, the resulting solution is in the form of an infinite series instead of the closed-form approach of the present paper.
A more general approach to deal with nMMEs is based on the Nakajima-Zwanzig projection operator technique [33][34][35]. In essence, one first eliminates the environment dynamics (corresponding to q ν of Eq. (24)) to obtain a nMME, and then performs an expansion in the systemenvironment coupling to achieve a time-covolutionless (TCL) equation of motion for the system [26][27][28]. Again, however, while a number of refined expansion protocols have been developed over the recent years [36], a closedform time-local equation which is easy to work with numerically or analytically cannot be derived using this approach.

A. Analytical expression of the adiabatic current
As can be seen from Eqs. (A14) and (A15) in Appendix A, we can write the pumping current in terms of the Berry-Sinitsyn-Nemenman (BSN) curvature in the adiabatic limit: where the BSN curvature is written as The Berry vector potential A α with α = 1, 2 is defined in terms of the eigenvectors of the modified transition matrix L(θ, χ) in Eq. (A15). In the above equations, the rates k 1 and k 2 are treated as functions of the radial and azimuthal parameter which prescribe the area to be integrated over. Note that F is a second-rank antisymmetric tensor in general, but in the case of only two independent modulation parameters, this can readily be reduced to a scalar quantity. A surface plot of the curvature together with the circular modulation trajectory (for r = 0.9 and φ = π/2) is shown in Fig. 1.
The parametrization of Eq. (13) has the Jacobian | sin φ|r, so we have Let us first consider expanding around φ = π/2; a straightforward integration for N = 2 shows that It can be confirmed that the adiabatic current is negative for all frequencies. Next, we expand the current for N = 2 around φ = 0, which gives . (38) This indicates that in the absence of phase difference, no geometrical current is generated. While the above equations hold only for the adiabatic limit, the numerical results shown below support this conclusion also for the non-adiabatic regime.

B. Numerical results
To check the validity our analytical findings, we also numerically solve the Markovian ME, Eq. (24), of the N = 2 system. The computations were performed for frequencies below Ω = 20, and for a range of the system parameters w, r and the phase difference φ. It is found that to guarantee convergence and independence from initial values, 40 initial cycles, i.e. k 0 t 0 = 80π/Ω is sufficient. Equations (23) and (22) are used to compute the pumping current for each solution.
We first show plots of the pumping current for selected values of w with modulation amplitude r = 0.9 and phase difference φ = π/2 in Fig. 2. From these plots we see that the agreement between the numerical averaged current (solid lines) and the adiabatic current expressed by Eq. (37) (dashed lines) is good in the adiabatic regime Ω → 0.
We have also investigated the dependence of the pumping current on the modulation phase difference φ in the low frequency regime (Ω = 0.1). Results are shown in Fig. 3 , where the solid circles indicate the numerical result, obtained again by direct solution of Eq. (24), the dashed lines correspond to the analytical result for φ ≃ π/2 of Eq. (37) (including O(φ 3 )) and the solid lines show the result for φ ≃ 0 as shown in Eq. (38) (including O((φ − π/2) 3 )). We see that both expansions agree well with the numerical results, and together the two expansions capture the φ-dependence of the pumping current in the adiabatic regime sufficiently well.
Finally, we investigate the dependence of the pumping current on the modulation amplitude r and the internal rate factor w again for the adiabatic regime (Ω = 0.1), as shown in Figs in this regime, the agreement between the numerical and analytical results is good.
Comparing our results with those reported e.g. in Ref. [24], we can see that the qualitative features of the current are unchanged. We note, however, that, due to the fact that our system is fundamentally different from the two-state Markovian SN model characterized only by a 2 × 2 matrix, detailed comparison is not possible.
We also note that the above analytical calculations can in principle be extended to the non-adiabatic regime by using a more general form of Eq. (A13) that allows for transitions between different eigenstates. This approach is investigated in Appendix B. The higher order eigenvectors have a complex form and the required integrations rapidly become involved. As shown in Appendix B, we have confirmed the quantitative validity of our nonadiabatic analysis, at least at O(Ω 3 ). Compared to the ease of the above numerical treatment, however, these calculations do not yield any worthwhile insight into the physics, nor do they seem to offer a significant improvement over the adiabatic calculations.

IV. DISCUSSION
The most important finding of our research is that a large class of systems possessing memory can be embedded into a larger Markovian system. The computational cost incurred due to the increased number of dynamical variables is outweighed by not having to deal with the time convolution integral, at least for the model we analyzed (see Eq. (11)). This is certainly true for numerical calculations, but we have shown that it also applies to the analytical calculations in the adiabatic limit.
Furthermore, analyzing the system in the non- Markovian framework does not have any apparent merits.
Involving the time convolution integral in the discussion obscures the definition of the pumping current; indeed, it is not clear how to implement full counting statistics (FCS) in the presence of memory effects. On the other hand, since the larger Markovian system is not merely an approximation, but rather an exact equivalent of the non-Markovian system, reliability is not compromised by limiting oneself to the Markovian system. Thus, direct treatment of an equation of the type of Eq. (31) would be required only when the embedding into a larger Markovian system is not feasible.

V. CONCLUSION
In this paper, the Sinitsyn-Nemenman (SN) model, a periodically modulated N -state system, was generalized. It was shown that the non-Markovian SN model, governed by a non-Markovian master equation (nMME) which includes a time convolution integral, can be converted to a larger Markovian system, which is easier to analyze. Thus a method of solving the system dynamics at least numerically without needing to resort to any perturbative expansions was presented, yielding an approach to deal with this type of nMMEs which is in principle exact. In addition to solving the master equation numerically and using the results to calculate the pumping current for different values of the control parameters of the system, such as modulation frequency, phase and amplitude of the external modulation, we also presented an analytical calculation of the geometrical current in the adiabatic limit. It was found that the agreement between the numerical results and the analytical calculation was good.
Prospects of further research into this problem in-clude the following: detailed microscopic derivation of the nMME in the framework of the present model, studying how the fluctuation theorem is affected by memory effects, and applying the method developed here to periodically driven quantum mechanical models. In particular, it would be interesting to see whether the dynamically modeled environment could be interpreted as the diagonal part of the density matrix of a two-state quantum system, thus leading to a connection between quantum coherence and non-Markovian time evolution. We are also interested in the extension of this analysis to a quantum mechanical system with charge current induced by the BSN curvature. It should be noted that for such systems, the pumping charge current becomes zero if the repulsive interaction between electrons is not considered, but it has a finite value in the presence of such interactions [16,37]. We also note that a quantum master equation can be derived by using Green's function technique [16].

ACKNOWLEDGMENTS
The authors would like to thank Kazutaka Takahashi, Kazunari Hashimoto, Yuki Hino, Kiyoshi Kanazawa and Hiroyasu Tajima for fruitful discussions. This work is partially supported by Grant-in-Aid of MEXT for Scientific Research (Grant Nos. 16H04025 and 21H01006), a MEXT scholarship and ISHIZUE 2020 of Kyoto University Research Development Program.
Appendix A: Details of the 6-state model In this Appendix, we consider the 6-state system (i.e. with N = 2) obeying the Markovian ME of Eq. (24) in more detail, and outline the results used in the main text.
We begin by obtaining the BSN curvature tensor using the machinery of full counting statistics (FCS) [38]. Firstly, we note that in vector notation, after adding the counting field χ to the transition rates between the right adjacent states and the right environment, the Markovian ME reads where the full modified transition matrix L(θ, χ) is written as (A2) Figure 6 shows a schematic representation of this system: as described in the main text, there is no coupling between the two states at each site, except at the left and right reservoirs. Furthermore, we assume that only k ν 1 and k ν 2 depend on time. We note here that, together with the observation that an alternative formula for the pumping current is it can be seen from the above explicit form that unless k R 2 (θ) = k R 1 (θ), the pumping current vanishes, as mentioned in the main text below Eq. (13). In more detail, having no phase difference between these hopping rates means that the cumulant generating function becomes O(χ 2 ), resulting in the first derivative at χ = 0 evaluating to zero.
(A4) In general, the above quantities cannot be obtained exactly, but as we are only interested in their first derivatives with respect to χ here, we find approximate expres- sions Since λ(θ, 0) = 0, we have the approximate eigenvalue λ(θ, χ) ≃ λ 1 (θ)χ, which we wish to compute explicitly. We first note that the characteristic equation of the modified transition matrix for the eigenvalue λ has the form det[L(θ, χ) − λ1] = a n λ n + ...
where a n are functions of the matrix elements of L(θ, 0). This form guarantees that when χ = 0, λ = 0 is a solution of the characteristic equation. Thus, substituting in λ = αχ and ignoring all terms of O(χ 2 ) or higher, we obtain λ 1 = −b 0 /a 1 . This eigenvalue approximation can then be used to obtain an approximation for the right and left eigenvectors. Writing down the eigenvalue equation to first order in χ, we have and similarly for the left eigenvector. These systems of linear equations can be solved up to normalization algebraically. Next, we need to impose the normalization ℓ † (θ, χ)r(θ, χ) = 1 up to first order in χ, achieved by the normalization factor r(t, χ) → N (θ, χ)r(θ, χ), N (θ, χ) ≃ N 0 (θ) + N 1 (θ)χ.
(A11) We clearly must have which shows that N 0 Finally, we apply the above results to the calculation of the pumping current in the adiabatic limit. Assuming the normalization ℓ † (θ, χ)r(θ, χ) = 1, as long as Ω is sufficiently small, we can show that where v(θ, χ) := ℓ † (θ, χ) ∂ ∂θ r(θ, χ) gives the geometrical contribution. Note that with our choice of modulation protocol, the dynamical term λ(θ, χ) averages to zero, and we are left with just v(θ, χ). Using this approximation to compute the generating function of the current, and then picking out the term linear in χ (that is, taking the derivative with respect to χ and setting χ to zero), we obtain where we performed a change of variables from time to the modulation parameter space spanned by k 1 and k 2 . This is equivalent to Eq. (23). Here the Berry vector potential is defined as Finally, employing the generalized Stokes' theorem, which in the case of two-parameter modulation reduces to Green's theorem, we indeed obtain Eq. (35). This means that to obtain the adiabatic current, we need to compute λ(θ, χ), r(θ, χ) and ℓ † (θ, χ) only, and only up to first order in χ.
Appendix B: Eigenvalue expansion beyond the adiabatic limit In this section we consider an approach based on the full eigenvector expansion of the N = 2 system [24]. We begin with obtaining all 6 eigenvalues and the corresponding right and left eigenvectors. Note that in this Appendix, the subscripts refer to the index of the eigenvalue or eigenvector, and series expansions are indexed by bracketed superscripts. Plots of the eigenvalues and each component of the right eigenvectors as a function of θ are shown for w = 2, φ = π/2, and r = 0.9 in the figures below. We note that while for this choice of parameters, all the eigenvectors remain well-behaved for all θ, exceptional points where the dimension of the vectors space spanned by the eigenvectors drops can arise for some modulation parameters. This may limit the usability of the approach presented here.
We normalize all eigenvectors as follows, [ FIG. 9. The 6 components of the eigenvector r2 of the N = 2 system plotted against θ, with r = 0.9, φ = π/2 and w = 2. and hence, necessarily we have Transforming the eigenvectors according to whereṙ := ∂ ∂θ r, leaves orthonormality and completeness unchanged, but results in the convenient propertỹ  Noting that the terms multiplying the constants of integration C i decay exponentially and can thus be neglected for long times, the expansion becomes The numerical integrations required by Eq. (B11) are rather slow, so let us obtain δ i recursively instead. Defining δ (0) i := Ωl † iṙ 1 λi , we obtain To investigate the usefulness of this approach, we show the pumping current obtained for n = 0, i.e. the adiabatic limit (solid line) and n = 2, i.e. up to O(Ω 3 ) (note that only odd powers of the frequency contribute to the current) compared to the numerical result from the main text in Fig. 14. We see that the adiabatic result is correctly reproduced by the eigenvalue analysis, and the addition of the O(Ω 3 ) offers an improvement over the linear approximation.