Relaxation Dynamics of Meso-Reservoirs

We study the phenomenology of maximum-entropy meso-reservoirs, where we assume that their local thermal equilibrium state changes consistently with the heat transferred between the meso-reservoirs. Depending on heat and matter carrying capacities, the chemical potentials and temperatures are allowed to vary in time, and using global conservation relations we solve their evolution equations. We compare two-terminal transport between bosonic and fermionic meso-reservoirs via systems that tightly couple energy and matter currents and systems that do not. For bosonic reservoirs we observe the temporary formation of a Bose-Einstein condensate in one of the meso-reservoirs from an initial nonequilibrium setup.


Introduction
Usually, a reservoir is treated as constant and inert to all systems that are coupled to it [1,2]. In contrast, realistic experimental implementations always deal with finite-sized reservoirs [3,4,5]. These are often too large to be treated exactly, but too small to neglect their dynamics in good faith, which has triggered research on thermodynamics with finite-size reservoirs [6,7,8]. Especially noteworthy in this context are experimental setups which utilize ultra cold atoms embedded in optomagnetical traps or optical latices [9,10,11,12,13].
In particular in nonequilibrium setups (e.g. realized by periodic driving or by locally different thermal states), also small systems may in the long-run transfer a significant amount of heat, and it may no longer be applicable to talk about a constant reservoir temperature or chemical potential [14,15,16,17,5,18]. For such reservoirs, we will use here the term meso-reservoir, with which we simply want to indicate that some sort of system back-action has to be taken into account, and that the state of a meso-reservoir is allowed to change in time.
We assume that the system of interest can transfer entropy to the meso-reservoir in form of heat (in this paper, we will only consider matter and energy currents). Moreover, we suppose that the meso-reservoir is subject to further processes that may potentially increase its entropy without additional heat transfer. One possible microscopic example for such a process are interactions with a larger super-reservoir that leave energy and particle number invariant. In Appendix A we discuss this in detail for an energy-conserving interaction. These processes generally lead in the eigenbasis of Hamiltonian H and number operator N to a fast decay of off-diagonal matrix elements (pure dephasing [19,20]), while the diagonals are by construction constant. Since for any density matrix ρ, the entropy associated to the diagonal elements only, S D = − i ρ ii ln ρ ii is larger than the entropy S = −Tr {ρ ln ρ} of the density matrix S D ≥ S [21], pure dephasing processes will increase the entropy [22] without injecting additional heat into the meso-reservoir. Furthermore, it is well-known that almost all states appear thermal when sufficiently many degrees of freedom are traced out, a statement known as canonical typicality [23,24,25,26,27,28]. In usual derivations of master equations [29] the reservoir is therefore always assumed in thermal equilibrium. From the perspective of the meso-reservoir, the presence of the system will induce transfers between diagonal elements of the density matrix together with the exchange of matter and energy. Additional elastic scattering processes within the mesoreservoir may support these equilibration processes. It should be noted that this will always also generate off-diagonal matrix elements in the meso-reservoir density matrix, see Appendix B. We assume that these are quickly dephased.
In consequence, we do here as usual assume that the meso-reservoir is always kept at local equilibrium, i.e., at its maximum entropy state where β(t) = 1/T (t) and µ(t) are now however time-dependent inverse temperature and chemical potential of the meso-reservoir, and H M R and N M R denote Hamiltonian and particle number operator of the meso-reservoir. At this maximum entropy state, the internal entropy production of the meso-reservoir vanishes [30], and the change of its entropy is solely governed by the heat transferṠ →Ṡ D = β(t)[J E − µ(t)J M ], quantified by the energy current J E and matter current J M entering the meso-reservoir via the system. These energy and matter currents can be quantified for a large number of models [31]. We note that energy contained in the interaction e.g. between system and meso-reservoir may in principle also affect its energy balance, but in the framework of our weak-coupling scenario we neglect these contributions in the long-term limit.
In this paper, we will consider the induced change of the meso-reservoir, which we compute self-consistently from the currents through the system. The system will only provide the dependence of the currents on temperatures and chemical potentials. Therefore, we implicitly assume that the fastest timescale is the equilibration of the meso-reservoir to a thermal state (1). Mainly for simplicity, we will also assume that the system quickly relaxes to its (possibly non-thermal) steady state, such that the current through the system has no signature of its initial state. In this paper, we will be interested in the slow changes of the meso-reservoir parameters β(t) and µ(t).
This paper is organized as follows. In Sec. 2 we derive the differential equations determining the evolution of temperatures and chemical potentials in a general way. In Sec. 3 we make these findings explicit for two fermionic meso-reservoirs coupled by single quantum dots. In Sec. 4 we show how to treat bosonic meso-reservoirs including the possibility of Bose-Einstein condensation. Finally, in Sec. 5 we compare efficiencies of converting thermal gradients to chemical work.

