Slow dynamics and thermodynamics of open quantum systems

TL;DR

This paper introduces a perturbation theory for open quantum systems with slowly changing parameters, enabling analysis of finite-time thermodynamic processes and deriving formulas for efficiency at maximum power.

Contribution

It develops a perturbation framework for master equations with slowly varying parameters, applicable to finite-time thermodynamics and efficiency calculations.

Findings

Derived a general formula for efficiency at maximum power of finite-time Carnot engines.

Clarified conditions under which Curzon-Ahlborn efficiency is recovered.

Provided a method to evaluate work and heat perturbatively in slowly driven systems.

Abstract

We develop a perturbation theory of quantum (and classical) master equations with slowly varying parameters, applicable to systems which are externally controlled on a time scale much longer than their characteristic relaxation time. We apply this technique to the analysis of finite-time isothermal processes in which, differently from quasi-static transformations, the state of the system is not able to continuously relax to the equilibrium ensemble. Our approach allows to formally evaluate perturbations up to arbitrary order to the work and heat exchange associated to an arbitrary process. Within first order in the perturbation expansion, we identify a general formula for the efficiency at maximum power of a finite-time Carnot engine. We also clarify under which assumptions and in which limit one can recover previous phenomenological results as, for example, the Curzon-Ahlborn efficiency.

Figures

Equations

$$\dot{\rho }(t)={\mathcal{L}}_t[\rho (t)].$$

$$Assumption1:$$

$$\rho (t)\simeq {\rho }_0(t).$$

$${\tilde{\mathcal{L}}}_{t^{\prime }}$$

$$\dot{\tilde{\rho }}(t^{\prime })=\tau {\tilde{\mathcal{L}}}_{t^{\prime }}[\tilde{\rho }(t^{\prime })],$$

$$\tilde{\rho }(t^{\prime })={\tilde{\rho }}_0(t^{\prime })+\frac{{\tilde{\rho }}_1(t^{\prime })}{\tau }+\frac{{\tilde{\rho }}_2(t^{\prime })}{{\tau }^2}+\dots .$$

$$\mathrm{tr}[{\tilde{\rho }}_j(t^{\prime })]=0,\forall j>0.$$

$$0$$

$${\dot{\tilde{\rho }}}_j(t^{\prime })$$

$${\tilde{\rho }}_j(t^{\prime })={\left[({\tilde{\mathcal{L}}}_{t^{\prime }}\mathcal{P}{)}^{-1}\frac{d}{d{t}^{\prime }}\right]}^j{\tilde{\rho }}_0(t^{\prime }).$$

$$\rho (t)=\frac{1}{1-({\mathcal{L}}_t\mathcal{P}{)}^{-1}\frac{d}{dt}}{\rho }_0(t).$$

$$Q={\int }_0^{\tau }\mathrm{Tr}[\dot{\rho }(t)H(t)]dt={\int }_0^1\mathrm{Tr}[\dot{\tilde{\rho }}(t^{\prime })\tilde{H}(t^{\prime })]d{t}^{\prime },$$

$$W={\int }_0^{\tau }\mathrm{Tr}[\rho (t)\dot{H}(t)]dt={\int }_0^1\mathrm{Tr}[\tilde{\rho }(t^{\prime })\dot{\tilde{H}}(t^{\prime })]d{t}^{\prime },$$

$$U_1(t)$$

$$S_1(t)$$

$$Q_1$$

$$W_1$$

$$Assumption3:$$

$$P$$

$$\eta$$

$${\eta }_{\mathrm{C}}=1+Q_0^C/Q_0^H=1-T_C/T_H.$$

$${\eta }^{\ast }={\left[\frac{2}{{\eta }_{\mathrm{C}}}-\frac{1}{1+\sqrt{Q_1^C/Q_1^H}}\right]}^{-1}.$$

$$Q_1^C/Q_1^H=(T_C/T_H{)}^{1-\alpha },$$

$${\eta }^{\ast }={\left[\frac{2}{1-T_C/T_H}-\frac{1}{1+(T_C/T_H{)}^{(1-\alpha )/2}}\right]}^{-1}.$$

$${\eta }^{\ast }={\eta }_{\mathrm{C}}/2+{{\eta }_{\mathrm{C}}}^2/8+{{\eta }_{\mathrm{C}}}^3(2-\alpha )/32+O({{\eta }_{\mathrm{C}}}^4).$$

$$H=\frac{\hslash \omega }{2}{\sigma }_z+V(t),$$

$$\dot{\rho }(t)$$

$$\displaystyle+\gamma(N+1)\left(\sigma_{-}\rho(t)\sigma_{+}-\frac{1}{2}\big{\{}\sigma_{+}\sigma_{-},\rho(t)\big{\}}\right)$$

