Non-adiabatic effect in quantum pumping for a spin-boson system

We clarify the role of non-adiabatic effects in a quantum pumping for a spin-boson system. When we sinusoidally control the temperatures of two reservoirs with \pi/2 phase difference, we find that the pumping current strongly depends on the initial condition, and thus, the current deviates from that predicted by the adiabatic treatment. We also analytically obtain the contribution of non-adiabatic effects in the pumping current proportional to \Omega^3 where \Omega is the angular frequency of the temperature control. The validity of the analytic expression is verified by our numerical calculation. Moreover, we extend the steady heat fluctuation theorem to the case for slowly modulated temperatures and large transferred energies.


Introduction
A pump converts an external bias into work. We need the average bias to get the work from a macroscopic mechanical pump, but it is known that the average bias to get pumping current is not necessary in a mesoscopic systems. When a mesoscopic system, thus, is slowly and periodically modulated by several control parameters such as chemical potentials, gate voltages and tunneling barriers, there exists a net average current without dc bias. This phenomenon is known as the adiabatic pumping, and has been observed in various processes such as quantized charge transport [1][2][3][4][5][6][7][8][9][10][11], spin pumping [12][13][14][15][16][17][18] and quibit manipluation [19]. The first proposal of the adiabatic pumping is given by Thouless[1] for a closed quantum system. The idea of the quantum pumping for closed systems has been extended to for open systems [9,11,[20][21][22][23][24][25][26][27][28]. Such adiabatic pumping processes have been experimentally realized in mesoscopic transport processes [4-8, 18, 29-32]. It is recognized that the mechanism of adiabatic pumping is originated from the geometrical effect of Berry phase in quantum mechanics [33], where a circular operation in a parameter space creates a non-zero geometrical quantity associated with the pumping current.
Most of the previous studies, however, assume that the pumping process is only modulated adiabatically, where the validity of the approximation is ensured if the modulation speed is zero. This situation is practically useless, because pumping current under the adiabatic modulation is zero in strict sense. It is, thus, important (i) to clarify the limitation of the adiabatic approximation and (ii) to analyze the pumping process without the introduction of the adiabatic approximation to get finite pumping current under a finite speed modulation.
Although there exist some papers to discuss non-adiabatic pumping effects based on a stochastic equation with a weak noise [53], the master equation [54], the Floquet scattering theory [55], and Green function [56,57], it is unclear how non-adiabatic effects affect the pumping current. Indeed, it is known that a non-adiabatic process can cause a phase transition through the analysis of a simple quantum mechanical model [58].
We may ask another non-trivial question associated with a non-adiabatic pumping process besides the pumping current. Although there exists the heat fluctuation theorem [59][60][61][62][63][64], at least, for adiabatic dynamics of open Markovian processes, the heat fluctuation theorem seems to be violated under some situations such as the dynamics under modulated external fields [42], or a non-Gaussian noise [65] or the dry friction [66]. We have to clarify the reason why the heat fluctuation theorem seems to be violated.
In this paper, we systematically study non-adiabatic pumping effects within the framework of the quantum master equation under the Markovian approximation. For this purpose, similar to Ref. [54], we analyze the simplest spin-boson model under the weak coupling condition between surrounding environments and the system. We continuously control the temperatures in the environments with the modulation frequency Ω/2π, we clarify initial condition dependence of the pumping current and the essential non-adiabatic effects to the pumping current. We also extend the steady heat fluctuation theorem to the cases for slowly modulated temperatures and large amount limit of the transferred energy.
The organization of this paper is as follows. In Sec. 2, we introduce a model of a spin-boson system and the methods of the generalized quantum master equation with the full counting statistics (FCS). Section 3 is the main part of this paper, which consists of three part. In Sec. 3.1, we derive general expressions of non-adiabatic pumping current. In Sec. 3.2, we apply our formulation to the spin-boson system introduced in Sec. 2, and present the results for the pumping current to clarify the non-adiabatic effects. In Sec. 3.3, we discuss whether the heat fluctuation theorem is still valid. Finally, we discuss and summarize our results, in Sec. 4. In Appendix A, we briefly summarize the properties of the cumulant generating function and the first moment. In Appendix B, we derive the master equation with parameters modulation in the context of FCS. In Appendix C, we reproduce the adiabatic Markovian pumping current obtained in Ref. [42] within our framework. In Appendix D, we summarize the relationship between our formulation and that in Ref. [42]. In Appendix E, we derive the asymptotic expansion of the density matrices and the non-adibatic pumping current. In Appendix F, we explain the detailed derivation of the extended heat fluctuation theorems showed in Sec. 3.3.

Model and Method
In this section, we introduce our model and the method to be used in our analysis. We analyze a spin-boson system, and adopt the generalized quantum master equation with the full counting statistics (FCS) as a basic equation for our analysis.

2/26
The spin-boson system is a simple two-level system {|0 , |1 } coupled with two environments (denoted as L and R) characterized by the inverse temperatures β ν where ν = L or R. We modulate the temperatures periodically with the angular frequency Ω under the condition that environments are always equilibrium. The system Hamiltonian H S and the environmental Hamilitonian H ν E (ν = L or R) are, respectively, given by where b k,ν and b † k,ν are, respectively, bosonic annihilation and creation operators at the wave number k for the environment ν, ǫ n and ω k,ν are the energy for the level n(= 0, 1) and the angular frequency to characterize the bosonic environment ν, respectively. We introduce the characteristic frequency ω 0 from the relation ω 0 ≡ ǫ 1 − ǫ 0 . The interaction Hamilitonian H ν SE is given by with the coupling strength g k,ν , which is characterized by the spectral density function Γ ν (ω) = 2π k g 2 k,ν δ(ω − ω k,ν ). We assume that the environments are always characterized the equilibrium operator ρ eq E (β ν ) = e −βνH ν E /Z. To calculate the average of energy transfer ∆q t from a reservoir to the system during the time interval t, we use the method of FCS. When the two-point projective measurement on a quantity Q is performed at times 0 and t, the corresponding outcomes are q 0 and q t respectively. Thanks to the method of FCS, we can calculate the cumulant generatig function S(χ, t) ≡ ln P (∆q t )e iχ∆qt d∆q t , where P (∆q t ) is the probability distribution function of ∆q t = q t − q 0 and χ is the counting field. Once we know S(χ, t), we can get n-th cumulant of P (∆q t ) from n-th derivative of S(χ, t) at χ = 0. Therefore the average energy transfer is given by ∆q t c = ∂S(χ, t)/∂(iχ)| χ=0 . The detailed method of the calculation of the cumulant generating function S(χ, t) is explained in Appendix A. In this method, the cumulant generating function is given by S(χ, t) = Trρ tot (χ, t), where ρ tot (χ, t) is the generalized density matrix for the total system defined in Eq. (A9).
In the weak coupling limit g k,ν ≪ ω k,ν , ǫ n / , it is straightfoward to obtain the quantum master equation for the reduced density matrix ρ(χ, t) ≡ Tr E ρ tot (χ, t) (see Appendix B). According to Appendix B, the correlation time scale τ C of environments is characterized by the symmetrized time correlation function, which is, for the operator of environtments B ν = k g k,ν b k,ν in our model, given by (see Ref. [67]) where we have used the Bose distribution b † k,ν b k,ν = (e βν ωk,ν − 1) −1 and the Ohmic spectral density Γ ν (ω) = gωe −ω/ωc,ν with the cutoff ω c,ν . Here, ℜA represents the real partof A and ψ(x) ≡ Γ ′ (x)/Γ(x) is the digamma function. For our setting of parameters in this paper, the characterized time scale τ C in Eq. (3) satisfies τ C ω 0 ∼ β ν ω 0 ∼ O(1). On the other hand, the relaxation time scale τ R of the system is estimated as τ R ω 0 ∼ 10 3 and we 3/26 consider τ R Ω −1 . Therefore, if the condition τ C ≪ τ R Ω −1 is satisfied, we can derive the Markovian quantum master equation where [H, A] ≡ H χ A − AH −χ for an arbitrary operator A and H χ ≡ e iχQ/2 He −iχQ/2 , and β(t) is the vector representation of {β ν (t)}. In the Markovian case, spectral density Γ ν (ω) is reduced to the constant tunneling rate Γ ν ≡ Γ ν (ω 0 ). Because we consider the identical environments, let us introduce Γ ≡ Γ L = Γ R , which characterizes the relaxation time scale τ R ∼ Γ −1 of the system. Let |ρ(χ, t) be the vector |ρ(χ, t) ≡ T ( 0|ρ(χ, t)|0 , 1|ρ(χ, t)|1 ) consisting of the diagonal element of ρ(χ, t) with the notation of the transverse T A of an arbitrary vector A. Note that the diagonal part of Eq. (4) can be independet of off-diagonal part in our model. Thus, the quantum master equation (4) can be written as where the evolution matrix K χ M (β(t)) is given by Here, we have introduced where Here, we explicitly write the control parameters, inverse temperatures β, and · β represents the average over the bosonic field in the environment characterized by β. Namely, we have 4/26 assumed that the environments are always in thermal equilibrium even if we modulate β(t). Thus, it is not appropriate to apply our formulation to too fast modulations. This means that we cannot use our theory for the case of abrupt change of the temperature. We also assume that the time evolution of β(t) satisfies where T 0 and T A are, respectively, the average temperature and the amplitude of the modulation.

General Expression
It is straightfoward to extend the adiabatic approximation used in Appendix C which is reduced to that used in Ref. [42]. At first, let us decompose the average current into two parts: which is a natural extension of Eq. (C5). It should be noted that the contribution from the dynamical phase ∆q t d is invariant even in the non-adiabatic treatment, while the adiabatic geometrical current in Eq.(C6) is now replaced by ∆q t E na : where l ′ + | is the χ-derivative at χ = 0 of l χ + |, which is the left eigenvector of K χ M for the eigenvalue λ χ + with the maximum real part. Namely, the right eigenvector |λ 0 + for the steady state used for the adiabatic process is relaced by the density matrix |ρ(0, t) . This result can be interpreted as follows. Because the dynamical phase depends only on the average bias for the symmetric cyclic modulation, it is reasonable that the pumping current through the dynamical phase is unchanged, even if we consider the non-adiabatic effects. On the other hand, the adiabatic transfer ∆q g a depends on modulation speed and the path on the parameter space. In the non-adiabatic case, hence, the excess energy transfer corresponding to the contribution from the geometrical phase in the adiabatic limit has to be replaced by ∆q t E na . Now, let us prove the expressions (15) and (16). The formal solution of Eq. (5) is where we have introduced the time-ordering product from left to right as ds n K(s 1 )K(s 2 ) · · · K(s n ). Thus, we obtain the expression of the energy transfer ∆q t as : where we have used ∆q 0 = 0 and the following deformation under the condition 1|K 0 M (β(τ )) = 0 : Equation (18) can be rewritten as The first term on the right hand side (RHS) of this equation is equal to ∆q t d by using λ 0 + = 0, and the second term on the RHS of (20) is reduced to Eq.(16) with the aid of Eq.(5). Thus, we reach Eq.(15).

Application to Spin-Boson System
We, now, apply our formulation to the spin-boson system. Let us introduce a χ is given in Eqs. (7)- (10). We consider the case that the measured quantity is the Hamiltonian in the right environment Q = H R E . Thus, the explicit form of each a χ i (t) in this case is given by: where n ν (t) = (e βν (t) ω0 − 1) −1 . The eigenvalues and the eigenstates of K χ M (β(t)) are explicitely written as where we have introduced in Eq. (26). These expressions satisfy the orthonormal condition l χ i |λ χ When the counting field χ is absent, the above results reduce to |λ 0 where λ(t) ≡ a 1 (t) + a 4 (t), and .
To derive Eqs. (29)-(31) we have used the trivial relations a 0 2 (t) = −a 4 (t) and a 0 3 (t) = −a 1 (t). By solving Eq. (5) under the condition χ = 0, one of the diagonal components of the density matrix is given by where ρ 11 (t) represents 1|ρ(0, t)|1 . The time evolution of ρ 11 (t) at The initial relaxation of ρ 11 (t) against Γt for Ω/Γ ≤ 0.1. All numerical data are collapsed on the red line given by Eq. (34). (ii) The long time behavior of ρ 11 (t) as the plot against Ωt/2π, where the red line represents the adiabatic form ρ ad (t).
Because we consider a system that has symmetric junctions between the system and the environments under the no-average bias, it is easy to show that ∆q t d is zero (see (D7)).
We, thus, plot the non-adiabatic pumping current which is defined by ∆q t E na /τ p , and the adiabatic one against the frequencies of modulation Ω with the numerical calculation and the asymptotic expansion (Fig.2). It should be noted that large deviation between the adiabatic current and the obtained current is mainly originated from the initial condition dependence. In other word, if we start the measurement of the current after t ≫ λ(0) −1 , the adiabatic approximation gives a reasonable result in wide range of Ω.  2 The plot of the pumping current against Ω/Γ. The set of parameters is equivalent to that used in Fig.1. The difference between the analytic result obtained from the first order of the asymptotic expansion of J E na and the numerical result represented by the red dots is invisible. 8/26 As shown in Appendix E, the asymptotic expression of the pumping current is given by where Thus, J E na in the lowest order in Ω is reduced to the adiabatic pumping J g a if we begin with ρ 11 (0) = ρ ad (0). If we begin with ρ 11 (0) = ρ ad (0), however, the expression of the adiabatic current does not give the correct result of the pumping current. To verify the results in Eqs. (41)-(44), we explicitly plot how the pumping current depends on the initial condition (Fig.3), where the analytic result (solid line) perfectly reproduces the numerical results. In the case ρ 11 = 0.5, we plot the result of the pumping current in Fig.4(a). To remove the initial condition dependence, we alseo plot the result of the pumping current for the initial measurement starting from t ≫ λ(0) −1 in Fig.4(b), where we start the measurement 9/26 from τ p , which initial condition corresponds ρ 11 (0) = ρ ad (0). It is visible that the adiabatic current gives the reasonable result for Ω/Γ < 0.2, but there exists a little systematic deviation between the linear or adiabatic result and the numerical result for Ω/Γ > 0.2. To clarify the non-adiabatic contribution up to Ω 3 , we plot the J − B 0 Ω (see Fig.5(a) and (b)). It is obvious that our analytic non-adiabatic expression in Eqs. (41)- (44) gives the reasonable result even if the linear expression in J E na is no longer valid.  (41) and (42). The set of parameters is equivalet to that used in Fig.1  nu eri al up to 3rd order in Fig. 5 The non-adiabatic pumping current obtained from the subtraction of the first order term B 0 Ω to clarify the term up to the third order for (a) ρ 11 (0) = 0.5 (b) ρ 11 (0) = ρ ad (0). We use the identical parameters to those used in Fig.1. The non-adiabatic pumping current obtained from the analytic calculation (blul line), asymptotically reproduces the numerical result (red dots) in the small region of Ω/Γ. 10/26

Extended Fluctuation Theorem
In this subsection, we discuss whether the fluctuation theorem for the heat currents exist in the our system under the existence of dc bias, i.e. ∆β ≡ β R − β L . In this sense, the set up on the heat fluctuation theorem in this subsection differs from the case without dc bias discussed in the other parts.
First, we consider the case without the temporal change of the temperatures, i.e. β ν (t) = β ν . Then, we readily obtain J g a = J E na = 0 because ρ(χ, t) decays with negative eingenvalues in Eq. (25). Therefore, the steady cumulant generating function is reduced to S(χ) ≡ lim τ →∞ 1 τ S(χ, τ ) = λ χ + , and the product of a χ 2 and a χ 3 satisfies the Gollavotti-Cohen (GC) symmetry [69]: where α = ln nL(1+nR) nR(1+nL) / ω 0 = ∆β. If the cumulant generating function satisfies GC symmetry, the steady fluctuation theorem is held: Nevertheless, the simple GC symmetry for the cumulant generating function is no longer held when the temperature varies with time even in the adiabatic limit because the constant α is replaced by α(t) and the contribution of the geometrical term exists. This may give the basis of the violation of the heat fluctuation theorem in Ref. [42], and the geometrical entropies introduced in Refs. [43,52]. In the case of the periodic change of temperatures, there is room to choose a variable of large deviation to obtain the correlation of the fluctuation theorem.
On the other hand, for N = ∆q/ ω 0 → ∞, the fluctuation thorem for the transferred energy is given by where u ν (θ) = (5 + 2n ν (θ))dβ ν (θ)/dθ and P B is the time-reversal distribution against P F as shown in Appendix F.2. Numerical verification of the extended fluctuation theorems will be reported elsewhere. Nevertheless, we believe that the extend fluctuation theorem in Eq. (47) is, at least, universal for the slowly modulated case. Indeed, the derivation of Eq. (47) does not contain any specific feature of the spin-boson system.

Discussion and Conclusion
We have successfully extended the theory for the adiabatic pumping to the non-adiabatic pumping for finite speed modulations within the framework of the Markovian quantum master equation. We have applied our formulation to the spin-boson system and found that (i) the pumping current strongly depends on the initial condition, and (ii) the contribution of the non-adiabatic pumping current is relevant for relatively large Ω/Γ if the contribution of the initial relaxation is eliminated. (iii) The contribution of non-adiabatic effect has is analytically reproducible in terms of the technique of the asymptotic expansion. (iv) The extended fluctuation theorems for slowly modulated temperatures and large transferred energies are derived.
Our master equation in a weak coupling limit does not have any contribution from the off-diagonal elements of the density matrix. Therefore, our master equation is reduced to the classical rate equation [42]. To extract the off-diagonal contributions, we may consider a strong coupling regime or more complicated model such as a three level system. This will be our future work.
Although our system is equivalent to that analyzed in Ref. [54], there are various differences in the analysis between two papers. Indeed, Ref. [54] uses the discretized time evolution under the finite interval τ p /41 when the master equation is solved, we have obtained both numerical and analytic solutions under the continuous time evolution. Moreover, we have explicitly obtained the analytic form for the pumping current whose validity is quantitatively verified through the comparison with the numerical calculation. Furthermore, the discussion on the extended fluctuation theorems in Sec. 3.3 is completely new. Therefore, we believe that there exist several merits for the publication of our paper besides Ref. [54].
We have derived the analytic expression for the non-adiabatic pumping current J E na corrensponding to the geometrical phase in the adiabatic limit. We should note, however, that J E na is no longer geometric quantity as in the adiabatic case because the curvature depends on time. Such time dependent quantity may be interpreted by the method of Aharonov-Anandan phase [70].
We have also derived the fluctuation theorem with the temporal change of parameters in the case of τ p ω 0 → ∞ or N → ∞. This is a natural extension of the steady fluctuation theorem to the time dependent fluctuation theorems.
We expect that Floquet theory is applicable to our system because we study periodic modulations to the system. It is our future work to compare our analysis with that based on Floquet theory.
We only analyze the case of the continuous modulation of the temperature (14) under the assumption that both environments are always in equilibrium. It is straightforward to apply our formulation to Fermion systems such as the impurity Anderson model [51].
It should be noted that our theory cannot be applied to neither discontinuous changes of the temperatures nor fast modulation. Although it is possible to apply our formulation to non-Markovian processes, we are suspicious whether the analysis within this framework is meaningful, because non-Markovian processes may affect the state of the environments. We also note that the distinction between two current terms is no longer valid for non-Markovian processes, because λ(t) sometimes becomes positive positive values. We will discuss the non-Markovian pumping process elsewhere.

A. Cumulant generating function
In this Appendix, we briefly summarize the relationship between FCS and the cumulant generating function. Let us perform a projection measurement of Q at 0 and t, where their measeured values are set to be q 0 and q t , respectively. Here, we assume that Q satisfies [Q, ρ tot (0)] = 0, where ρ tot (t) is the total density matrix at time t without the counting field. The probability to measure q 0 and q t is given by where P q0 is the projection operator onto the eigenstates corresponding to the eigenvalue q 0 and U (t, 0) is the unitary time evolution operator of the total system, which is defined by is the total Hamiltonian. Thus, the probability of the current ∆q t at t is given by Let us introduce the characteristic function as From the identities P 2 q0 = P q0 , q 0 P q0 = QP q0 , q0 P q0 = 1, we can rewite (A4) as where U χ/2 (t, 0) is defined by and H tot,χ ≡ e iχQ H tot (t)e −iχQ . Equation (A5) automatically satisfies G(χ, 0) = 1 because of the relation U χ (t, 0) → 1 in the limit t → +0. Thus, the cumulant generating function satisfies S(χ, 0) = 0. Hence, all of the cumulants at t = 0 satisfy We introduce the total modulated densitiy matrix and the mudulated density matrix of the system is ρ(χ, t) ≡ Tr E ρ tot (χ, t).
Note that the argument presented here is still valid even for non-Markovian case. 13/26

B. Derivation of the quantum master equation with parameters modulation
In this appendix, we derive the FCS quantum master equation under the modulation of parameters though the derivation of the master equation without modulation is well known [28,59]. The total Hamiltonian is given by where H S is the Hamiltonian of the target system, H E is the Hamiltonian of environments, H SE is the interaction between the system and environments characterized by coupling constant g.
The total system with the full counting statistics is expressed by the modified von-Neumann equation from Eq. (A9) where the modified Liouvillian L χ (t) is defined by and The formal solution of Eq. (B3) can be written as To trace out the degree of freedom of environments, we introduce the Nakajima-Zwanzig projection operator where ρ eq E (π t ) is the equilibrium operator of environments and under the modulation of parameters π t such as temperatures and chemical potentials. We assume that environmets are always at equilibrium during the modulation of parameters π t . Even if we adopt such a simplification, the derivation of the master equation is non-trivial, because π t depends on time. The projection operators satisfy P(t)P(t ′ ) = P(t), P(t)Q(t ′ ) = 0, P(t)L 0 = L 0 P(t). From the relation Tr E (ρ eq E H SE ) = 0, we, thus, obtain Let us introduced the projected time evolution operator The time-evolutions of X and Y can be written as The formal solution of Eq. (B16) is given by where we have introduced Eq. (B17) can be rewritten as Substituting Eq. (B20) into Eq. (B15), we obtain d dt To perfome the perturvation for the small coupling constant g, we rewrite e −iL χ t and Eqs. (B18) as where we have introduced and L χ SE (τ ) = e iL0τ L χ SE e −iL0τ .

15/26
In the small g limit, Eqs. (B22) and (B23) are reduced to Thus, S(t) reduces to where we have used Eq. (B11) and defined M(t) as From the relation {1 − S(t)} −1 ≃ 1 + S(t) under the weak coupling limit, Eq. (B12) and Eq. (B21) can be rewritten as When we operate Eq. (B30) on ρ tot (χ, 0), we obtain the master equation where we have introduced and the initial correlation term which vanishes if ρ tot (χ, 0) = ρ(χ, 0)ρ eq E (π 0 ). There exist several characteristic time scales in Eq. (B31): the time scale τ S for the energy level of the system H S , the relaxation time scale τ R of the system, the correlation time scale τ C of the environments and the time scale Ω −1 of the modulation of parameters. τ C is the time scale that characterizes the symmetrized time correlation fuction Tr E [{B i (t) † , B j }ρ eq E ] where B i is the operator of environments when H SE can be expressed by H SE = i S i B i where S i is the operator of a system. 16/26 We apply the Markovian approximation τ C ≪ τ R , τ C ≪ Ω −1 to Eq. (B31). By using τ = τ C s, t = τ R u, s ∼ O(1), u ∼ O(1), the integration in Eq. (B31) become and m(t) is negligible because of, by using dP (τ )/dτ = ΩdP (v)/dv| v=Ωτ , becomes much smaller than the unity. Therefore, if we set the initial condition to ρ tot (χ, 0) = ρ(χ, 0)ρ eq E (π 0 ), we obtain the Markovian master equation

C. Adiabatic Markovian pumping: general expressions
In this section, we briefly review the adiabatic Markovian pumping process under the condition Ω/Γ ≪ 1. The argument in this section is parallel to that in Ref. [43]. Under this approximation, we can express the density matrix by the zero eigenvector which characterizes the steady state as |ρ(χ = 0, t) ≃ |λ 0 + (β(t)) where the subscript + represents the zero eigenvector and the upper script 0 represents the state without the counting field, i.e. χ = 0. Thus, the density matrix with the counting field χ can also be approximated by where we have introduced the proportional constant which satisfieṡ where we have used l χ + |λ χ + = 1. Note that l χ + | is reduced to 1| for χ = 0, which means trace.
Equation (C2) is readily solvable as In the second line we have introduced the total differentiation d. Substituting Eq.(C3) into Eq.(C1) we obtain the cumulant generating function S(χ, t) ≡ ln 1|ρ(χ, t) (see (A10)): 17/26 where the first, the second and the last terms on the right hand side (RHS), respectively, correspond to the geometrical phase, the dynamical phase and the surface term. Let us consider the energy transfer ∆q t from the right reservoir to the system during time t. The average of ∆q t can be calculated from the cumulant generating function as ∆q t = ∂S(χ, t)/∂(iχ)| χ=0 . Therefore, we obtain where ∆q g a represents the adiabatic pumping current in terms of the geometrical phase: and the ∆q t d is the adiabatic pumping current in terms of the dynamical phase: where ′ denotes the differentiation with respect to χ. Therefore, the adiabatic pumping current during the period τ p can be written as where we have introduced where S is the surface integral with the perimeter C and the integrand is called Berry curvarture. As shown in Appendix D, our adiabatic approximation is equivalent to that in Ref. [42] if we apply this formulation to the spin-boson system.

D. Adiabatic pumping for the spin-boson model
In this Appendix, we apply the general framework in the previous section to the spin-boson system (1) and (2) to verify whether we can reproduce the results in Ref. [42]. In this case Eq. (6) in Eq. (5) is given by Eqs. (21)- (24). Furthermore, we also introduce To avoid heavy notations, we replace the parameter dependence through β(t) by t.

18/26
From the differentiations of (25) and (27) we obtain . (D4) Substituting Eqs. (21), (24),(D1),(D2) into (D3) we can rewrite Subsituting this into Eq.(C9) we obtain the dynamical current which is equivalent to Eq.(13) in Ref. [42]. In the case of a symmetric junction under the environments Γ = Γ L = Γ R without average bias, Eq.(D6) can be rewritten as To derive the final equality of Eq. (D7) we use that n L and n R are sinusoidal functions of time, and thus, n L /(1 + n L + n R ) sweeps an identical area to that by n R /(1 + n L + n R ) during a period. On the other hand, let us rewrite the integrand in Eq.(C10) as Because of (D4) the only relevant term is the second component in the above equation. From the straightforward calculation, we can rewrite Eq.(D8) as Thus, we reproduce Eqs. (14) and (15) in Ref. [42].

F. Derivation of fluctuation theorems
In this Appendix, we explain the detailed derivation of two types of extended heat fluctuation theorems for slowly modulated temperatures. In Appendix F1, we discuss the extended heat fluctuation theorem in the limit τ p ω 0 → ∞. In Appendix F2, we discuss the extended fluctuation theorem for large amount of the transferred energy.