Consistent equilibrium states
Together with Eq. (1), the basic assumption of our framework is that particle number and energy content of meso-reservoirs α ∈ {L, R} (left, right) are given by integrals over densities of states D α (ω) and occupation numbers n α (ω) -supplemented by a few states that are separately treated, e.g. the ground state (thereby complementing previous work [32]) Here, the density of states depends on dimensionality and character of the meso-reservoir -but not on its thermodynamic parameters µ α and β α = 1/T α (we omit the explicit notion of time-dependence for brevity). In contrast, the occupation number depends explicitly on these where n + α corresponds to the Fermi-Dirac distribution in the fermionic and n − α to the Bose-Einstein distribution the bosonic case. Furthermore, we note that for bosons we have µ α < 0 and denotes the occupation of the ground state and thus allows the treatment of Bose-Einstein condensation [33,34]. For fermions, we set n g α = 0 (including it however would not substantially change the dynamics, since the Fermi-Dirac distribution is bounded). The change of particle numbers and energy content in every meso-reservoir has to balance the currents into each meso-reservoir, which gives -when the currents are known -rise to implicit ordinary differential equations for the thermodynamic potentials where C α is a 2 × 2 capacity matrix for reservoir α. It can be split into a continuum contribution and a ground state contribution A direct observation is that in contrast to classical quantities such as the geometric charge capacitance or heat capacity, the capacity matrix combines both temperatures and potentials to currents. Its matrix elements will in general depend on temperatures and potentials, too, which e.g. is not the case for the geometric charge capacitance. The integral term of the 11-component is a continuum version of what is usually called quantum capacitance [35]. It implicitly depends on the geometry via the chosen density of states, and is expected to approach the conventional geometric capacitance in the limit of large N α . In particular the separate treatment of the ground state in our case may however retain quantum features also in the macroscopic limit. After fixing the density of states one can explicitly calculate the capacity matrix C α -for which it is helpful to realize that the derivatives can be written -and the differential equation system (5) can be made explicit by inverting C α . Since the inverse of C α also depends on µ α , T α and the specifics of D α (ω), one thereby obtains a coupled set of highly nonlinear differential equations.
When one now considers two meso-reservoirs α ∈ {L, R} coupled indirectly via the same system, such a two-terminal transport setup will obey conservation of total energy and matter which will lead to two conservation laws. These imply that we can, in principle, eliminate two of the four thermodynamic variables to obtain two coupled nonlinear differential equations for e.g. the potential differences V = µ L − µ R and temperature differences . The shown tunneling rates Γ i α enable for the exchange of energy and particles. The individual currents through the two channels obey tight-coupling conditions, i.e., their matter and energy currents are proportional, but the combined current is not, J E ∝ J M (unless ε 1 = ε 2 ). Together with the entropy increase due to pure dephasing (wavy lines, see Appendix A), these processes are assumed to lead to fast local equilibration of the meso-reservoirs.
Since the conserved quantities may be quite complex, we have however technically found it more convenient to evolve all four variables according to and use the conservation laws as a numerical sanity check instead. We see immediately that at configurations with vanishing currents the chemical potentials and temperatures will remain stationary. Normally, this can only be fulfilled at global equilibrium (µ L = µ R and T L = T R ), as the vanishing of both currents imposes two independent conditions. However, in the tight-coupling regime (J E = εJ M ), these conditions are not independent, and in consequence, stationary states may arise that are not global equilibrium states. We note further that these equations could in principle be further simplified in the linear-response regime, where the currents are linear in potential and temperature differences [36,37]. However, also far away from this equilibrium regime, the currents must obey the second law (recall that J M and J E count positive when entering the right meso-reservoir), stating that the entropy production of the systeṁ is non-negative.
In the following, we will make this explicit for fermionic and bosonic reservoirs coupled via simple model systems, where we will assume that the energy and matter currents are in general not tightly coupled J E ∝ J M . This is rather generic for realistic systems, but to keep the analysis simple, we consider coupling the meso-reservoirs via two non-interacting systems that -when considered separately -exhibit tight coupling [38] (see Fig. 1). In conventional master equation derivations [29], our results can be obtained by performing the usual Born approximation for the full density matrix as ρ(t) = ρ S (t) ⊗ ρ M R (t) with Eq. (1), and for the Markov approximation assuming that ρ M R (t) changes even slower than the system density matrix ρ S (t), see Appendix C.