$$\displaystyle+\gamma N\left(\sigma_{+}\rho(t)\sigma_{-}-\frac{1}{2}\big{\{}\sigma_{-}\sigma_{+},\rho(t)\big{\}}\right),$$

$$\gamma$$

$$N$$

$$V_I=\frac{\hslash \Delta }{2}{\sigma }_x,$$

$$\rho=\left[\begin{array}[]{cc}\rho_{11}&\rho_{12}\\ \rho_{21}&\rho_{22}\end{array}\right]\rightarrow||\rho\rangle=[\rho_{11},\rho_{12},\rho_{21},\rho_{22}]^{\top},$$

$$P=\frac{1}{2}\left[\begin{array}[]{cccc}1&0&0&-1\\ 0&0&0&0\\ 0&0&0&0\\ -1&0&0&1\end{array}\right],$$

Full text

A central result in the study of open quantum systems Breuer2002 ; Rivas2012 is the Markovian Master Equation (MME) approach which, under realistic assumptions, describes the temporal evolution of a system of interest $A$ induced by a weak coupling with a large external environment $E$. This consists in a first order linear differential equation $\dot{\rho}(t)=\mathcal{L}[{\rho}(t)]$ where $\rho(t)$ is the density matrix of $A$ and where the generator of the dynamics is provided by a quantum Liouvillian superoperator $\mathcal{L}$ that can be casted in the so called Gorini-Kossakowski-Sudarshan-Lindblad form KOS ; LIN ; GO . For autonomous systems the latter does not exhibit an explicit time dependence and the dynamics of $A$ exponentially relaxes to a (typically unique) equilibrium steady state $\rho_{0}$ identified by the null eigenvector equation $\mathcal{L}[{\rho}_{0}]=0$. MMEs can also be employed to describe the temporal evolution of $A$ when it is tampered by the presence of slow varying, external driving forces. Indeed, as long as these operate on a time scale which is much larger than the characteristic bath correlations times and the inverse frequencies of the system of interest, the effective coupling between $A$ and $E$ adapts instantaneously to the driving control, resulting on a MME governed by a time-dependent Liouvillian generator, i.e.

[TABLE]

An explicit integration of this equation is in general difficult to obtain. Yet, if the control forces are so slow that their associated time scale is also larger with respect to the relaxation time of the system induced by the interaction with $E$, one expects $\rho(t)$ to approximately follow the instantaneous equilibrium state $\rho_{0}(t)$ that nullifies ${\mathcal{L}}_{t}$. Our aim is to estimate quantitative deviations from this ultra-slow driving regime. For this purpose we develop a perturbation theory valid in the limit of slowly varying Liouvillians ${\cal L}_{t}$ and derive a formal solution of Eq. (1) which allows one to evaluate non-equilibrium corrections up to arbitrary order.

The main motivation of our analysis is to model thermodynamic processes and cycles beyond the usual reversible limit which is strictly valid only for infinitely long quasi-static transformations. Finite-time thermodynamics Andresen2011 ; Benenti2016 is a well established research field which is focused on this issue and in particular on the tradeoff between efficiency and power of realistic heat engines. Several results in this context have been derived from the geometrical notion of thermodynamic length Salamon1983 , from non-equilibrium identities known as fluctuation theorems Campisi2011 ; Cavina2016 , or from phenomenological models of heat engines Benenti2016 ; Curzon1975 . The latter approach led to the identification of a quite general value for the efficiency at maximum power which is the celebrated Curzon-Ahlborn (Chambadal-Novikov) efficiency Curzon1975 ; Chambadal1957 ; Novikov1957 . Our framework is complementary to previous approaches since it allows to explicitly express irreversible thermodynamic quantities (e.g. heat and work) in terms of the Liouvillian operator that governs the system dynamics. In this way we identify a general link between the frequency scaling of the spectral density and the efficiency of finite-time Carnot heat engines, clarifying for which kind of thermal baths the Curzon-Ahlborn result or other particular limits can be recovered.