Fermionic transport
For fermionic meso-reservoirs, we can relate the thermodynamic potentials with the currents using three standard integrals In many solid state models one usually has only positive single particle energies (D α (ω < 0) = 0), and the integrals can be evaluated in this case too. However, to illustrate the method we consider the simpler case of the complete wideband limit D α (ω) = D α (normally corresponding e.g. to a 2d free electron gas) also for negative frequencies. Then, it is straightforward to show that the integrals become I = π 2 2 D α T 3 α (the same results would follow when only positive frequencies were allowed and the additional constraint µ α ≫ k B T α was imposed). Then, we obtain the capacity matrix which can easily be inverted. We see that in this particular case (negative energies or only positive energies with µ α ≫ k B T α ) the 11-component does not depend on temperatures or potentials -just as the geometric capacitance. Furthermore, the 22component is linear in the temperature, which is well-known for the electronic heat capacity. We also note that the capacity matrix becomes singular at zero temperature, the positivity of the system's entropy production (9) however ensures that the heat flow into a low temperature reservoir is always non-negative, and therefore the extreme zero-temperature limit cannot actually be reached. Finally, we note that matter and energy conservation imply the conserved quantities

Quantum-Dot Coupling
Our simplest example is the single-electron transistor in weak-coupling approximation.
Here, the current triggered by a single quantum dot hosting at most one electron is given by [31] where the constant γ depends on the details of the coupling between system and meso-reservoir and ε denotes the dot level. These expressions also arise from the Landauer current formula [39] when considering a strongly peaked transmission function. Obviously, the currents will vanish when ∆T = 0 and V = 0, but for the specific example it is also possible to obtain a vanishing current whenever The stationary state will therefore when plotted in the V − ∆T -plane depend on the initial condition.
When we consider two quantum dots with on-site energies ε i that connect the mesoreservoirs in parallel but do not interact directly as sketched in Fig. 1, the individual currents just add and we see that the tight-coupling condition is not obeyed, i.e., J E ∝ J M , when ε 1 = ε 2 . As before, the constants are given by the coupling details between system and meso-reservoirs (compare Fig. 1). Since each quantum dot can host at most one electron, the currents remain finite at infinite external bias, (where n + L (ε i ) → +1 and n + R (ε i ) → 0). In the loose coupling regime, the currents will in general only vanish when all thermodynamic parameters are equal, i.e., when T L = T R and µ L = µ R .

Meso-Reservoir Dynamics
In Fig. 2 we show the relaxation dynamics of the temperatures (dashed curves) and chemical potentials (solid curves) for two fermionic reservoirs in the wideband limit coupled via two non-interacting quantum dots. Whereas for the tight-coupling configuration ε 1 = ε 2 the stationary state of the complete system is a non-thermal nonequilibrium steady state (thin curves in lighter colors), the generic situation without tight-coupling (thick curves) yields an equilibrium state with T L = T R and µ L = µ R . However, for a near-tight-coupling configuration ε 1 ≈ ε 2 , one observes an intermediate pseudo-steady state before relaxation to complete equilibrium sets in at a much later time. The lifetime of the pseudo-steady state increases as one approaches the tightcoupling configuration, e.g., the smaller the difference |ε 1 − ε 2 | becomes. During the evolution into this pseudo-steady state, the device has created a potential difference, thereby performing chemical work. This is only possible with an initial temperature temperatures (dashed) of the hot (red) and cold (blue) reservoirs for one fermionic transport channel with energy ε 1 = ε (thin curves in lighter colors) and for two fermionic transport channels with energies ε 1 = ε and ε 2 = 1.1ε (thick curves). At tight coupling (thin curves in lighter colors), a stationary nonequilibrium state is found. Choosing the dot energies as different destroys the tight-coupling condition and leads to long-term equilibration (thick curves), but for intermediate times one observes a nonequilibrium pseudo-steady state where a potential bias is built up using the thermal gradient. The inset shows the time evolution of the particle number (solid) and internal energy (dashed) of the reservoirs for two transport channels. Other parameters: γ 1 = γ 2 = γ/2, µ 0 L = µ 0 R = −ε, T 0 L = 2ε, T 0 R = 0.24ε, capacity coefficients adjusted to εD L = εD R = 10000 such that initially N 0 L ≈ 9482 and N 0 R ≈ 37.
difference between the meso-reservoirs, and we will consider the energetic efficiency of this process in Sec. 5.