Similar questions and problems have been addressed in the literature with different aims and methods. Finite-time quantum thermodynamics Alicki1979 ; Kosloff2013 ; Brandner2016 and Brownian quantum engines Schmiedl2007 ; Dechant2016 were studied using the formalism of open quantum systems. In particular single-qubit heat engines subject to Markovian dissipation were considered in Szczygielski2013 ; Esposito2009 ; Esposito2010bis ; Wu2006 . The impact of the bath spectral density on the efficiency of quantum engines was also noticed in the context of autonomous heat pumps Correa2014 ; Palao2016 and single-qubit minimal machines Gelbwaser2013 . Universal features and bounds for the efficiency at maximum power of finite-time Carnot cycles were identified in Esposito2010 ; Schmiedl2007 ; Vandenbroeck2005 ; Johal2016 for generic heat engines. Other results were obtained combining MMEs and linear response theory Davies1978 ; Bonanca2014 ; Benenti2016 and similar approaches were used to demonstrate the universality of heat engines in the limit of infinitesimal cycles Uzdin2015 . Outside the field of thermodynamics, our theory of slowly driven open quantum systems also contributes to the current research activity on quantum adiabatic driving. Several generalizations of the adiabatic theorem to open quantum systems have been already proposed Davies1978 ; Carollo2006 ; Oreshkov2010 ; Albash2012 . Our aims and results are however different: we are not interested in the derivation of a modified MME for time dependent Hamiltonians but only on the dynamical evolution of $A$ for an assigned time dependent Liouvillian ${\cal L}_{t}$ that, for all $t$, admits a unique (instantaneous) equilibrium state. This kind of adiabatic approach has recently proved to be effective in the study of general time dependent Liouvillians with higher dimensional kernels Avron2012 ; Venuti2016 ; Albert2016 . We finally would like to stress that, by replacing ${\rho}(t)$ with a probability vector all the results that we are going to present in the following are directly applicable also to classical continuous-time Markov processes Norris1997 .

*Slow Driving Perturbation theory.— * Let us consider the case of an open quantum system $A$ evolving as in Eq. (1) under the action of a quantum Liuovillian ${\cal L}_{t}$ which exhibits an explicit temporal dependence induced by the external modulation of some control parameters, say the value of a magnetic field or the intensity of a laser which are gradually changed according to some assigned protocol. In what follows we shall assume that for all $t$, ${\cal L}_{t}$ admits a unique zero (instantaneous) eigenstate $\rho_{0}(t)$ and that all the other eigenvalues have a strictly negative real part (in this case the map is said to be relaxing or mixing Rivas2012 ; Burgarth2013 ). This causes the system to exponentially converge to the instantaneous steady state, for a fixed value of the external modulation:

[TABLE]

Under these conditions, one can easily verify that for an infinitely slow modulation of ${\cal L}_{t}$, the system $A$ will be forced to follow quasi-statically the trajectory determined by the time-dependent density matrix $\rho_{0}(t)$, i.e.

[TABLE]

Physically this follows from the fact that, in this regime, there is enough time for $A$ to track the instantaneous equilibrium states defined by ${\cal L}_{t}$. This solution well approximates the dynamics of realistic configurations where, for instance, the system of interest is in thermal contact with a reservoir while being subject to a quasi-static external control, continuously relaxing to the instantaneous Gibbs ensemble. Notice however that at this stage of the analysis $\rho_{0}(t)$ can be an arbitrary quantum state, covering more general open evolutions including also systems in contact with engineered non-thermal baths (e.g. squeezed environments, negative temperatures, artificial dissipative maps, etc.). To characterize deviations from the quasi-static solution Eq. (3) we find convenient to introduce the following time-rescaled quantities

[TABLE]

where $\tau$ is the duration of the protocol, i.e. the total time interval on which the system evolves under the influence of the external control. With this choice Eq. (1) can be expressed as

[TABLE]

the dynamics being confined now into the unit interval $t^{\prime}\in[0,1]$. In this way the total duration of the protocol appears only as a multiplicative factor while all the information about its “shape” is contained in $\tilde{\mathcal{L}}_{t^{\prime}}$. Notice also that while $\tilde{\mathcal{L}}_{t^{\prime}}$ is independent of $\tau$, the time-rescaled solution $\tilde{\rho}(t^{\prime})$ of Eq. (5) is not. In particular the quasi-static solution (3) can be recovered from Eq. (5) in the asymptotic limit $\tau\rightarrow\infty$. Therefore we look for a perturbation expansion of the solution of Eq. (5) in powers of $1/\tau$:

[TABLE]

where normalization implies that all perturbations are traceless

[TABLE]

A rigorous mathematical analysis of the convergence properties of the series (6) is beyond the aim of this work. For our purposes it is sufficient that Eq. (6), truncated up to a finite order, provides a good approximation of the dynamics. The validity of this approach is verified in several numerical examples presented in the Supplemental Material (Appendix A,B and C).

Substituting Eq. (6) in both sides of Eq. (5), and equating the terms proportional to the same powers of $1/\tau$, we get the following set of recursive relations:

[TABLE]

As expected, Eq. (8) implies that $\tilde{\rho}_{0}(t^{\prime})$ is the (time-rescaled) instantaneous steady state of the Liouvillian. Moreover from Eq. (9), all finite-time perturbations can be recursively obtained. Indeed, exploiting Eq. (7), Eq. (9) can be univocally inverted by introducing the projector $\mathcal{P}$ on the traceless subspace, i.e. $\mathcal{P}[X]=X-{\rm tr}(X)\mathbb{I}/d$ with $\mathbb{I}$ being the identity operator on $A$ and $d$ the dimension of its Hilbert space. Now since ${\rm tr}[\rho_{0}(t)]=1\neq 0$, the Liouvillian within the subspace of traceless operators is invertible and we also have $(\mathcal{P}\mathcal{L}_{t}\mathcal{P})^{-1}=(\mathcal{L}_{t}\mathcal{P})^{-1}$. Accordingly for all $j$ we can write Eq. (9) as $\tilde{\rho}_{j+1}(t^{\prime})=(\tilde{\mathcal{L}}_{t^{\prime}}\mathcal{P})^{-1}\dot{\tilde{\rho}}_{j}(t^{\prime})$, which by iteration yields explicit formulas for each perturbative term:

[TABLE]

Switching from time-rescaled variables back to the original ones, the full solution of the MME can be compactly expressed as an operator geometric series:

[TABLE]

It is worth stressing that the above solution is unique and independent on the initial conditions. This is due to the fact that, at this level of approximation, we are neglecting any term exponentially decaying in $\tau$. In other words, the perturbations that we are considering are intrinsic to the Liouvillian operator without any influence from the initial state. In practice this means that all exact solutions of the master equation, after an exponentially short transient depending on the initial conditions, will converge to the asymptotic solution (11). This is a fundamental difference with respect to the Hamiltonian adiabatic theorem in which finite-time corrections depend on initial conditions and decay exponentially in $\tau$. In the next sections we present a couple of relevant applications of the results presented above to the context of finite-time thermodynamics.

*Finite-time isothermal process.— * Let us consider the case of a thermal MMEs (1) describing the dynamical evolution of $A$ induced by an external driving that slowly modifies its Hamiltonian $H(t)$ while the system is constantly kept in thermal contact with a bath of fixed inverse temperature $\beta$. In this scenario the instantaneous equilibrium state of the problem can be identified with the Gibbs density matrix $\rho_{0}(t)=\exp[-\beta H(t)]/Z(t)$, with $Z(t)={\rm tr}\{\exp[-\beta H(t)]\}$ being the associated partition function. Exploiting the derivation of the previous section we wish to determine how departures from the associated quasi-static trajectory (3) influence the thermodynamic properties of the process. For this purpose we remind that the mean energy and the von Neumann entropy of $A$ are given by $U(t)={\rm tr}[H(t)\rho(t)]$ and $S(t)=-{\rm tr}[\rho(t){\rm log}(\rho(t))]$, respectively. Moreover, following a common approach Alicki1979 ; Kosloff2013 ; Anders2013 we identify the mean heat absorbed by the system during the time interval $[0,\tau]$ with

[TABLE]

and the mean work done on the system with

[TABLE]

such that the first law of thermodynamics is obtained as $\Delta U=U(\tau)-U(0)=W+Q$. Now, since the time-rescaled Hamiltonian $\tilde{H}(t^{\prime})$ does not depend on $\tau$ but only on the shape of the protocol, all the previous quantities are influenced by $\tau$ only through the solution of the master equation $\rho(t)$ for which we know how to evaluate each perturbative term of the series (6). Therefore we can write $X=X_{0}+X_{1}/\tau+X_{2}/\tau^{2}+\dots$ with $X=U,S,W,Q$, and we can easily evaluate each term using Eq. (10). For example, at the zeroth order approximation we recover the standard results of equilibrium thermodynamics: $U_{0}(t)=-\frac{\partial}{\partial\beta}\log Z(t)$, $S_{0}(t)=(1-\frac{\partial}{\partial\beta})\log Z(t)$, $Q_{0}=\Delta S_{0}/\beta$ and $W_{0}=\Delta U_{0}-\Delta S_{0}/\beta$. The corresponding first order irreversible corrections are instead

[TABLE]

with $\tilde{\rho}_{1}(t)$ as in (10). It is worth noticing that, independently from the selected form of the MME the first law of thermodynamics is valid at the level of each perturbative coefficient $\Delta U_{j}=W_{j}+Q_{j}$ and that, for an initial Gibbs state $\rho(0)=\rho_{0}(0)$, the second law can be expressed as $Q\leq Q_{0}$ or equivalently as $W\geq W_{0}$. Taking the limit of large $\tau$, this implies that $Q_{1}\leq 0$ and $W_{1}\geq 0$. Moreover, if we consider the time-reversed process $\tilde{\mathcal{L}}^{\leftarrow}_{t^{\prime}}=\tilde{\mathcal{L}}_{1-t^{\prime}}$, then it is easy to check that odd-order perturbations are invariant while even-ored perturbations change sign $Q_{j}^{\leftarrow}=(-1)^{j+1}Q_{j}$ and $W_{j}^{\leftarrow}=(-1)^{j+1}W_{j}$.

*Finite-time Carnot cycles.— * As a second application of our slow driving perturbative approach consider the case where $A$, initialized in a Gibbs state with Hamiltonian $H_{A}$, evolves following a Carnot cycle composed by:

Isothermal expansion: the system is put in contact with a hot bath of temperature $T_{H}$ and the Hamiltonian is slowly changed from $H^{A}$ to $H^{B}$, in a time interval $\tau_{H}$. 2. 2.

Adiabatic expansion: the Hamiltonian is suddenly changed from $H^{B}$ to $(T_{C}/T_{H})H^{B}$. 3. 3.

Isothermal compression: the system is put in contact with a cold bath of temperature $T_{C}$ and the Hamiltonian is slowly changed from $(T_{C}/T_{H})H_{B}$ to $(T_{C}/T_{H})H_{A}$, in a time interval $\tau_{C}$. 4. 4.

Adiabatic compression: the Hamiltonian is suddenly changed from $(T_{C}/T_{H})H_{A}$ back to $H_{A}$.

The scaling factor $T_{C}/T_{H}$ characterizing the adiabatic operations is chosen to ensure that the stationary (Gibbs) state $\rho_{0}(t)$ evolves continuously in time during the cycle such that, for infinitely long processes, the system remains always in equilibrium without irreversible jumps. In addition to this standard requirement we also assume that $\rho_{0}(t)$ is sufficiently smooth, so that all the derivatives appearing in our perturbation theory are well defined along the cycle. Specifically, since in what follows we are going to consider only first order perturbations, we assume:

*Assumption 2: $\rho_{0}(t)$ is continuous and differentiable.

We also assume that, apart from the Hamiltonian scaling factor $T_{C}/T_{H}$ and its time length, the cold isothermal protocol is the time-reversal of the hot one: $\tilde{H}^{C}(t^{\prime})=(T_{C}/T_{H})\tilde{H}^{H}(1-t^{\prime})$. In terms of the associated time-rescaled Gibbs states, this is equivalent to:

[TABLE]

The latter assumption is common in many realistic heat engines and can be relaxed if we are free to optimize the shape of the two isothermal processes. Indeed in this case, within first order in the perturbation expansion that we are going to present later, the maximum output power is obtained when the driving protocol respects the time-reversal symmetry condition (18).

Now we are ready to analyze the performances of our finite-time Carnot engine. In the limit of many cycles the system evolution becomes periodic ($\Delta U=0$) and the work per cycle depends only on the heat exchanged in the hot and cold isothermal processes: $W=-Q^{H}-Q^{C}$. The power is the ratio between the work extracted and the time length of the cycle $P=-W/(\tau_{H}+\tau_{C})=(Q^{H}+Q^{C})/(\tau_{H}+\tau_{C})$, while the thermodynamic efficiency is defined as $\eta=-W/Q^{H}=1+Q^{C}/Q^{H}$. If we substitute the perturbative expansion, keeping only first order terms in $\tau_{H}$ and $\tau_{C}$, we obtain

[TABLE]

In the quasi-static limit $\tau_{C},\tau_{H}\rightarrow\infty$, the power tends to zero and the efficiency reaches the Carnot limit

[TABLE]

The power can be maximized with respect to $\tau_{H}$ and $\tau_{C}$, taking into account the physical constraints $Q_{0}^{H}>0$, $Q_{0}^{C},Q_{1}^{H},Q_{1}^{C}<0$. The corresponding efficiency at maximum power can be easily computed (see e.g. Schmiedl2007 ):

[TABLE]

We also notice that, in the particular case in which heat corrections are proportional to the temperature, Eq. (22) reduces to the efficiency derived in Ref. Esposito2010 .

We can now make use of Eq. (16) and (10) to express both $Q_{1}^{H}$ and $Q_{1}^{C}$ in terms of the Liouvillians $\mathcal{L}_{t^{\prime}}^{H}$ and $\mathcal{L}_{t^{\prime}}^{C}$ describing the two isothermal processes. Both irreversible heat corrections depend on the particular choice of the protocol however we find that their ratio has the following scaling:

[TABLE]

where $\alpha$ is the frequency exponent of the bath spectral density $J(\omega)\propto\omega^{\alpha}$ which is assumed to be the same for both the hot and the cold reservoirs. The proof of Eq. (23) can be found in the Supplemental Material (Appendix D) and follows from the time-reversal condition (18) and a general scaling property characterizing all MME obtained from standard microscopic models. Substituting (23) in (22), we get our final result (see Fig. 1 for a plot):

[TABLE]