Bosonic transport
For bosons, we have to take into account some important differences: First, the singleparticle energies of the Hamiltonian must all be positive to bound the spectrum of the Hamiltonian. Second, the chemical potentials must be negative to bound the occupation of the individual modes. Finally, we allow for the possibility of a macroscopic occupation of the ground state n g α = e −βαµα − 1 −1 , which however does not significantly contribute to the total energy, cf. Eq. (2). Then, we can relate the currents with the change of the thermodynamic parameters using just three integrals where the integrals are given by When the chemical potentials are negative one can obtain under additional assumptions on the density of states analytic expressions for these integrals. We have also considered bosonic transport for a flat density of states (corresponding to 2d massive bosons, not shown), but here Bose-Einstein condensation will not occur. Therefore, we consider an ohmic density of states instead. With D α (ω) = J α ω (2d massless Bose gas supporting condensation [40]), the integrals become where Li n (z) = ∞ k=1 z k /k n denotes the polylog function [41]. In the high-temperature limit, the 22-component of the capacity matrix simplifies to C 22 α → 6J α ζ(3)T 2 α + (−µ α T α )J α π 2 /3, with Riemann Zeta-function ζ(3) and has thus simple linear and quadratic contributions in the temperature. For bosonic transport, the conserved quantities are given by We define condensation by assuming that half of all meso-reservoir particles are in the ground state at negligible chemical potential, which defines with Li 2 (1) = π 2 /6 the condensation temperature Since N α (t) is a dynamic variable, this also transfers to the condensation temperature, such that one may also consider the condensate fraction (number of particles in the ground state of the reservoir versus total number of particles in each meso-reservoir) instead.

Boson-Boson Transport model
When the two bosonic meso-reservoirs are coupled via non-interacting harmonic oscillators, the master equation currents can be written as which, similar to the fermionic case, can alternatively be obtained from the Landauer formula for heat transport [42,43] in case of a strongly peaked transmission function. In contrast to the fermionic case however, we observe that an infinite thermal bias (e.g. n − L (ε i ) → ∞ and n − R (ε i ) → 0) will let the currents diverge, too. This essentially arises since the carrying capacity of the system between the reservoirs is not bounded.  for the hot (red) and cold (blue) meso-reservoirs versus dimensionless time for two transport channels with energies ε 1 = ε and ε 2 = 1.1ε. The inset shows the corresponding number of particles (solid) and the internal energy (dashed) of the reservoirs. Other parameters: 5ε, and capacity coefficients ε 2 J L = ε 2 J R = 1000, such that initially N 0 L ≈ 577829 and N 0 R ≈ 35.