The peculiarity of the universal expression (24) is that, by changing the parameter $\alpha$, it interpolates between several results which were previously obtained in the literature. For example, for a flat spectral density ($\alpha=0$) we recover the Curzon-Ahlborn Novikov1957 ; Chambadal1957 ; Curzon1975 efficiency $\eta^{*}|_{\alpha=0}=1-\sqrt{T_{C}/T_{H}}$, while for an Ohmic bath $\alpha=1$ we get $\eta^{*}|_{\alpha=1}=2\eta_{\rm C}/(4-\eta_{\rm C})$ which is the efficiency obtained by Schmiedl and Seifert for a specific Brownian engine Schmiedl2007 . Finally, consistently with known lower and upper bounds Schmiedl2007 ; Esposito2010 ; Benenti2016 , by taking the two limits of infinitely super-Ohmic or sub-Ohmic bath, we get $\eta^{*}|_{\alpha\rightarrow\infty}=\eta_{\rm C}/2\leq\eta^{*}\leq\eta_{\rm C}/(2-\eta_{\rm C})=\eta^{*}|_{\alpha\rightarrow-\infty}$.

In Ref.s Gelbwaser2013 ; Correa2014 ; Palao2016 , the efficiency of minimal models of heat engines and refrigerators was also linked to the spectral density. Such models are quite different from the Carnot cycles considered in this work but still one can find a qualitative agreement with our results.

A final remark should be made on the range of applicability of Eq. (24). Its derivation follows from the first order expansion performed in Eq.s (19) and (20) which is not guaranteed to be always valid in the regime of maximum power, especially for $T_{C}/T_{H}\simeq 0$ where the optimal values of $\tau_{H}$ and $\tau_{C}$ could become too small. This is confirmed by the single-qubit exact simulation shown in Fig. 1 where we observe deviations from (24) for small values of $T_{C}/T_{H}$. It is interesting then to expand Eq. (24) around the opposite regime, i.e. for small values of $\eta_{\rm C}=1-T_{C}/T_{H}$, obtaining

[TABLE]

We notice that the first and second order coefficients (i.e. $1$ and $1/8$) are independent of $\alpha$ and correspond to the same values in the Taylor expansion of the Curzon-Ahlborn Curzon1975 ; Esposito2009 and of the Schmiedl-Seifert Schmiedl2007 efficiencies. This also implies that, up to second order in $\eta_{\rm C}$, our results are in agreement with previous analyses based on linear response thermodynamics Johal2016 ; Benenti2016 .

Conclusions and outlook.— We have derived a perturbation theory for the solution of generic master equations with slowly varying coefficients. We focused in particular on finite-time thermodynamic processes beyond the reversible limit. Our analysis allows to analytically derive finite-time thermodynamic quantities in terms of the Liouvillian operator. We have also shown that, for a Carnot cycle, the efficiency at maximum power can be reduced to a universal formula depending only on the temperature ratio and on the scaling of the bath spectral density. Other implications of the perturbation theory presented in this work could also be studied in the future, in particular within the general context of quantum adiabatic driving.

Acknowledgments.— We thank L. C. Venuti, M. Fraas, A. Carollo, V. V. Albert, L. A. Correa, D. Alonso and D. Gelbwaser, for fruitful discussions.

SUPPLEMENTAL MATERIAL

Appendix A Slow driving of a two-level system

In this section we compare the perturbation theory presented in the main text with the exact open dynamics of a two-level system. Consider a system with Hamiltonian

[TABLE]

where $\omega$ is the energy splitting, $\sigma_{z}$ is the Pauli matrix and $V(t)$ is a generic perturbation. Moving to interaction picture with respect to $\hbar\omega\sigma_{z}/2$, we can model the contact of the system with a bosonic heat bath through a standard master equation Breuer2002A ; Rivas2012A :

[TABLE]

where

[TABLE]

are the damping rate and the mean excitation number respectively, associated to a thermal bath with spectral density $J(\omega)\propto\omega^{\alpha}$ and inverse temperature $\beta$.

Now, just to make an example in which off-diagonal coherences are relevant in the dynamics, consider the following perturbation

[TABLE]

and for the moment assume that $\omega$ and $\Delta$ are constant parameters. Recasting the $2\times 2$ density matrix into a $4\times 1$ complex vector

[TABLE]

it is easy to check that the projector on the traceless subspace $\mathcal{P}$ is represented by the matrix

[TABLE]

while the Liouvillian superoperator $\mathcal{L}$ defined by the right-hand-side of Eq. (27) has the following matrix representation:

[TABLE]

The steady state for which $\mathcal{L}(\rho_{0})=0$ is given by the unique null eigenvector of $L$:

[TABLE]

with

[TABLE]

Now assume that $\Delta$ is slowly modulated in the time interval $t\in[0,\tau]$. For example, let us take

[TABLE]

In the quasi-static limit $\tau\rightarrow\infty$ we expect the system to follow the instantaneous steady state (34), which now depends on time through the parameter $\Delta(t)$.

For a finite $\tau$ we can apply the perturbation theory presented in the main text. We want to compare the exact numerical solution with the corresponding [math]th, $1$st and $2$nd order approximations which are given respectively by:

[TABLE]

where the pseudo-inverse of the projected Liouvillian can be explicitly computed. Omitting the time parameter to simplify the notation, the latter is given by:

[TABLE]

We can finally compare the numerical solution of the master equation (27) (for a maximally mixed initial state) with the corresponding analytical approximations (37,38,39). The results are reported in Fig. 2 where we observe the, after an initial transient depending on the initial conditions, the dynamics is well captured by the perturbative approximations.

Appendix B Finite-time isothermal process of a two-level system

In this section we explicitly compute the finite time corrections to the heat absorbed by a two-level system during an isothermal process in which the energy splitting $\omega(t)$ is changed in time within the interval $t\in[0,\tau]$. The system can be described by the previous master equation (27) setting $V_{I}(t)=0$ (the density matrix remains diagonal). In interaction picture the Hamiltonian evolution is canceled and the master equation depends on the modulated frequency only via the damping rate and the mean excitation number:

[TABLE]

Setting $\Delta=0$ in Eq. (33) and dropping the time variable in the parameters $\gamma(t)$ and $N(t)$, the Liouvillian matrix simplifies to the following block form

[TABLE]

whose instantaneous eigenvector is the Gibbs state

[TABLE]

where

[TABLE]

is the $z$-coordinate of the Gibbs state in the Bloch sphere. The inverse of the projected Liouvillian (A) reduces in this case to:

[TABLE]

and it is easy to check that $[L(t)P]^{-1}$ is diagonal in the basis of vectorized Pauli matrices, in particular:

[TABLE]

Using the previous equations, the perturbative solution (11) presented in the main text becomes:

[TABLE]

where we have reintroduced density matrix representation. In terms of time rescaled variables $\tilde{X}(t^{\prime})=X(t^{\prime}\tau)$, the perturbation terms defined in Eq. (10) of the main text are given by:

[TABLE]

The heat absorbed by the system during the process can be expanded with respect to $1/\tau$ as $Q=\sum_{j=0}^{\infty}Q_{j}\tau^{-j}$, where

[TABLE]

Analogously, for the work done on the system we have:

[TABLE]

Appendix C Finite-time Carnot cycle of a two-level system

In this appendix, we consider a Carnot cycle performed on a two-level system with Hamiltonian $H=\sigma_{z}\hbar\omega(t)/2$. The protocol is defined as follow:

Isothermal expansion: the system is put in contact with a hot bath of temperature $T_{H}$ and the frequency is reduced according to the function $\omega^{H}(t)=\omega_{0}(\cos(\pi t/\tau_{H})+1)$, in a time interval $t\in[0,\tau_{H}]$. 2. 2.

Adiabatic expansion: the system is decoupled from the bath and the frequency is shifted to $\omega^{C}(0)=(T_{C}/T_{H})\omega^{H}(\tau_{H})$. Actually this step is useless in this case since $\omega^{H}(\tau_{H})=0$, but we keep it just for clarity. 3. 3.

Isothermal compression: the system is put in contact with a cold bath of temperature $T_{C}$ and the frequency is increased according to the function $\omega^{C}(t)=(T_{C}/T_{H})\omega_{0}(\cos[\pi(1-t/\tau_{C})]+1)$, in a time interval $t\in[0,\tau_{C}]$. 4. 4.

Adiabatic compression: the system is decoupled from the bath and the frequency is shifted to back to $\omega^{H}(0)=(T_{C}/T_{H})\omega^{C}(\tau_{C})$.

Notice that the driving protocol is such that the Bloch coordinate $z(t)=\tanh(\beta(t)\omega(t)/2)$ of the quasi-static Gibbs state $\rho_{0}(t)$ remains continuous and differentiable during the full cycle. In particular, at the extreme of the two isothermal processes, $\dot{\omega}(t)=0$ ensuring the continuity of $\dot{z}(t)$ also at such critical points. Moreover the shape of the hot and of the cold isothermal processes respect the time reversal symmetry, such that

[TABLE]

Thus we can apply the perturbation theory presented in the main text for estimating the heat exchanged in one cycle. From (60) we have

[TABLE]

where the function

[TABLE]

depends on the particular shape of the process only through $\tilde{z}(t)$. Thus, the time reversal symmetry (63) implies that the two integrals in (66) and (67) are equal and we get, consistently with Eq. (23) of the main text,

[TABLE]

Optimizing the output power with respect to $\tau_{H}$ and $\tau_{C}$ and using Eq. (69), one obtains

[TABLE]

which is the same formula for the efficiency at maximum power that we derived for general systems in the main text.

For the same optimal values of $\tau_{H}$ and $\tau_{C}$, it is interesting now to compare Eq. (70) with the corresponding exact result

[TABLE]