Meso-Reservoir Dynamics
Similar to the fermionic situation, the tight-coupling scenario will lead to a stationary non-thermal steady state (not shown). However, for slight modifications of the tightcoupling scenario, an intermediate nonequilibrium state will emerge with a lifetime defined by the deviation from tight coupling (see Fig. 3). Initially starting with a hot reservoir filled with many particles (red) and a cold reservoir with just a few particles (blue), one clearly observes that the initial thermal and particle gradients are used to dynamically induce condensation in the cold reservoir. Eventually, the condensate evaporates again and global equilibrium is reached.
To evaluate the quality of the induced condensate, we have also investigated the condensate fraction for different transport channel configurations in Fig. 4. For the case of a near tight-coupling configuration (solid) we observe a high quality condensate with about 80% of the particles occupying the ground state. This effect occurs due to the circumstance that the density in the cold reservoir grows faster than its temperature such that the condensation temperature is increased (inset) and Bose-Einstein condensation eventually sets in. Further away from the tight-coupling configuration (dashed) the condensate quality as well as its lifetime is reduced. Additionally, considering a near tight-coupling configuration with increased dot energies (dotted), we find that the condensate quality is further reduced, however it persists over longer times. We observe a temporarily decrease below a critical temperature leading to a macroscopic occupation of the ground-state energy level in the respective reservoir. Other parameters: 5ε, and capacity coefficients ε 2 J L = ε 2 J R = 1000, such that initially N 0 L ≈ 577829 and N 0 R ≈ 35. , and decay to zero for large times. In contrast, the cumulative efficiencies (dashed) may remain finite. We observe a temporarily negative chemical power output resulting from an inversion of the matter current direction. For the fermionic reservoirs we set T 0 L = 84ε, T 0 R = 0.24ε, and capacity coefficients εD L = εD R = 10000, such that initially N 0 L ≈ 577259 and N 0 R ≈ 37. For the bosonic reservoirs we set T 0 L = 20ε, T 0 R = 0.5ε, and capacity coefficients ε 2 J L = ε 2 J R = 1000, such that initially N 0 L ≈ 577829 and N 0 R ≈ 35. Other parameters: γ 1 = γ 2 = γ/2 and µ 0 L = µ 0 R = −ε.

Efficiency
In analogy to the intensively studied electronic solid-state setups, the transport setups suggested within this paper might be put to use as thermo-electric or thermo-chemical generators [10,16,44]. As a useful measure for the quality of such devices we consider the efficiency with which they generate power from an incoming heat current. The internal energy of meso-reservoir α changes according to the fundamental equation (we have no volume change in the reservoirs) dE α = T α dS α + µ α dN α . Here, the term T α dS α corresponds to the heat flow into the meso-reservoir, and the term µ α dN α represents the chemical work [45]. To define an energetic time-local efficiency, we consider the chemical power instead P = J M (µ R − µ L ). When the current flows from left to right J M > 0 although µ R > µ L , the power becomes positive, and the corresponding energetic efficiency is obtained by dividing by the heat flow entering the system from the hot (left) reservoir, i.e., for T L > T R and µ R > µ L one has for the efficiency [46] η where the bound by the time-dependent Carnot efficiency follows from the second law (9). It is actually only reached in the tight-coupling case [38,31] (not shown).
In contrast, when one considers the cumulative efficiency, defined as ratio of total chemical work performed and total heat influx from the hot (left) meso-reservoir up to time t upper bound is given by the initial Carnot efficiency η(t) ≤ η CA (0). In Fig. 5, we show the resulting efficiencies for a fermionic setup in the wideband limit (brown) and for a bosonic setup with an ohmic density of states (green). In both cases, we observe an increase of the time-local efficiencies (solid) with time, getting ever closer to the respective Carnot efficiencies (dotted-dashed). However, at some specific time (around tγ ≈ 10 6 ) the power output becomes negative and, hence, the time-local efficiencies vanish. Only for the bosonic setup this behavior is reversed for even later times, leading to a finite time-local efficiency again. Contrary, the cumulative efficiencies (dashed) are finite over a rather large time interval, and, moreover, they can have finite values even for arbitrary long times as can be seen for the bosonic setup.