where the heat exchanged in the two processes can be computed from the exact numerical solution of the master equation (evolved for sufficiently many cycles such that any memory of the initial conditions is suppressed). The comparison between the two values of the efficiency at maximum power is reported in Fig. 3. Notice that exact numerical results are missing for $T_{C}/T_{H}$ below a certain threshold. This is due to the fact that for too small values of $T_{C}/T_{H}$, the exact dynamics deviates so much from the idealized Carnot cycle that the system does not produce positive work anymore and cannot be interpreted as an engine.

Appendix D Scaling properties of thermal Liouvillians and heat fluxes

In this final appendix we derive a universal scaling property of thermal Liouvillians which implies the general validity of Eq. (23) presented in the main text and already derived, for the particular case of a two-level system, in Eq. (69). Let us start from a typical microscopic model used in the derivation of Markovian master equations Breuer2002A ; Rivas2012A , in particular we closely follow the continuous bath notation used in Rivas2012A . Consider a system of Hamiltonian $H_{S}$ coupled to a heat bath of harmonic oscillators such that the full Hamiltonian can be written as

[TABLE]

where

[TABLE]

is the Hamiltonian of the bath, $a_{\omega}$ being the annihilation operator of a mode with frequency $\omega$ obeying the bosonic commutation rule $[a_{\omega},a_{\omega}^{\dagger}]=\delta(\omega-\omega^{\prime})$. For the interaction we consider the following potential

[TABLE]

where $A$ is an arbitrary system operator while the coupling function $h(\omega)$ simultaneously contains information about the coupling strength and the degeneracy of the bath modes. The square of $h(\omega)$ is known as the bath spectral density, which is usually modeled as a power law (with a suitable cut-off at large frequencies):

[TABLE]

The exponent $\alpha$ is a characteristic parameter of the system-bath interaction: e.g. for $\alpha=0$ the bath is usually called flat, while for $\alpha=1$ the bath is said to be Ohmic. The operator $A$ can always be decomposed in terms of the eigenoperators of $H_{S}$:

[TABLE]

where

[TABLE]

Moving in interaction picture with respect to $H_{S}+H_{B}$, within the usual weak-coupling and secular approximation, one can derive the following Markovian master equation for the reduced state of the system Breuer2002A ; Rivas2012A ,

[TABLE]

where $D_{X}(\rho)$ is the Gorini-Kossakowski-Sudarshan-Lindblad dissipator defined as

[TABLE]

and

[TABLE]

One can also show Breuer2002A ; Rivas2012A that the steady state $\mathcal{L}[\rho_{0}]=0$ of the master equation is the Gibbs state $\rho_{0}\propto\exp{(-\beta H_{S})}$ associated to the system Hamiltonian $H_{S}$.

Consider now another system-bath configuration in which we perform the rescaling: $H_{S}\rightarrow\lambda H_{S}$ and $\beta\rightarrow\beta/\lambda$ , for some $\lambda>0$. Clearly the steady state remains the same but what can be said about the master equation? The quantum Liouvillian defined by the right-hand-side of (78) is completely determined by the system Hamiltonian $H_{S}$, the inverse temperature $\beta$ of the bath, the damping constant $\gamma_{0}$ and the spectral density exponent $\alpha$. Now, by simply following all the steps of the previous microscopic derivation and observing that the rescaling of $H_{S}$ implies that all energy gaps change as $\omega\rightarrow\lambda\omega$, it is easy to check that:

[TABLE]

Basically, Eq. (81) states that a rescaling of $H_{S}$ and of $\beta$ with inverse factors leaves the Liouvillian invariant up to a renormalization of the damping rate. In interaction picture with respect to $H_{S}$, one may also rewrite Eq. (81) simply as

[TABLE]

The validity of Eq. (82) also extends to more general interactions of the form $V=\sum_{l}A_{l}\otimes B_{l}$, conditioned on the requirement that all the spectral densities associated to the different bath operators $B_{l}$ should scale with the same frequency exponent, i.e. $J_{l}(\omega)=\gamma_{l}\omega^{\alpha}/2\pi$.

Now, consider the first order correction to the heat exchanged during an arbitrary isothermal process which, according to the perturbation theory presented in the main text is given by

[TABLE]

How does the quantity (83) changes if we preform the rescaling $H\rightarrow\lambda H$ and $\beta\rightarrow\beta/\lambda$ that we have previously discussed? We know that the Gibbs state $\rho_{0}$ is obviously invariant, moreover exploiting Eq.(82), we get:

[TABLE]

We can finally consider the Carnot cycle discussed in the main text. By definition, up to an irrelevant time-inversion, the cold isothermal compression is related to the hot isothermal expansion by exactly the same kind of rescaling that we have just considered above with $\lambda=T_{C}/T_{H}$. Therefore from (84) we get

[TABLE]

independently on the particular shape of the driving protocol, corresponding to Eq. (23) of the main text.