Summary
We have demonstrated that with a simple phenomenological approach conservation laws may be used to track the dynamical evolution of thermodynamic parameters of mesoreservoirs. Our approach is applicable to a rather wide range of models, although we have exemplified it only for two-terminal fermionic and bosonic transport setups and although we have for simplicity neglected the energy and particle content of the system. The generalization to other systems is straightforward, it is however necessary that the currents through the system obey the first law (conservation of matter and energy currents) and the second law (to prevent unphysical temperatures and chemical potentials). Naturally, when considering reservoirs of infinite capacities (formally e.g. by considering the limit D α , J α → ∞), temperatures and chemical potentials remain fixed and we recover the usual weak-coupling master equation results.
In general, we have observed global equilibration of both meso-reservoirs in the long-term limit, except for the highly idealized tight-coupling scenario, where the nonequilibrium stationary state is frozen due to vanishing currents. For situations close to tight-coupling, the system assumes a temporary pseudo-steady-state, and the dwell time of the system in this nonequilibrium state roughly depends on the deviation from the tight-coupling scenario. We note that the dynamical generation of such nonequilibrium pseudo-steady state may be desirable in many experimental contexts, and we have sketched the efficient use of such phases as a thermo-electric generator and for preparing a Bose-Einstein condensate.
where ρ denotes the full density matrix. We choose to represent its most general form by energy eigenstates of A where we use the convention that Tr ρ (B,ij) = 1 for all i and j, such that in the energy eigenbasis of A we have the representation ρ A = ij ρ ij A |i j|. Inserting this decomposition and also the decomposition of the interaction Hamiltonian (which now need not commute with H A ), we obtain For the dynamics of the reduced density matrix elements this implieṡ We note that this equation is non-perturbative in the interaction. While it is probably useless for practical calculations, one can see that there is no direct coupling between different diagonal elements. This implies that to transfer population between different diagonal elements ρ ii A and ρ jj A one has to populate also off-diagonal elements ρ ij A as an intermediate step, too.
The basic assumption behind our Eq. (1) is that -with A describing the mesoreservoir and B taking the role of the system -additional pure dephasing processes as described in Appendix A quickly eliminate the off-diagonal matrix-elements in mesoreservoirs after local equilibration has been reached. In contrast to usual derivations of master equations, we do however not neglect the energy (and particles) injected into the meso-reservoir.

Appendix C. Derivation of a master equation
Here, we will follow the usual derivation of a master equation in the weak-coupling limit for time-dependent chemical potentials and temperatures. We will perform the derivation only for a single meso-reservoir B -onto which in absence of stationary transport a small system A would have negligible effect -and highlight the changes arising from its time-dependence. In the interaction picture (defined by bold-written operators) A(t) = e +i(H A +H B )t Ae −i(H A +H B )t , the complete density matrix follows the von-Neumann equationρ = −i[H AB (t), ρ(t)]. Re-inserting the formal solution in the right-hand side, one obtainṡ We now insert the Born approximation with a time-dependent reservoir density matrix ρ(t) = ρ A (t) ⊗ ρ B (t) and trace out the reservoir degrees of freedom ρ A (t) = Tr B {ρ(t)}. We note that due to trace conservation we have Tr B {ρ B } = 0. Furthermore, we assume that Tr B {B α ρ 0 B } = 0, which is fulfilled for many microscopic models from the start but can always be achieved by a suitable transformation. Then we can insert the decomposition of the interaction Hamiltonian H AB (t) = α A α (t) ⊗ B α (t) to obtain an integro-differential equation (non-Markovian master equation) for the system density matrixρ Next, we use the invariance of the trace under cyclic permutations and introduce the reservoir correlation function C αβ (τ, t ′ ) = Tr B {B α (τ )B β ρ B (t ′ )}. This requires to make use of [H B , ρ B (t)] = 0, cf. Eq. (1). After the substitution τ = t−t ′ , the master equation becomesρ We apply the Markov approximation by assuming that the reservoir correlation function decays with respect to its first argument much faster than ρ A (t − τ ) changes. In fact, one can for many microscopic models explicitly confirm that the correlation function has a Dirac-δ-function-type behavior near τ = 0. Since ρ B (t − τ ) changes even slower (the time-dependence in the second argument of the correlation function only refers to the change in temperatures and chemical potentials), this allows to replace ρ A (t − τ ) → ρ A (t) and C αβ (±τ, t − τ ) → C αβ (±τ, t) under the integral and to extend its upper bound to infinity, yielding a Markovian master equatioṅ Finally, we represent the system coupling operators in terms of eigenvectors of the system Hamiltonian A α (t) = ij A ij α e +i(E i −E j )t |i j| and neglect for large times all terms that oscillate in t (secular approximation), i.e., exp[i( ∞ 0 dτ e +i(E i −E j )τ C αβ (+τ, t) .
(C. 5) To see that this master equation is of Lindblad type, we can insert the even γ αβ (ω, t) = C αβ (τ, t)e +iωτ dτ and odd σ αβ (ω, t) = C αβ (τ, t)sgn(τ )e +iωτ dτ Fourier transforms of the reservoir correlation functions with respect to their first argument with which we can replace the half-sided Fourier transforms to yielḋ This is exactly the same Lindblad master equation as one would have obtained when assuming a constant reservoir and afterwards inserting the time-dependent reservoir parameters [31]. The used approximations do not go beyond those normally used in the derivation of master equation, except that some back-action onto the reservoir is taken into account.