Non-Markov Enhancement of Maximum Power for Quantum Thermal Machines

TL;DR

This paper investigates how non-Markovian dynamics can enhance the maximum power output of quantum thermal engines, showing that reservoir correlations can be exploited for improved control and performance.

Contribution

It introduces a model incorporating reservoir correlations via Hamiltonian coupling and demonstrates how non-Markovian effects can accelerate control and boost power output.

Findings

Non-Markovianity can accelerate control processes.

Reservoir correlations improve engine power.

Enhanced performance in Carnot and Otto cycles.

Abstract

In this work we study how the non-Markovian character of the dynamics can affect the thermodynamic performance of a quantum thermal engine, by analysing the maximum power output of Carnot and Otto cycles departing from the quasi-static and infinite-time-thermalization regime respectively, introducing techniques for their control optimization in general dynamical models. In our model, non-Markovianity is introduced by allowing some degrees of freedom of the reservoirs to be taken into account explicitly and share correlations with the engine by Hamiltonian coupling. It is found that the non-Markovian effects can fasten the control and improve the power output.

Equations

$$\dot{\hat{\rho }}(t)={\mathcal{L}}_t[\hat{\rho }(t)]:=-i[{\hat{H}}_t,\hat{\rho }(t){]}_{-}+{\mathcal{D}}_t^{(j)}[\hat{\rho }(t)],$$

$${\mathcal{D}}_t^{(j)}={\mathcal{D}}^{(j)}({\hat{H}}_t).$$

$${\hat{\Omega }}_{{\hat{H}}_t}^{(j)}:=\frac{e^{-{\beta }_j{\hat{H}}_t}}{\text{Tr}[e^{-{\beta }_j{\hat{H}}_t}]},$$

$${\mathcal{D}}_t^{(j)}[\hat{\rho }]=0\iff \hat{\rho }={\hat{\Omega }}_t^{(j)}.$$

$${\mathcal{L}}_t[\hat{\rho }]=0\iff \hat{\rho }={\hat{\Omega }}_{{\hat{H}}_t}^{(j)}.$$

$$\underset{t\to \infty }{\lim }\hat{\rho }(t)={\hat{\Omega }}_{\hat{H}}^{(j)},$$

$$E(t):=\text{Tr}[\hat{\rho }(t){\hat{H}}_t].$$

$$dW(t):=\text{Tr}[\hat{\rho }(t)d{\hat{H}}_t],$$

$$d{Q}_j(t)$$

$$=$$

$$F(\hat{\rho }(t),{\hat{H}}_t):=E(t)-S(\hat{\rho }(t))/{\beta }_j,$$

$$F(\hat{\rho }(t),{\hat{H}}_t)-F({\hat{\Omega }}_{{\hat{H}}_t}^{(j)},{\hat{H}}_t)=S(\hat{\rho }(t)\parallel {\hat{\Omega }}_{{\hat{H}}_t}^{(j)})/{\beta }_j,$$

$$dF(\hat{\rho }(t),{\hat{H}}_t)\le dW(t)⟺{\beta }_{j}d{Q}_j(t)\le dS(\hat{\rho }(t)),$$

$$\hat{\rho }(t)\simeq {\hat{\Omega }}_{{\hat{H}}_t}^{(j)}.$$

$$\hat{\rho }(t)={\hat{\rho }}^{(0)}(t)+{\hat{\rho }}^{(1)}(t)+...,$$

$${\hat{\rho }}^{(1)}(t)=({\mathcal{L}}_t\mathcal{P}{)}^{-1}[{\dot{\hat{\rho }}}^{(0)}(t)],$$

$$d{Q}_j(t)\simeq d{Q}_j^{(0)}(t)+d{Q}_j^{(1)}(t),$$

$$\displaystyle dQ_{j}^{(0)}(t):=\text{Tr}\Big{[}\hat{H}_{t}d\hat{\rho}^{(0)}(t)\Big{]}=dS^{(0)}(t)/\beta_{j}\;,$$

$$\displaystyle dQ_{j}^{(1)}(t)=\text{Tr}\Big{[}\hat{H}_{t}d\hat{\rho}^{(1)}(t)\Big{]}\;,$$

$${\hat{H}}_t=ϵ(t){\hat{\sigma }}^z/2,$$

$${\hat{\Omega }}_{{\hat{H}}_t}^{(j)}={\hat{\Omega }}_{ϵ(t)}^{(j)}:=p_j(ϵ(t))|0⟩⟨0|+(1-p_j(ϵ(t)))|1⟩⟨1|,$$

$$p_j(ϵ):=\frac{1}{1+e^{-{\beta }_jϵ}},$$

$${ϵ}_3:={ϵ}_2{\beta }_{\mathrm{C}}/{\beta }_{\mathrm{H}},$$

$${ϵ}_4:={ϵ}_1{\beta }_{\mathrm{C}}/{\beta }_{\mathrm{H}},$$

$$S({\hat{\Omega }}_{ϵ}^{(j)})=-p_j(ϵ)\ln p_j(ϵ)-(1-p_j(ϵ))\ln (1-p_j(ϵ)),$$

$$d{Q}_{\mathrm{C}}(t)\simeq d{Q}_{\mathrm{C}}^{(0)}(t)$$

$$d{Q}_{\mathrm{H}}(t)\simeq d{Q}_{\mathrm{H}}^{(0)}(t)$$

$$Q_{\mathrm{ABS}}$$

$$Q_{\mathrm{REL}}$$

$$p_{\mathrm{H}}({ϵ}_3)=p_{\mathrm{C}}({ϵ}_2)p_{\mathrm{H}}({ϵ}_4)=p_{\mathrm{C}}({ϵ}_1),$$

$$\Delta Q_{\mathrm{H}}^{(0)}=-\frac{{\beta }_C}{{\beta }_H}\Delta Q_{\mathrm{C}}^{(0)},$$

Full text

I Introduction

Quantum Thermodynamics Gemmer et al. (2009); Goold et al. (2016); Vinjanampathy and Anders (2016) was born and rapidly grew in the last decades. Fuelled by high experimental control of quantum systems and engineering at microscopic scales, one of the central goals of physicists is to push the limits of conventional thermodynamics, and the extension of standard models and cycles to include quantum effects and small ensemble sizes. Beyond the drive to clarify fundamental physical issues, these models may also turn out to be relevant from a more practical point of view: it is expected that industrial need for miniaturisation of technologies will benefit from the understanding of quantum thermodynamic processes. In both biology, for example, and nanotechnology, where the benefits from a cooling at the atomic scales are clear, refrigerators models Linden et al. (2010); Skrzypczyk et al. (2011) based on quantum thermal machines could find actual application. Moreover, proposals for experimental realisations of quantum engines were made considering various physical platforms, and many were actually realised Koski et al. (2014); Pekola (2015); Batalhão et al. (2014); An et al. (2015); Roßnagel et al. (2016); Zhang et al. (2014a); Abah et al. (2012); Ronzani et al. (2018); Chen (1994); Rezek and Kosloff (2006); Watanabe et al. (2017); Scully et al. (2011); Correa et al. (2013); Dorfman et al. (2013); Brunner et al. (2014); Zhang et al. (2014a); Campisi and Fazio (2016); Brandner et al. (2017).

Thermodynamics is, par excellence, a theory involving non-isolated systems, and it must take into account the interaction and evolution induced by external degrees of freedom on a working medium. The description of open quantum systems Breuer and Petruccione (2002) needs however, especially in cases where the number of degrees of freedom of the surroundings is big, an effective description on the local degrees of freedom by means of some approximation or assumption. The most important class of simplified dynamics of open systems goes under the name of Markovian dynamics. From the physical point of view, Markovianity is associated to systems interacting with large, unperturbed environments that “spread away the information” contained in the system, while on the formal side different definitions of quantum Markovianity Rivas et al. (2014); Breuer et al. (2016) were introduced in the literature. We stand by the approach (although the model we will consider is non-Markovian even for stronger definitions of quantum Markovianity Rivas et al. (2014); Breuer et al. (2016)) which identifies the Markovian character of a quantum process with its CP-divisibility Rivas et al. (2010) hence admitting a first order Master Equation (ME) that can be casted in the Gorini-Kossakowski-Sudarshan-Lindblad form (GKSL) Gorini et al. (1976); Lindblad (1976).

Recent works have started to investigate how the breaking of the Markovianity in quantum dynamics can affect control and performance of quantum thermodynamic systems, motivated both by the necessity to overcome the approximation on very small systems, and by the speculation of non-Markovianity possibly being an actual resource in practical tasks, see e.g. Refs. Bhattacharya et al. (2018); Bylicka et al. (2016); Mirkin et al. (2017); Mukherjee et al. (2015); Raja et al. (2018); Reich et al. (2015); Thomas et al. (2018); Zhang et al. (2014b); Basilewitsch et al. (2017); Pezzutto et al. (2018). We contribute here considering two archetypical classes of thermal engines, i.e. the quantum Carnot cycle and the quantum Otto cycle Quan et al. (2007); Karimi and Pekola (2016); Kosloff and Rezek (2017); Watanabe et al. (2017); Rezek and Kosloff (2006); Abah et al. (2012) which use as working medium a two-level (qubit) system coupled to two thermal reservoirs while being externally driven. For these models we simulate non-Markovian effects by splitting the degrees of freedom of the system environmental baths into a local contribution, which we treat dynamically, and a remote component which instead is described in terms of an effective GKSL Master Equation that tends to drive the rest of the model into thermal equilibrium. In this configuration it can be shown that the coupling with the local bath components ignites the non-Markovian behaviour of the model whose effects can then be tested in terms of the engine performance. In particular, performing an optimization on the external driving, we show that, both in the Carnot and Otto scheme, the maximum power extractable improves with respect to the Markovian limit. To do this we first discuss both cycles in the finite-time regime; to solve the dynamics and optimize the control for the Carnot cycle, we use the powerful technique introduced in Cavina et al. (2017a) (Slow-Driving approximation, or S-D), which efficiently solves the approximate dynamics of a system slowly perturbed from thermalization. For the Otto case we use exact solutions.

The article is structured as follows:

Sec. II and Sec. III are devoted to introduce the technical tools we use to derive the results of Sec. IV. Specifically in Sec. II we discuss the physics of an externally controlled, quantum thermal machine introducing the notation in Sec. II.1, drawing general thermodynamic considerations in Sec. II.2 and reviewing some basic facts about the S-D approximation method Cavina et al. (2017a) in Sec. II.3. In Sec. III instead we analyze the performances of some thermodynamic cycles. In particular Sec. III.1 is devoted to study the quantum Carnot cycle in the quasi-static approximation and its first order S-D corrections, recovering some known results in a slightly broader context. Sec. III.2 instead focuses on the Otto cycle. In Sec. IV we finally introduce the specific non-Markovian model. Using the preceding section results, we then show how the power output of the cycles gets affected both for the Carnot machine (Sec. IV.2) and for the Otto machine (Sec. IV.3). In Sec. IV.4 an argument is presented to interpret the results obtained, focusing on why the information flow induced by non-Markovianity can fasten the speed of thermalization. Comments and conclusions are presented in Sec. V while the Appendix contains some technical derivations.

II Quantum thermal machines in the Markovian regime

In this section we review some basic facts about quantum thermal machines in the Markovian regime, setting the notation and developing the tools that we shall later employ for analysing the non-Markovian case.

II.1 The setup

Consider a quantum working medium ${\cal S}$ characterized by a time-dependent internal Hamiltonian $\hat{H}_{t}$ which can be externally controlled via some classical pulses. As schematically shown in Fig. 1 ${\cal S}$ is coupled to a collection of external thermal baths $\{{\cal E}_{j}\}$ characterized by temperatures $T_{j}$, which are also externally controlled to allow selective activation and deactivation. In particular we shall assume at each time $t$ only one of the baths is actively coupled with the working medium. Accordingly, enforcing the Markovian character in the system-bath interactions, we describe the evolution of ${\cal S}$ in terms of a Master Equation Breuer and Petruccione (2002) associated with a step-continuous generator ${\cal L}_{t}$ which, on the time interval ${\cal I}_{j}$ where only the $j$-th bath interaction is active, writes

[TABLE]

where $\hat{\rho}(t)$ is the density matrix of ${\cal S}$ at time $t$, $[\cdots,\cdots]_{-}$ is the commutator symbol, and where finally $\mathcal{D}_{t}^{(j)}$ is the GKSL dissipator Gorini et al. (1976); Lindblad (1976) mimicking the interaction with ${\cal E}_{j}$ (hereafter for easy of notation we set both the Plank and the Boltzmann constant equal to one, i.e. $\hbar=k_{\rm B}=1$).

As indicated by the notation the $\mathcal{D}_{t}^{(j)}$s exhibit an explicit time dependence which, in a weak-coupling regime, we assume to be a direct consequence of the modulations affecting the system Hamiltonian, i.e.

[TABLE]

Furthermore, in order to impose proper thermalization conditions on the scheme we require $\mathcal{D}_{t}^{(j)}$ to admit the instantaneous Gibbs state

[TABLE]

with $\beta_{j}:=1/T_{j}$ being the associated inverse temperature, as unique fixed point, i.e.

[TABLE]

Notice that the functional dependence of $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ with respect to $\hat{H}_{t}$, ensures that the requirement Eq. (4) is fully compatible with (2) and it implies that for $t\in{\cal I}_{j}$, $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ is also the unique fixed point of the full generator ${\cal L}_{t}$, i.e.

[TABLE]

Explicit examples of dissipators $\mathcal{D}_{t}^{(j)}$ obeying the above constraints are presented in Appendix A, here we only remark that they have been extensively used in the characterization of equilibration processes induced by fermionic or bosonic baths, see e.g. Refs. Cavina et al. (2017b); Gardiner et al. (2004); Breuer and Petruccione (2002); Esposito et al. (2010). In the absence of Hamiltonian modulations (i.e. for $\hat{H}_{t}=\hat{H}$ constant), Eqs. (4) and (5) ensure that if ${\cal S}$ is left in contact with the $j$-th bath, it will be forced by (1) to asymptotically reach thermal equilibrium at temperature $T_{j}$, i.e.

[TABLE]

irrespectively from the initial condition of the problem.

II.2 Energy exchanges and thermodynamic consistency

Within the above theoretical framework the internal energy $E(t)$ of ${\cal S}$ can be identified with the expectation value of $\hat{H}_{t}$ on $\hat{\rho}(t)$, i.e.

[TABLE]

Its infinitesimal variation comprises two terms which, following the canonical approach of Refs. Alicki (1979); Anders and Giovannetti (2013); Kieu (2004); Vinjanampathy and Anders (2016), are associated respectively with a work (performed on ${\cal S}$) contribution

[TABLE]

and with a heat (absorbed by ${\cal S}$) contribution

[TABLE]

where in the second identity we make explicit use of Eq. (1), ${\cal E}_{j}$ being the only bath that is coupled with ${\cal S}$ at time $t$. It is worth stressing that the consistency of the above identifications is explicitly justified by the Markovian character of the thermalizing process we are considering. To see this let us introduce the functional Parrondo et al. (2015); Esposito and den Broeck (2011)

[TABLE]

where $S(\hat{\rho}(t)):=-\mbox{Tr}[\hat{\rho}(t)\ln\hat{\rho}(t)]$ is the von Neumann entropy of $\hat{\rho}(t)$. Exploiting the formal connection between informational and thermodynamical entropy, the quantity (11) can be identified with the counterpart of the free energy functional of classical equilibrium thermodynamics. One can easily verify that it obeys the identity

[TABLE]

where $S(\hat{\rho}_{1}\parallel\hat{\rho}_{2}):=-S(\hat{\rho}_{1})-\mbox{Tr}[\hat{\rho}_{1}\ln\hat{\rho}_{2}]$ is the relative entropy functional Holevo (2012). The latter is know to be decreasing when the same completely positive mapping acts on both its argument: accordingly, given that the dynamical generator ${\cal L}_{t}$ of Eq. (1) is guaranteed to grant complete positive evolution and using the invariance (5) of $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ we can claim that $S({\cal L}_{t}[\hat{\rho}(t)]\parallel\hat{\Omega}^{(j)}_{\hat{H}_{t}})\leq 0$. Inserting this into (12) we can establish that the time derivative of the l.h.s. must be upper bounded by the quantity $S(\hat{\rho}(t)\parallel\frac{d}{dt}\hat{\Omega}^{(j)}_{\hat{H}_{t}})/\beta_{j}$, which after proper reordering of the various terms leads to the inequality

[TABLE]

that is an instance of the 2nd Law of thermodynamics providing an operational justification for the definitions (8) and (9).

II.3 Thermodynamic cycles

Integrating Eq. (1) we can now analyze the work production rates, their associated efficiencies, and the corresponding heat fluxes, of thermodynamic cycles where the system ${\cal S}$ is externally driven by an assigned modulation of the Hamiltonian $\hat{H}_{t}$ while being put in selective contact with the baths ${\cal E}_{j}$s – see below. Unfortunately the presence of Hamiltonian modulations makes typically Eq. (1) hard to solve. Yet assuming the time scale at which (6) takes places to be short enough, one expects ${\cal S}$ to have enough time to adiabatically follow the instantaneous fixed points of Eq. (3), obtaining

[TABLE]

This is the standard quasi-static regime where the working medium is always at thermal equilibrium with one of the baths. Departing from this scenario one enters the regime of Finite Time Thermodynamics (FTT) Andresen et al. (1984), where the time-scales on which the external controls responsible for the modulations of $\hat{H}_{t}$ occur, begin to compete with the thermalization times. In what follows we shall study this complex regime by adopting the Slow-Driving (S-D) approximation technique introduced in Ref. Cavina et al. (2017a). The latter is a perturbative approach which can be applied to study deviations from Eq. (14) in the limit of slow variation of $\mathcal{L}_{t}$. As we detail in Appendix B, the S-D approximation can be used as a way for putting on firm ground some of the assumptions typically adopted in FFT analysis. It accounts in expressing the solution of Eq. (1) as an expansion series with a perturbation parameter given by the ratio $\tau_{R}/\tau$ between the typical timescale $\tau\sim\|\dot{\mathcal{L}}_{t}/\mathcal{L}_{t}\|$ associated with the variation of the dynamics generator, and the typical relaxation time $\tau_{R}$ governing the convergence of the limit (6). At the lowest orders one has

[TABLE]

with $\hat{\rho}^{(0)}(t):=\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ being the zero-th order term, while the first order correction $\hat{\rho}^{(1)}(t)$ is obtained as Cavina et al. (2017a)

[TABLE]

where $\mathcal{P}$ is the projector on the null-trace subspace of linear operators (its presence being required to make $\mathcal{L}_{t}$ invertible, under the assumption of unique null eigenstate). Therefore, by direct substitution in Eq. (10) we get

[TABLE]

where

[TABLE]

is the quasi-static contribution which, by using the fact that $\hat{\rho}^{(0)}(t)$ is the Gibbs state $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$, we expressed in terms of the infinitesimal increment the von Neumann entropy $S^{(0)}(t):=-\mbox{Tr}[\hat{\rho}^{(0)}(t)\ln\hat{\rho}^{(0)}(t)]$ of the latter, and where

[TABLE]

is the first order correction term.

III Thermodynamic cycles optimisation

In this section we will show how it is possible to optimize the control on a quantum engine in order to maximize is performance, i.e. its power output, addressing the paradigmatic case of Quantum Carnot and Otto cycles performed on a two-level (qubit) system ${\cal S}$ which evolves under the influence of a hot bath H and a cold bath C, the modulation of its Hamiltonian being associated with control pulses that act on its energy gap $\epsilon(t)\geq 0$, i.e.

[TABLE]

with $\hat{\mathcal{\sigma}}^{z}$ being the third Pauli matrix, with eigenstates $|0\rangle$ and $|1\rangle$. It’s not difficult to generalise these cycles (in the quasi-static regime) to more general Hamiltonians. While in deriving the above considerations we shall make explicit reference to the expressions we developed in Sec. II for the Markovian regime, we stress that the results we obtain also hold for non-Markovian dynamics, as we shall use them later in Sec. IV to analyse the non-Markovian model we present.

III.1 Quantum Carnot Cycle

A Quantum Carnot cycle is identified with a 4 steps process inspired directly by its classical counterpart, that is two isothermal strokes where the Hamiltonian of ${\cal S}$ is modulated while keeping the system in thermal contact with one of the two baths, alternated with two iso-entropic (adiabatic) strokes, where instead the Hamiltonian undergoes to instantaneous sudden switches (quenches). In the ideal quasi-static limit (14) the operations are performed slowly enough to allow the system to be in thermal equilibrium at every instant, i.e. states which for the Hamiltonian (20) can be expressed as

[TABLE]

with

[TABLE]

being the associated ground state population.

In this case the 4 steps of the cycle are as in Figure 2:

while being coupled to the cold reservoir C, the energy gap is modified continuously and monotonically, from the initial value $\epsilon_{1}$ to $\epsilon_{2}\geq\epsilon_{1}$ (more precisely we require $\epsilon(t)$ to be continuous and differentiable with first order derivative which is not negative);

• 2)

with the system isolated from the reservoirs, a quench is now performed by suddenly taking the gap from $\epsilon_{2}$ to

[TABLE]

which by construction is larger than $\epsilon_{2}$, i.e. $\epsilon_{3}\geq\epsilon_{2}$;

• 3)

while being coupled to the hot reservoir H, the energy gap is then modified continuously, and monotonically, from $\epsilon_{3}$ to

[TABLE]

that automatically fulfils the constraint $\epsilon_{1}\leq\epsilon_{4}\leq\epsilon_{2}{\beta_{\rm C}}/{\beta_{\rm H}}=\epsilon_{3}$ (again, more precisely we require $\epsilon(t)$ to be continuous and differentiable with first order derivative that is non-positive);

• 4)

finally isolating the system a quench is performed to restore the gap at the initial value $\epsilon_{1}$.

It is worth pointing out that the continuity requirement of $\epsilon(t)$ during the steps 1) and 3) is inserted in order to make sure that one could later on apply the S-D expansion which needs to have a zero-th order contribution of term differentiable – see Eq. (15). More specifically in what follows we shall require $\epsilon(t)$ to have null first order derivative at the extrema of the isotherms. This is a technical assumption which we introduce in order to ensure the solution of the dynamics (15) to be continuous and differentiable also in proximity of the quenches (where it coincides with the Gibbs state (21)), which in turn implies that no first order correction (16) at the extremal points of the isothermal strokes has to be expected. The monotonicity behaviour of $\epsilon(t)$ during the steps 1) and 3) is instead motivated by energetic considerations. As a matter of fact having set the gap to evolve monotonically from $\epsilon_{1}\geq 0$ to $\epsilon_{2}\geq\epsilon_{1}$, we can ensure that at each instant of step 1) the system always releases heat to the cold bath without absorbing it: this can be easily verified by observing that the von Neumann entropy of a Gibbs state (21) writes

[TABLE]

which is monotonically decreasing with $\epsilon$, and from the fact that at the lowest order in the expansion (18) the associated incremental heat can be expressed as

[TABLE]

Similarly having ensured that in step 3) the value of the gap decreases monotonically from $\epsilon_{3}$ to $\epsilon_{4}\geq\epsilon_{3}$, we can guarantee that the heat in the process is always absorbed from the bath H, i.e.

[TABLE]

Thanks to these properties, and by the observation that of course no heat is exchanged between ${\cal S}$ and the baths during the steps 2) and 4), the total heat absorbed by the working medium in a cycle can be obtained by integrating (27) over the full duration of step 3), i.e.

[TABLE]

while the total released heat is given by

[TABLE]

Notice also that the constraints (23) and (24) impose

[TABLE]

which implies $\hat{\Omega}^{({\rm H})}_{\epsilon_{3}}=\hat{\Omega}^{({\rm C})}_{\epsilon_{2}}$ and $\hat{\Omega}^{({\rm H})}_{\epsilon_{4}}=\hat{\Omega}^{({\rm C})}_{\epsilon_{1}}$. Accordingly by direct inspection of (III.1) and (III.1) we obtain the fundamental identity

[TABLE]

which we expressed in terms of the simplified notation $\Delta Q^{(0)}_{j}:=\int_{{\cal I}_{j}}dQ^{(0)}_{j}(t)$. Now, since the work produced by ${\cal S}$ on a cycle can be identified with $Q_{\rm ABS}+Q_{\rm REL}$ by invoking the internal energy conservation, the efficiency (work done over heat absorbed) of the process can be shown to correspond to the Carnot efficiency $\eta_{\text{c}}:=1-\frac{\beta_{\rm H}}{\beta_{\rm C}}$. Indeed

[TABLE]

the last identity following directly from (31). It is worth stressing that Eqs. (26) and (32) are universal results that do not depend on the specific structure of the generators $\mathcal{D}_{t}^{(j)}$ entering the system ME. This is a consequence of the quasi-static approximation (14) in which, as in classical thermodynamics, complete thermalization is allowed at any time in contact with a thermal source: in this regime no explicit dynamics as in Eq. (1) is needed to describe the thermodynamics of the engine, neither the exact temporal dependence of the control $\epsilon(t)$, except the properties of the equilibrium state (3) and the knowledge of the Hamiltonian at the turning points of the protocol. All this of course holds true as long as we can neglect the first-order contributions in the S-D expansion (15). To account for them we now use (17) to refine Eqs (III.1) and (III.1), writing $Q_{\rm REL}\simeq\Delta Q^{(0)}_{\rm C}+\Delta Q^{(1)}_{\rm C}$ and $Q_{\rm ABS}\simeq\Delta Q^{(0)}_{\rm H}+\Delta Q^{(1)}_{\rm H}$ with

[TABLE]

obtaining

[TABLE]

where in the last identity we employed (17) to express the ratio $\Delta Q^{(0)}_{\rm C}/\Delta Q^{(0)}_{\rm H}$ in terms of the Carnot efficiency and for $j\in\{{\rm H,C}\}$ introduced the parameter

[TABLE]

to gauge the ratio between the first and the zero-th order heat contributions associated with the $j$-th bath. In a similar fashion we can also express the power $P$ associated with the work production per cycle. Indicating hence with $\tau_{\rm H}$ and $\tau_{\rm C}$ the durations of the transformations 1) and 3) (the only being time-consuming given that step 2) and 4) are assumed to be instantaneous), we write

[TABLE]

where we used Eqs. (31) and (35).

III.1.1 Performance optimization in the S-D regime

To proceed with our analysis we need to provide some details on the system ME and in particular on the GKSL dissipators which define it. As a preliminary step, however we observe that thanks to our choice (20) we can express Eq. (10) as

[TABLE]

where for $k,k^{\prime}=0,1$, ${\rho}_{kk^{\prime}}(t):=\langle k|\hat{\rho}(t)|k^{\prime}\rangle$ are the matrix elements of $\hat{\rho}(t)$ with respect to the eigenbasis of $\hat{H}_{t}$ and where in the second identity we use the normalization condition $\mbox{Tr}[\hat{\rho}(t)]=1$ to write $d{\rho}_{11}(t)=-d{\rho}_{00}(t)$. Due to linearity Eq. (37) applies to all orders of the S-D expansion (17), implying in particular that Eqs. (18), (19) take the form

[TABLE]

where ${\rho}^{(0)}_{00}(t)$ and ${\rho}^{(1)}_{00}(t)$ are respectively the zero-th and first order contribution to the population of the ground state of ${\cal S}$. The first of these two terms is nothing but the function (22), i.e. ${\rho}^{(0)}_{00}(t)=p_{j}(\epsilon(t))$. The second instead can be determined exploiting Eq. (16). In particular due to the linearity of operators in Eq. (16) and the one-parameter dependence of $\hat{\rho}^{(0)}$ it is possible to draw, in full generality, the following formal connection between ${\rho}^{(1)}_{00}(t)$ and the function $p_{j}(\epsilon(t))$ which, effectively, becomes the real control parameter of the setting. Specifically we get

[TABLE]

where $A_{j}$, which we dub the S-D amplitude of the problem, quantifies how large is the first order correction determining the relaxation timescale of the setup. In general, besides depending on the the parameters of the model, the S-D amplitude is an explicit functional of ${p}_{j}(\epsilon(t))$, e.g. as in the case of dissipators $\mathcal{D}_{t}^{(j)}$ associated with Bosonic baths defined by Eq. (95) with rates as in (A) for which we get $A_{j}=(2p_{j}(\epsilon(t))-1)/\Gamma_{j}$. When considering instead as dissipators $\mathcal{D}_{t}^{(j)}$ the super-operators defined in Eq. (94) or those associated with fermionic baths defined by Eq. (95) with rates as in (A), one gets an S-D amplitude which is constant, i.e.

[TABLE]

with $\Gamma_{j}$ being a fundamental constant of the model. In what follows, for the sake of simplicity we shall focus on this special case: our finding however can be approximatively applied to all those configurations where, for all $t\in{\cal I}_{j}$, $A_{j}$ is a slowly varying functional of $p_{j}(\epsilon(t))$.

With the help of the above identities we can hence cast (33) as

[TABLE]

where in the first identity we used Eq. (22) to write $\epsilon(t)$ in terms of $p_{j}(t):=p_{j}(\epsilon(t))$, i.e. \epsilon(t)=\frac{1}{\beta_{j}}\ln\big{(}\tfrac{p_{j}(t)}{1-p_{j}(t)}\big{)}, and in the second we adopted integration by parts exploiting the fact that at the extrema of the isotherms steps the control functions have been set to have null first order derivative. Equation (41) should be compared with the zero-th order term $\Delta Q_{j}^{(0)}$ which we have already computed in the previous section and which, expressed in terms $p_{j}(t)$, results to be the integral of an exact differential that depends only on the initial and final values $p_{j}^{(in)}$ and $p_{j}^{(fin)}$ assumed on the interval ${\cal I}_{j}$, i.e.

[TABLE]

the last identity being an alternative way of expressing the entropy increment of the Gibbs state (21).

Our next problem is to determine which choices of $\epsilon(t)$, or equivalently of $p_{j}(t)$, can be used in order to guarantee better performances with respect to the quasi-static regime. To begin with it is worth stressing that from Eq. (41) it follows that for all choices of the control functions the first order correction term to the heat is always negative semi-definite, i.e.

[TABLE]

which in turn implies

[TABLE]

due to the positivity of $\Delta Q_{\rm H}^{(0)}$ and the negativity of $\Delta Q_{\rm C}^{(0)}$ (incidentally we observe that  (43) continues to hold by the same argument even if $A_{j}(p_{j})$ is not constant but explicitly dependent on the control $p_{j}(t)$). The first consequence of Eq. (44) is the fact that the efficiency $\eta$ of Eq. (34) cannot be larger than $\eta_{\text{c}}$, as one expects from the second principle of thermodynamics (formally speaking to show that $\eta\leq\eta_{\text{c}}$ we also need $|\alpha_{j}|\ll 1$ which however is always implicit assumed by the perturbative character of the S-D approach). At the level of the power (36) we notice instead that first order corrections explicitly depend on features which one may try to optimize with proper choices of the controls. For this purpose looking at the expression (41) we can isolate different contributions:

• •

Control speed: keeping the same shape (and extrema) for the driving protocol, we can modify its duration via the mapping $\tau_{j}\rightarrow\lambda\tau_{j}$ with $\lambda>0$. By a simple change of variable $t\rightarrow t/\lambda$ in Eq. (41), it is immediate to find that this induces the following rescaling

[TABLE]

while, of course, the zero-order terms $\Delta Q_{j}^{(0)}$ are unaffected;

• •

Control shape: over a fixed time length, we can clearly optimize with respect to the shape of the function $\epsilon(t)$, i.e. with respect to the function $p_{j}(t)$ under the constraint i), ii) and iii). Once more this will induce a modification of $\Delta Q_{j}^{(1)}$ while leaving unaffected the zero-order contribution terms;

• •

S-D amplitude selection: this is the main figure of merit after control optimisation. It merely consists in selecting different kind of bath-system interactions in order to influence the value of $\Delta Q_{j}^{(1)}$ via its dependence upon the S-D amplitude $A_{j}$ (this optimization will be specifically analyzed in the study of non-Markovian models).

Let us first analyze how the power $P$ is affected by Speed Control optimization. Using the scaling relations (45) we find that Eq. (36) changes as

[TABLE]

while the associated efficiency

[TABLE]

where we used (44) to rewrite $-\alpha_{\rm H}=|\alpha_{\rm H}|$, and where $\lambda_{\rm C},\lambda_{\rm H}>0$ represent the stretching of the intervals ${\cal I}_{\rm C}$ and ${\cal I}_{\rm H}$, respectively. A simple analytical study reveals that the function (46) admits a maximum for

[TABLE]

Replacing these values into (46) and (47) it is possible then to express the maximum power $P_{\max}$ and the correspondent efficiency at maximum power (EMP) that we indicate with the symbol $\eta^{*}$. For the sake of simplicity let us now suppose $\alpha_{\rm C}=|\alpha_{\rm H}|=\alpha$ and $\tau_{\rm C}=\tau_{\rm H}=\tau$, a regime attained for instance under symmetric bath couplings and driving assumptions, i.e. posing $A_{\rm C}=A_{\rm H}$ and requiring $\epsilon(t)$ during the hot isotherm to be the time-reversal of the cold isotherm – see however Ref. Cavina et al. (2017a) for an explicit treatment of the cases where this hypothesis is relaxed. In this case, using Eq. (31) and the quasi-static relation $\Delta S_{j}^{(0)}=\Delta Q_{j}^{(0)}/T_{j}$ by direct integration of (18), we find

[TABLE]

where we introduced the adimensional functional

[TABLE]

the variable $x:=t/\tau$ being a rescaled temporal coordinate and $q_{x}:=p(x\tau)$. Regarding the EMP instead we get

[TABLE]

which is the Curzon-Ahlborn efficiency Curzon and Ahlborn (1975), see also Appendix B.

Equation (50) implies that the maximum power $P_{\max}$ is inversely proportional to S-D amplitude, hence the larger values of $A_{j}$ is, the worse the effects on thermodynamic performance are (note that also the efficiency (34) worsen for larger $A_{j}$s). Regarding the shape-pulse optimization instead, remembering that the zero-th order terms $\Delta Q_{\rm H}^{(0)}$, and $\Delta Q_{\rm C}^{(0)}$ are not affected such choice, we observe that larger values of $P_{\max}$ are attained by minimizing the term ${\cal F}[q_{x}]$ appearing at the denominator for all possible choices of a monotonic, continuous, differentiable function $q_{x}\in[0,1]$, i.e.

[TABLE]

For each assigned initial and final values $q_{0}$ and $q_{1}$ of $q_{x}$, the problem can be solved by a variational study of the integrand, leading to solutions of the form $q_{x}=\frac{1+\cos({\Omega}(x+\varphi))}{2}$, cf. Appendix D. The resulting value of $P_{\max}$ obtained with such a driving reaches the maximal performances for $q_{1}=q_{0}+\varepsilon$, for which it is possible to obtain an analytic expression valid in the $\lim\varepsilon\rightarrow 0$, that is

[TABLE]

which we expressed in terms of the bath temperatures $T_{\rm H}$ and $T_{\rm C}$, with $\xi$ being the numerical constant

[TABLE]

the maximum being reached for $q_{0}\simeq 0.92$, which thus corresponds to the optimal thermal ground state population around which the cycle shall be performed, or in terms of the energy gap $\epsilon/T\simeq 2.4$. Note that this result accounts in taking $q_{0}\sim q_{1}$ which formally corresponds to performing a quasi-Otto cycles Abiuso and Perarnau , as it has been found for the exact optimal control of Carnot cycle in Ref. Cavina et al. (2017b)(with a specific dissipator) and Erdman et al. (2018).

III.2 Otto cycle

Again taking inspiration by the classical version translated in our setup, the Otto cycle is composed by two isoentropic (adiabatic) strokes alternated with two thermalizations (classically isochores). Considering the same qubit engine used for the description of the Carnot Cycle, the 4 steps can be summarized as in Fig. 3:

starting from an initial state $\hat{\rho}_{1}$, keeping fixed the gap $\epsilon_{1}$ the system is let thermalize in contact with the cold reservoir C;

• 2)

after isolating the system from the bath, a quench is performed taking $\epsilon_{1}\rightarrow\epsilon_{2}(>\epsilon_{1})$;

• 3)

while the gap is fixed, the system is let thermalize in contact with the hot reservoir H;

• 4)

a final quench restores $\epsilon_{2}\rightarrow\epsilon_{1}$.

Unless considering infinitesimal transformations where $\epsilon_{2}\simeq\epsilon_{1}$, it is clear that at variance with the Carnot cycle, in the Otto cycle the working medium ${\cal S}$ is always in a out-of-equilibrium state. Accordingly in this case the S-D approximation technique Cavina et al. (2017a) cannot be applied.

As for the Carnot cycle, the system exchange heat with the baths only during the steps 1) and 3). In particular exploiting the fact that now during the thermalization the Hamiltonian is kept constant we have

[TABLE]

where for $j={\rm H,C}$, $\Delta\hat{\rho}_{j}$ represent the increment experienced by the system density during the associated step. Equations (55) are valid for any Otto cycle, but we can specify them for our case $\hat{H}_{1,2}=\frac{\epsilon_{1,2}}{2}\hat{\sigma}^{z}$. If we also allow infinite time for the thermalization stages (ITT limit), the states at the end of the steps 1) and 4) are the thermal states $\hat{\Omega}^{({\rm C})}_{\epsilon_{1}}$ and $\hat{\Omega}^{({\rm H})}_{\epsilon_{2}}$, respectively, but in general they do not need to. In this case we have

[TABLE]

which replaced into (55) yields the identities

[TABLE]

where we use the upper index “$(0)$” to indicate that these are the heat exchanged in the IIT regime, and we used asymptotic ground state probabilities for the two isochores, defined as in Eq.(22)

[TABLE]

If we further assume the constraint

[TABLE]

$\Delta Q^{(0)}_{\rm H}$ turns out to be positive while $\Delta Q^{(0)}_{\rm C}$ is negative. The absorbed and released heat contributions can hence be identified as

[TABLE]

leading to an efficiency

[TABLE]

which thanks to (60) is smaller than the corresponding Carnot efficiency (32). Departing from the ITT regime, corrections can be computed analogously to what done for the Carnot cycle when considering non quasi-static cycles. Specifically we can write

[TABLE]

where now the $\Delta Q^{(1)}_{j}$s refer to first order corrections associated with the finite thermalization times, while the $\alpha_{j}$s are the associated ratios  $\Delta Q^{(1)}_{j}/\Delta Q^{(0)}_{j}$ analogous to those introduced in Eq. (35) for the S-D corrections of the Carnot cycle. In a similar way the power of the cycle can be expressed as in Eq. (36) yielding

[TABLE]

where $\tau_{\rm C}$ and $\tau_{\rm H}$ are the finite temporal durations of the steps

1. and 3) respectively, which for $\alpha_{j}=0$ gives the quasi-static result

[TABLE]

III.2.1 Exact FTT Otto Cycle

An application of the perturbative analysis (65) in the case of a general engine evolving under the action of the dissipation model (94), is presented in Appendix C . Due to the simplicity of the scheme however, this approach can be replaced by the exact finite-time solution of the problem, which we are going to present in the following.

Departing from the ITT regime the two thermalizations (isochores) of the Otto cycle become inevitably partial. To account for this effect, we represent the ground state populations of the working medium ${\cal S}$ at time $t$ after the beginning of the isochore with the $j$-th bath,

[TABLE]

where $p_{j}$ is the equilibrium probability (59) one would get in the strict ITT regime, and where $\Delta_{j}$ quantifies how out of equilibrium is the system at the beginning of the isochore. In this expression $f_{j}(t)$ is a function of $t$ that depends on the explicit details of the dynamics and which, by construction must fulfil the conditions $f_{j}(0)=1$ and $\lim_{t\rightarrow\infty}f_{j}(t)=0$ to ensure proper thermalization in the ITT regime. The parameters $\Delta_{\rm H}$ and $\Delta_{\rm C}$ are not completely independent and can be connected via the temporal durations, $\tau_{\rm C}$ and $\tau_{\rm H}$, of the two isochore. Indeed by invoking continuity conditions for the density matrix of $\mathcal{S}$ between the two isothermal strokes, we obtain

[TABLE]

which in particular imply

[TABLE]

From Eq. (55) it follows now that the relative heat exchanged during the isochore can now be expressed as

[TABLE]

which yields a power equal to

[TABLE]

where in the second line we used (69) and the expression (66) for the power of the cycle under ITT conditions (notice that as expected when $f_{j}(\tau_{j})=0$ then $P$ reduces to the value $P^{(0)}$ of Eq. (66)). This is the exact expression for $P$ which depends on the explicit form of $f_{j}$. It is worth observing that the first numerator $(\epsilon_{2}-\epsilon_{1})(p_{\rm C}-p_{\rm H})$ of Eq. (72) depends only on the model temperatures and gaps, hence fixing the efficiency it is possible to maximize the remaining independently, choosing the optimal length of the strokes.

IV Quantum Thermal Machines in the Non-Markovian regime

As we have explicitly discussed in Sec. II.2 the Markovian character of the system dynamics guarantees that the 2nd Law of thermodynamics is satisfied in the open quantum systems setting, i.e. the unavoidable loss of free energy of systems interacting with large baths. With this in mind, it is easy to realise that modelling the coupling of an engine with a non-Markovian bath may result as a pumping of free energy from the environment. For this reason, any acclaimed boost of performance in such a setting can be considered trivial, or even meaningless if not justified physically (e.g. using non-equilibrium baths). Therefore we choose to model the dynamics of the reservoir coupling in an overall Markovian framework, picturing an environment which contains some degrees of freedom who share correlations and interact with the working medium, such as to make its local dynamics non-Markovian; in this way we avoid any unjustified external resource draining. In such a set up indeed, Refs. Wilming et al. (2016); Lekscha et al. (2018); Perarnau-Llobet et al. (2018) show how this mechanism can only have detrimental effects from the point of view of quasi-static Thermodynamics: that is, without external free energy injections, the 2nd Law assures the Carnot efficiency is the maximal one. Nothing instead, has been stated from the point of view of FTT: even if lowering the maximal efficiency, there is still question on the effects on power and EMP, which are, from the practical point of view, much more interesting than pure maximal efficiency. Here we try to fill this gap, finding indeed that non-Markovian dynamics may have positive effects.

IV.1 The model

In order to account for non-Markovian effects we consider a modification of the set-up of Sec. II.1, with the one schematically sketched in Fig. 4 where both the hot reservoir H and the cold reservoir C include a local and a remote component. The first, represented by the qubit ancillary subsystems ${\cal A}_{\rm H}$ and ${\cal A}_{\rm C}$ of the figure, corresponds to degrees of freedom characterized by local Hamiltonian terms

[TABLE]

which are directly connected with $\mathcal{S}$ through dedicated coupling (energy exchanging) Hamiltonians which we assume to have the form

[TABLE]

$\hat{\mathcal{\sigma}}_{{\cal A}_{j}}^{\pm}$ and $\hat{\mathcal{\sigma}}^{\pm}$ indicating the lowering/raising operators of ${\cal A}_{j}$ and ${\cal S}$ respectively. The remote components of the baths, instead, are associated with standard local GKSL dissipators ($\mathcal{D}_{t}^{(j)}$ for $\mathcal{S}$ , and $\mathcal{D}_{{\cal A}_{\rm C}}$ and $\mathcal{D}_{{\cal A}_{\rm H}}$ for ${\cal A}_{\rm H}$ and ${\cal A}_{\rm H}$, respectively) inducing local thermalization toward their associated Gibbs states, i.e. the usual canonical state $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ of (3) for ${\cal S}$, and for $j={\rm H,C}$,

[TABLE]

for ${\cal A}_{j}$. Once more, we shall consider cyclic operations where the gap of the local Hamiltonian of ${\cal S}$ is externally modulated as in Eq. (20), and the system, at each given time, is selectively coupled to one and only one of the two baths. Accordingly we describe the dynamics of joint density matrix $\hat{\mathbf{R}}(t)$ the compound formed by ${\cal S}$, ${\cal A}_{\rm C}$, and ${\cal A}_{\rm H}$, in terms of a standard Markovian evolution which has the same form of Eq. (1), the non-Markovian character of the local dynamics of $\mathcal{S}$ being obtained instead by tracing away the ancillas, i.e. $\hat{\rho}(t)=\mbox{Tr}_{\cal A}[\hat{\mathbf{R}}(t)]$ follows trajectories that no longer exhibit the divisibility condition that instead is granted to $\hat{\mathbf{R}}(t)$ – see Appendix F. Furthermore, in order to simplify the analysis we shall also enforce the approximation, that on the time intervals ${\cal I}_{\rm C}$ (resp. ${\cal I}_{\rm H}$) during which the working medium is coupled with the cold bath ${\rm C}$ (resp. ${\rm H}$), the other ancilla ${\cal A}_{\rm H}$ (resp. ${\cal A}_{\rm C}$), that is temporarily decoupled from ${\cal S}$, relaxes to thermal equilibrium with the remote counterpart of ${\rm H}$ (resp. ${\rm C}$), i.e.

[TABLE]

with $\hat{\mathbf{R}}_{\rm C}(t)$ and $\hat{\mathbf{R}}_{\rm H}(t)$ that describe the reduced density matrix of of ${\cal S}{\cal A}_{\rm C}$ and ${\cal S}{\cal A}_{\rm H}$, respectively – the assumption being consistent with the first order S-D approximation, where ultimately one only needs to determine the quasi-static trajectories of the system rel . The evolution of $\hat{\mathbf{R}}_{j}(t)$ is finally expressed as

[TABLE]

with Hamiltonian

[TABLE]

and dissipator

[TABLE]

whose local contributions on ${\cal S}$ and on ${\cal A}_{j}$ will be assumed to have the simple form (94), i.e.

[TABLE]

where $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ and $\hat{\Omega}^{(j)}_{\hat{H}_{{\cal A}_{j}}}$ are the Gibbs states of (3) and (75) while $\Gamma_{j}$ and $\Gamma_{{\cal A}_{j}}$ are rates. As evident from the above expressions, we are assuming control on the energy gap of the working medium but not on the one of ${\cal A}_{j}$ which formally is just an element of the bath. We finally stress that $\gamma_{j}$ in Eq. (74) is the parameter defining the non-Markovianity of the model: this follows from the fact that the Markovian regime is recovered in the limit $\gamma_{j}\rightarrow 0$ (separable dynamics) and from the fact that, as shown in Appendix F, the non-Markovian measure by Breuer, Laine, Piilo Breuer et al. (2016) is monotonously increasing in $\gamma_{j}$. As we shall see in the next section a similar dependence can be observed on the optimized power output of a Carnot and Otto engine providing hence a clear indication of the fact that the non-Markovian character of the dynamics can be beneficial to these figures of merit.

IV.2 Non-Markovian Carnot cycle performance

In what follows we focus on the quasi-resonant case where the gap modulations of ${\cal S}$ on the interval ${\cal I}_{j}$ are such that the system is almost at resonance with ${\cal A}_{j}$, i.e.

[TABLE]

(note that the optimal Quasi-Otto trajectories found in Section III.1 are obtained in the limit of $\epsilon(t)$ being infinitesimally modulated). In this regime the stationary state of $\mathcal{\bf L}^{(j)}_{t}$ is approximatively equal to the tensor product of the individual thermal states associated with the two local dissipators, i.e.

[TABLE]

ensuring thermodynamic consistency of the model and being in agreement with Eq. (76). We hence identify the heat absorbed by ${\cal S}$ from the $j$-th bath as in Eq. (9) we obtain

[TABLE]

where in the first identity we used the fact that $\hat{\rho}(t)$ is the partial trace with respect to ${\cal A}_{j}$ of $\hat{\mathbf{R}}_{j}(t)$. We then expand this quantity as in Eq. (17) by invoking the S-D approximation $\hat{\mathbf{R}}_{j}(t)\simeq\hat{\mathbf{R}}^{(0)}_{j}(t)+\hat{\mathbf{R}}^{(1)}_{j}(t)+\cdots$, where $\hat{\mathbf{R}}^{(0)}_{j}(t)$ is the quasi-static solution which according to (83) is the state $\hat{\Omega}^{(j)}_{\hat{H}_{t}}\otimes\hat{\Omega}^{(j)}_{\hat{H}_{{\cal A}_{j}}}$, while $\hat{\mathbf{R}}^{(1)}_{j}(t)$ is the first order correction term which according to Eq. (16) we identify with the operator

[TABLE]

Due to the factorization of the fixed point (83) the zero-th term contribution $dQ^{(0)}_{j}(t)$ of $dQ_{j}(t)$ is still provided by the increment of the von Neumann entropy of the Gibbs states $\hat{\Omega}^{(j)}_{\hat{H}_{t}}$ as in Eq. (18). On the contrary $dQ^{(1)}_{j}(t)$ can still be cast in the form (38) where now ${\rho}^{(1)}_{00}(t)=\langle 0|\mbox{Tr}_{{\cal A}_{j}}[\hat{\mathbf{R}}^{(1)}_{j}(t)]|0\rangle$ that in the limit (82) can be expressed as in Eq. (39) with an S-D amplitude $A_{j}$ that can be found (see Appendix E.2) equal to

[TABLE]

where we introduce the quantities

[TABLE]

To evaluate the effect of non-Markovianity on the maximum power associated with a Carnot cycle we can then follow the same analysis we performed in Sec. III.1.1. In particular under symmetrization of the bath couplings and driving conditions (i.e. choosing $A_{\rm C}=A_{\rm H}=A$ and imposing $\epsilon(t)$ along the cold isotherm to be the time reversal of the one along the hot isotherm), we can directly use Eq. (53), which makes it clear that to get higher power performance we should target the low values of the $A_{j}$s.

First of all we notice that for $\gamma_{j}=0$ (i.e. $y=0$) correctly reduces to $A_{j}={1}/{\Gamma_{j}}$, which is the value (40) one would obtained in the Markovian limit in the presence of the dissipator (80). In the strong coupling limit $\gamma_{j}\gg\Gamma_{j}$ (i.e. $y_{j}\gg 1$), instead we get

[TABLE]

which gets smaller than the non-Markovian limit for values of $c_{j}$ above the critical threshold $\sqrt{2}-1\simeq 0.414$. Another important value is $c_{j}=2$ that determines the sign of the second addend in the parenthesis in the r.h.s. of Eq. (86) (the third addend being always positive). In fact we find that for $c_{j}<2$, $A_{j}$ attains its minimum value, smaller than the Markovian $1/\Gamma_{j}$, at

[TABLE]

otherwise the optimal value is infinite, in the sense of $A_{j}$ monotonously decreasing with $y_{j}$. These results on the dependence of $A_{j}$ on the model parameters are summed up in the Fig. 5, where we plot $A_{j}^{-1}$ in adimensional units, which thanks to Eq.s (50) and (53) corresponds to the maximum power attainable by the Carnot cycle, normalized to its Markovian value. In each case we see how the presence of the coupling $\gamma_{j}$, i.e. of non-Markovian effects leads to the an improvement of the maximum power $P_{\max}$ of the Carnot cycle.

IV.3 Non-Markovian Otto cycle performance

To study the performance of an Otto cycle for the non-Markovian model introduced, we can apply the exact power result (72), which in turn is characterized by the function $f_{j}$ describing the relaxation of the ground state during the two isochores ($j=\rm H,C$) (cf. Eq. (67)). In this case we restrict ourself to case $\Gamma_{{\cal{A}}_{j}}=\Gamma_{j}$ in resonance conditions between the system $\mathcal{S}$ and the ancillas, i.e. $\epsilon_{1}=E_{\rm C}$ and $\epsilon_{2}=E_{\rm H}$. As shown in Appendix E.1, under these conditions (82), the model allows for simple analytical solution of the form

[TABLE]

where to simplify the notation the time has been expressed in units of $\Gamma_{j}^{-1}$ and where $\kappa_{j}:=\sqrt{1-16(\gamma_{j}/\Gamma_{j})^{2}}$. As for the previous example we enforce symmetric conditions in which the couplings to the thermal baths have the same Lindbladian form and strength, i.e. $\Gamma_{\rm H}=\Gamma_{\rm C}$ and $\gamma_{\rm H}=\gamma_{\rm C}$, implying $f_{\rm C}(t)=f_{\rm H}(t)=f(t,\kappa)$. Under this assumption it is not difficult to prove (see Appendix G) that the maximum value of the power obtainable from (72) is found on the bisector $\tau_{\rm C}=\tau_{\rm H}=\tau$. Specifically with this choice we get

[TABLE]

which for each value of $\kappa(y)$ has a maximum for a finite value of the duration $\bar{\tau}(y)$. In Fig. 6 we plot the obtained maximum as a function of $y=\gamma/\Gamma$, normalized to the $\gamma=0$ case, for which Eq. (91) is easily seen to have a maximum equal to $\frac{\Gamma}{4}(E_{\rm H}-E_{\rm C})(p_{\rm C}-p_{\rm H})$. As in the case of the Carnot cycle we see once more that increasing the strength of the non-Markovian coupling parameter $\gamma$ the power of the Otto engine dramatically increases.

IV.4 Free-energy analysis

Ruling out the possibility of using non-Markovian effects to improve the efficiency in the model  Wilming et al. (2016); Lekscha et al. (2018); Perarnau-Llobet et al. (2018), the power boost we reported above for the Carnot and Otto cycle can only be seen as a consequence of the latter in the reduction of thermalization timescales. We show here an argument to explain why it happens and how it is related to the non-Markovian building of correlations between the engine $\mathcal{S}$ and the $j$-th bath.

When attaching $\mathcal{S}$ to the $j$-th bath, the former is out of equilibrium, while according to Eq. (76) $\mathcal{A}_{j}$ is already thermalized: we can hence describe their initial state as $\hat{\rho}(0)\otimes\hat{\Omega}^{(j)}_{\hat{H}_{\mathcal{A}_{j}}}$. A possible way to quantify the rapidity of $\mathcal{S}$ thermalizing is to compute the relative entropy $S(\hat{\rho}(t)\parallel\hat{\Omega}^{(j)}_{\hat{H}_{t}})$ and see how fast it diminishes. We recall that according to (12) this quantity also measures the excess of free energy present in the system from the corresponding Gibbs state at temperature $T_{j}$. Indicating hence with $\hat{\mathbf{R}}_{j}(t)$ the joint state of ${\cal S}$ and ${\cal A}_{j}$ at time $t$ we notice that

[TABLE]

where $I(\mathcal{S}:\mathcal{A}_{j})\geq 0$ is the mutual information Holevo (2012) between the two systems at time $t$ (in the above derivation we explicitly use the fact that the mean value of the interaction term stays null due to the quasi-resonant conditions assumptions). Hence the variation of the free energy of $\mathcal{S}$ can be expressed as

[TABLE]

Now we observe that the last two terms provide a negative contributions to $\Delta S(\hat{\rho}(t)\parallel\hat{\Omega}^{(j)}_{\hat{H}_{t}})$. Indeed being the relative entropy and mutual information positive definite, and being the initial state of ${\cal S}$ and ${\cal A}_{j}$ factorized by hypothesis, they are initially null, meaning that they can only increase with time, thus bringing negative contribution to the r.h.s. of (IV.4), that is faster thermalization of ${\cal S}$. For a comparison, in case $\mathcal{S}$ and $\mathcal{A}_{j}$ do not interact (that is $\gamma=0$ in our model) mutual information would remain zero and $\hat{\rho}_{{\cal A}_{j}}$ would stay thermal making these two terms exactly null. In support of the above analysis we report in Fig. 7 an example of an evolution of the above quantities. It is possible to see how the free energy of $\mathcal{S}$ decreases faster than the total free energy, due to ”suction of free energy” by $\mathcal{A}_{j}$ and the correlations $I(\mathcal{S}:\mathcal{A}_{j})$ building.

V Conclusions

The purpose of the work is to assess the effects of non-Markovian dynamics on quantum thermal machines. We considered a simple class of models which allows shared correlation between the system and some degrees of belonging to the baths, preserving the global evolution as Markovian. This avoids resource pumping from the baths, and cannot induce advantages from the quasi-static point of view. Exploiting the S-D technique we studied the thermodynamic performance for a (finite-time) Carnot cycle: results indicate that the maximum power can indeed be boosted by the presence of this non-Markovian mechanism. Exact results obtained by studying Otto cycles confirm this trend. Noting that in general the S-D amplitude $A$ is related to the relaxation time of the system, we are naturally led to interpret this positive effect as an acceleration of the thermalization timescale of ${\cal S}$, in presence of its possible interaction with the local components of the baths (cf. Fig.4). Again this is intuitively reasonable, having $\mathcal{S}$ a new channel of thermalization which passes through $\mathcal{A}_{j}$, and we showed explicitly in Sec. IV.4 how this effect is related to the non-Markovian feature of the baths building correlations with the working medium.

Preliminary to the results, we also showed how to optimize the control for a 2-level engine performing a Carnot cycle (in the low-dissipation regime) or an Otto cycle (exactly).

A possible extension of the research would be to generalize these results for systems beyond qubits, with variable number of levels or in the geometrical picture introduced in Scandi and Perarnau-Llobet (2018); Abiuso and Perarnau .

Acknowledgements

The authors wish to thank V.Cavina, P.Erdman, A.Mari for useful discussions. P.A. is supported by the Spanish MINECO (QIBEQI FIS2016-80773-P, and Severo Ochoa SEV-2015-0522), Generalitat de Catalunya (SGR1381 and CERCA Programme), Fundacio Privada Cellex.

Appendix A The dissipators

Here we review few examples of dissipators $\mathcal{D}_{t}^{(j)}$ that obey the constraints (2) and (4).

The first, and simplest of such models, is provided by the super-operator Cavina et al. (2017b)

[TABLE]

with $\Gamma_{j}>0$ constant, which do not need any specification of the system Hamiltonian.

Assuming instead the Hamiltonian of ${\cal S}$ to be $\hat{H}_{t}={\epsilon}(t)\hat{\mathcal{\sigma}}^{z}/2$ (see Eq. (20)), another example is provided by the dissipator

[TABLE]

where $\hat{\sigma}_{+}$ and $\hat{\sigma}_{-}(=\hat{\sigma}_{+}^{\dagger})$ are, respectively, the raising and lowering operators of ${\cal S}$, $[\cdots,\cdots]_{+}$ is the anti-commutator, which exhibit the functional dependence (2) upon $\hat{H}_{t}$ through the rates $\Gamma^{(j)}_{\pm}(\epsilon)$ fulfilling the detailed balance equation condition

[TABLE]

which ensures (4). In particular taking

[TABLE]

with $\Gamma\geq 0$ and

[TABLE]

equation (95) can be used to describe the interaction of ${\cal S}$ with a Fermionic bath. Instead taking

[TABLE]

with $\Gamma\geq 0$ and

[TABLE]

it describes the interaction of ${\cal S}$ with a Bosonic bath.

Appendix B S-D approximation implies low dissipation

A virtue of the S-D approximation is that it provides a formal justification of the low-dissipation (L-D) assumption Esposito et al. (2010) which is typically introduced in FTT analysis as a phenomenological working hypothesis. To see this let us start observing that in the S-D theory, at the lowest order of the pertubative expansion (15) the von Neumann entropy of the density matrix $\hat{\rho}(t)$ can be expressed as

[TABLE]

where in the last step we used the fact that the term $\hat{\rho}^{(1)}(t)$ is traceless, i.e. $\text{Tr}[\hat{\rho}^{(1)}(t)]=0$, and the fact that $\hat{\rho}^{(0)}(t)$ is the instantaneous Gibbs state (3), i.e. $\hat{\rho}^{(0)}(t)=\hat{\Omega}^{(j)}_{\hat{H}_{t}}$. A close inspection reveals that the second contribution of $S(t)$ corresponds to the first order correction to the internal energy of the system defined in Eq. (7), i.e. $E_{1}(t):=\text{Tr}[\hat{\rho}^{(1)}(t)\hat{H}_{t}]$, allowing us to cast (99) as

[TABLE]

The temporal increment of this quantity can hence be computed as

[TABLE]

where we used Eq. (18) and wrote $dE^{(1)}(t)$ in terms of a work and heat contribution, i.e. $dE^{(1)}(t)=dQ^{(1)}(t)+dW^{(1)}(t)$ (first thermodynamics principle) with $dQ^{(1)}(t)$ as in Eq. (19) and

[TABLE]

Grouping together all the heat contributions we can hence finally write

[TABLE]

where $dS^{(irr)}_{j}(t):=dS(t)-\beta_{j}dQ(t)$ is the irreversible entropy production increment which quantifies the differences between information transfer rates and the heat transfer rate in the system. When integrated over a finite time interval $\tau_{j}$, Eq. (103) provides an estimation of the associated finite irreversible entropy production $\Delta S^{(irr)}_{j}$. In FTT under L-D assumption this term is postulated to be expressed as inversely proportional to $\tau_{j}$ via a constant term $\Sigma_{j}$ which only depends on the coupling constants to the bath, and the cycle endpoints, i.e. Esposito et al. (2010); Ma et al. (2018)

[TABLE]

Now a $1/\tau_{j}$ scaling as in Eq. (104) is exactly what one naturally get by computing $\Delta S^{(irr)}_{j}$ via direct integration of (103) due to the fact that in the S-D expansion the $\hat{\rho}^{(1)}(t)$ term has an explicit linear dependence upon $1/\tau_{j}$, whilst $\hat{\rho}^{(0)}(t)$ and the associated instantaneous Hamiltonian $\hat{H}_{t}$ are independent from such parameter Cavina et al. (2017a, 2018). According to this observation, on one side we can hence say that S-D provide a natural framework for discussing L-D assumption. On the other side instead we can conclude that the general results derived under FTT assumption Esposito et al. (2010); Ma et al. (2018) must apply in the characterization of system driven under S-D approximation, at least at the first order of the perturbative analysis. In particular it is not difficult to see that the condition $\alpha_{\rm C}=\alpha_{\rm H}$ we require in Sec. III.1 corresponds to set $\Sigma_{\rm C}=\Sigma_{\rm H}$, which in turn implies that the Curzon-Ahlborn efficiency is the EMP in this regime Esposito et al. (2010). Getting rid of this symmetry simply means to explore the different ratios $\Sigma_{\rm H}/\Sigma_{\rm C}$ and the relative results Esposito et al. (2010); Ma et al. (2018) are valid.

Appendix C Otto cycle beyond the ITT limit

To evaluate the correction terms appearing in Eqs. (64) and (65) for a generic engine we assume the dissipation model of Eq. (94). Let then indicate with $\hat{\rho}^{(in)}_{j,k}$ and $\hat{\rho}^{(out)}_{j,k}$ the states of ${\cal S}$ at the beginning and at the end of the time intervals ${\cal I}_{j}$ of the $k$-th Otto cycle. Due to the presence of the quenches at the steps 2) and 4), they must be related as follows

[TABLE]

meaning that input state of the $k$-th interval ${\cal I}_{\rm H}$ coincides with the output state of the $k$-th interval ${\cal I}_{\rm C}$, while the input of the $k$-th interval ${\cal I}_{\rm C}$ with the output of the $k$-th interval ${\cal I}_{\rm H}$. By direct integration of the ME (1) we get

[TABLE]

which, with the help of (105) can be equivalently cast in the following recursive expressions

[TABLE]

Now in the ITT limit $\tau_{j}\rightarrow\infty$ these yields $\hat{\rho}^{(out)}_{{\rm C},k}=\hat{\Omega}^{({\rm C})}_{\epsilon_{1}}$ and $\hat{\rho}^{(out)}_{{\rm H},k}=\hat{\Omega}^{({\rm H})}_{\epsilon_{2}}$ for all $k$, leading to (111) via Eq. (105). For finite $\tau_{j}$ instead, keeping only the most relevant order, we obtain

[TABLE]

for all $k$, which, exploiting once more (105), gives

[TABLE]

Inserting this into Eq. (55) we can express the first order corrections $\Delta Q^{(1)}_{j}$s as

[TABLE]

hence obtaining

[TABLE]

for $j={\rm H,C}$. From Eq. (64) then follows that the efficiency remains un-effected by the ITT corrections, i.e. $\eta=\eta_{\text{o}}$, while according to Eq. (66) the power becomes

[TABLE]

Appendix D Optimal protocol shape for the Quantum Carnot cycle

In this section we solve the minimization of the functional ${\cal F}[q(x)]$ of (51) (hereby ${\cal F}[q]$ for short) under the constraints

[TABLE]

which, according to Eq. (50) allow us to optimize the power production on the Quantum Carnot cycle. First of all we notice that it can be equivalently expressed as

[TABLE]

where the modulus has been replaced by a minus sign, due to the fact that integrand is guaranteed to be non-positive for all the allowed choices of the function $q(x)$ (same argument we used in Eq. (41) to establish the non positivity of $\Delta Q_{j}^{(1)}$).

For this purpose we consider the variation of the functional (51) under a small variation of the control $q\rightarrow q+\delta q$,

[TABLE]

where the last identity was obtained by integration by parts using the constraints 115. Imposing the latter to nullify under arbitrary variation we can then obtain the differential equation

[TABLE]

which can be solved using separation of variables and the substitution $q^{\prime}=q-\frac{1}{2}$, leading to optimal solutions of the form

[TABLE]

This class of solutions is parametrised by the two values $\{\omega,\varphi\}$ and is in general incompatible with the constraints 115 given at the extrema: this is a typical issue one meets in variational problems performed on given sets of functions that are not topologically closed; that is, it is possible to construct a sequence of functions $q_{k}(x)$ which decrease the functional toward an infimum which, however, is reached only for a function $\lim_{k\rightarrow\infty}q_{k}=\bar{q}$ that is outside the initial function space. We can build the sequence $q_{k}$ by simply stringing smoothly $\bar{q}(0)$ to $\bar{q}(\varepsilon)$ and $\bar{q}(1-\varepsilon)$ to $\bar{q}(1)$ for small $\varepsilon=\frac{1}{k}$,

[TABLE]

with sufficiently smooth functions $s_{k}$ such that $s_{k}(0)=\bar{q}(0),\ s_{k}(\frac{1}{k})=\bar{q}(\frac{1}{k}),\ \dot{s}_{k}(0)=0,\ \dot{s}_{k}(\frac{1}{k})=\dot{\bar{q}}(\frac{1}{k})$ and similarly for $x=1-\frac{1}{k},1$. The correction to the optimal functional $\delta_{k}{\cal F}={\cal F}[{q_{k}}]-{\cal F}[{\bar{q}}]$ will then be given from the contribution near the border $[0,\frac{1}{k}]$

[TABLE]

and the analogous term for $[1-\frac{1}{k},1]$. Using that $\dot{s}_{k}(\frac{1}{k})-\dot{s}_{k}(0)=\dot{\bar{q}}(\frac{1}{k})$ we can take the limit to obtain, adding the $[1-\frac{1}{k},1]$ contribution,

[TABLE]

and thus the optimal value

[TABLE]

It is however important to stress that the sequence $q_{k}$ will eventually break the S-D approximation, having a high second derivative near 0 and 1. Hence one should ”stop” to a $k$ which is not too big to achieve these approximate results. Numerical plots support the achievability of the limit.

Quasi-Otto limit.

When allowing the control to vary the initial and final point $\{q(0),q(1)\}$, numerical plots show that the choice that maximizes the value of the power (50), which is equal to

[TABLE]

is obtained in the limit of them being the same $q(1)-q(0)=\varepsilon\rightarrow 0$. In this limit $\bar{q}(x)$ (119) is essentially a line with $\dot{\bar{q}}\sim\frac{\bar{q}(1)-\bar{q}(0)}{1}-0=\varepsilon\;,$ hence the numerator of (124) is

[TABLE]

For the denominator ${\cal F}_{min}$ the contribution ${\cal F}[{\bar{q}}]$ nullifies while $\bar{\delta}{\cal F}$ can be estimated from (122) as

[TABLE]

where we use $\dot{q}=\varepsilon$ and the derivative of $\ln(q/(1-q))$, which is $1/(q(1-q))$. We can then write (124) in this limit

[TABLE]

It is possible plot this function to find the optimal value of $\bar{q}$ around which to perform the optimal control, as in Figure 8, where is evident that $q_{optimal}\sim 0.92$ and $P_{max}\sim 0.11\frac{(\sqrt{T_{H}}-\sqrt{T_{C}})^{2}}{A}$.

Appendix E Dynamical solution of the non-Markovian model

In this Appendix we show exact and approximate (Slow-Driving) solutions of the non-Markovian model of Section IV.1. We remind that we are interested in the analysis of the dynamics of the qubit $\mathcal{S}$ when coupled to one of the two baths, that without loss of generality can be considered to be the cold one, so that the state of the system and the ancillary qubit ${\cal{A}}_{\rm C}$ of the bath can be described by the density matrix $\hat{\mathbb{R}}_{\rm{C}}$, as well al the local states of $\mathcal{S}$ ($\hat{\rho}:=\text{Tr}_{\cal{A}_{\rm C}}[\hat{\mathbb{R}}_{\rm{C}}]$) and $\cal{A}_{\rm C}$ ($\hat{\rho}_{\rm C}:=\text{Tr}_{\cal S}[\hat{\mathbb{R}}_{\rm{C}}]$). As pictured in Fig.4 we remind the local Hamiltonians and coupling interaction

[TABLE]

as well as the thermalizing dissipators

[TABLE]

so that in the interaction picture the dynamical equation is

[TABLE]

with $\delta=\epsilon-E_{C}$.

Introducing the thermal ground state probabilities and a time dependent phase

[TABLE]

we can solve equation (130) by writing it in the computational basis $\hat{\mathbb{R}}_{\rm{C}}\equiv\sum_{\alpha,\beta,\mu,\nu=0,1}{\rho}_{\alpha\mu\beta\nu}|\alpha\rangle\!\langle\beta|_{{}_{\mathcal{A}_{\rm C}}}\otimes|\mu\rangle\!\langle\nu|_{{}_{\cal S}}\;.$ We get, in matrix form,

[TABLE]

where the inferior triangular part has been omitted to improve readability and can be filled by just noting $\hat{\mathbb{R}}_{\rm{C}}$ is hermitian. For each matrix element the 3 different lines represent the contributions from the ${\mathcal{A}_{\rm C}}$ dissipator ($\propto\Gamma_{{\cal A}_{C}}$), the $\mathcal{S}$ dissipator ($\propto\Gamma_{\rm C}$), and the Hamiltonian exchange ($\propto\gamma_{\rm C}$). Looking at the equation we can note that the time evolution generator is a sparse super-operator, which couples separately different subsets of components, namely the ones highlighted here with different colors

[TABLE]

Three different sets of equations can be then solved separately, but we will be interested in the Thermodynamics of the system, hence mainly the highlighted blue subset, because it contains the populations which determine thermodynamic variables (namely, the eigenstates of $\hat{H}_{t}+\hat{H}_{\mathcal{A}_{\rm C}}$). We can represent it as the vector

[TABLE]

where $q_{ab}={\rho}_{abab}$ is the population of $\mathcal{S}$ in the state $b$ and ${\mathcal{A}_{\rm C}}$ in $a$, while $k\equiv{\rho}_{0110}$ is the coherence between the states $\ket{01}$ and $\ket{10}$ which are the ones interacting by the exchange Hamiltonian $V_{\rm C}$.

E.1 Resonant case $(\epsilon=E_{\rm C})$

Consider the instance in which the gaps are fixed equal $\epsilon=E_{\rm C}=E$; in this case the algebra has some simplifications; indeed

[TABLE]

so that the value of the interaction Hamiltonian is conserved in absence of the dissipative dynamics (or decreases exponentially, see below). It is easy to check that the (only) stationary state is $\hat{\mathbb{R}}_{\rm{C}}=\hat{\Omega}^{(C)}_{\hat{H}_{t}}\otimes\hat{\Omega}^{(C)}_{\hat{H}_{{\cal A}_{\rm C}}}$ (or $\hat{\Omega}_{\rm C}\otimes\hat{\Omega}_{\rm C}$ for simplicity).Having the gap equal we can call

[TABLE]

We will also write $\mathcal{L}$, with a small abuse of notation, to indicate the Lindblad generator of the dynamics restricted to the different subsets of components. The equation (132) for the vector (133) can be written in this special case as

[TABLE]

Note that the real part of the coherence $\Re(k)$ satisfies  $\Re(\dot{k})+\Gamma\Re(k)=0$ , hence it is decoupled from the rest and it just dies exponentially111Note that $\langle V_{\rm C}\rangle=\text{Tr}[\hat{\mathbb{R}}_{\rm{C}}\gamma_{\rm C}(\sigma^{+}\otimes\sigma_{\mathcal{A}_{\rm C}}^{-}+h.c.)]=\gamma_{\rm C}(\rho_{{}_{0110}}+\rho_{{}_{1001}}=\gamma_{\rm C}2\Re(k)$ which is then always decreasing. This means for initial condition given by a product state $\rho\otimes\rho_{\rm C}$, $\Re(k)$ is constantly null, which in turn implies that the switch-on/switch-off work done to attach the system $\mathcal{S}$ to the baths is null and can be safely neglected in the performance analysis. $\sim e^{-\Gamma t}$. Calling $\Im(k)\equiv I$ the system can be thus be written

[TABLE]

the Lindblad generator being $\mathcal{L}=\tilde{\mathcal{L}}-(\Gamma_{\rm C}+\Gamma_{{\mathcal{A}_{\rm C}}})$ and

[TABLE]

which in the $\Gamma_{{\cal A}_{C}}=\Gamma_{\rm C}=\Gamma$ case becomes ($\gamma_{\rm C}^{\prime}=\gamma_{\rm C}/\Gamma$)

[TABLE]

To solve the dynamics one can find eigenvalues and eigenvectors of such a matrix. In the $\Gamma_{\rm C}=\Gamma_{{\cal A}_{C}}$ case the particular symmetry of the problem is reflected in the tractable form of the eigensystem of $\mathcal{L}$, which is (subtracting already $-2\Gamma$ to all eigenvalues and expressing in units of $\Gamma$)

[TABLE]

We can then solve completely the dynamics for an initial state of the form $\hat{\mathbb{R}}_{\rm{C}}(t_{0})=\hat{\rho}\otimes\hat{\Omega}_{\rm C}$, that is out of equilibrium on $\mathcal{S}$ and thermal on ${\mathcal{A}_{\rm C}}$, as requested by our model.

Suppose for the moment that also $\hat{\rho}$ is diagonal, that is

[TABLE]

We write $a=p_{0}+\Delta$ to quantify how much $\hat{\rho}$ is out of equilibrium. We decompose (141) as a combination of the eigenvectors (140), in order to write the solution which will be

[TABLE]

Summing the first two components we can obtain the time-dependent ground state population of $\mathcal{S}$ ​ that is, calling $\kappa_{\rm C}=\sqrt{1-16\gamma_{\rm C}^{\prime 2}}\ ,$

[TABLE]

having defined

[TABLE]

E.2 Non-resonant case $(\epsilon(t)\neq E_{\rm C})$ - Slow-Driving

In case the two qubits are not resonant equation (132) for the vector (133) takes the general form

[TABLE]

We note that the 2nd and 3rd equation here can be rewritten using222Remember $\phi(t)\equiv e^{i\delta t}$. $\hat{k}=\phi k$, which satisfies

[TABLE]

In this way we can write, calling $\Im(k)\equiv I$ and $\Re(k)\equiv R$,

[TABLE]

with $\mathcal{L}=\tilde{\mathcal{L}}-(\Gamma_{\rm C}+\Gamma_{{\cal A}_{C}})\mathds{1}\ ,$

[TABLE]

The null eigenvector of $\mathcal{L}$ (i.e. the stationary state $\rho^{(0)}$) is not in general simply the thermal state $\Omega_{\rm C}\otimes\Omega_{\rm C}$, but it reduces to it in the limit

[TABLE]

At first order333Note that $(p_{s}-p_{\rm C})\equiv\Delta_{p}$ is $\sim O(\delta)$. in $\delta$ we find

[TABLE]

the normalization being $N=(\Gamma_{\rm C}\Gamma_{{\cal A}_{C}}+2\gamma_{\rm C}^{2})(\Gamma_{\rm C}+\Gamma_{{\cal A}_{C}})\ .$ Note that in both approximations the real part of the coherence $\Re(\rho_{{}_{0110}})=[\vec{\rho}_{0}]_{6}=0$ is null; this allows us to neglect work contribution in the contacts and detachments from the baths, as in the resonant case444$\langle V_{\rm C}\rangle=\text{Tr}[\hat{\mathbb{R}}_{\rm{C}}\gamma_{\rm C}(\sigma^{+}\otimes\sigma_{\mathcal{A}_{\rm C}}^{-}+h.c.)]=\gamma_{\rm C}(\rho_{{}_{0110}}+c.c.)=2\gamma_{\rm C}\Re(k)$ which is therefore null at the beginning and ending of each isothermal stroke, when $\hat{\mathbb{R}}_{\rm{C}}(t)=\hat{\rho}_{0}(t)$..

Following the approach described in Section II.3 we can now compute the first order correction in Slow-Driving to the dynamics. Looking at the formal solution of the S-D technique (16) we need for the computation:

• •

the quasi-static solution $\hat{\rho}^{(0)}$ found in (150),

• •

the dynamics generator $\mathcal{L}$ we wrote explicitly,

• •

the projector on the null-trace subspace $\mathcal{P}$.

This last operator is easily found. The trace of the state is given from the sum of the 4 populations

[TABLE]

In order to project on the null-trace subspace we have to subtract to each population $\text{Tr}[\rho]/4$, that is $q_{ij}\rightarrow q_{ij}-\frac{1}{4}\sum_{a,b}q_{ab}$, while the coherences stay unchanged. This can be written as

[TABLE]

Now we have all the ingredients to compute (with the help of Wolfram Mathematica) the first correction \vec{\rho}_{1}=(\mathcal{LP})^{-1}\dot{\vec{\rho}}_{0}\. The resulting expression is too complicated to be reported here, however, it is possible to compute the correction to the ground state population of $\mathcal{S}$ ​, $q_{00}^{(1)}+q_{10}^{(1)}$, which is in the form

[TABLE]

with an amplitude $A_{\rm C}$ that admits a closed but unfortunately still very convoluted general expression, which for the sake of readability we do not report here in its full extension. Nevertheless, in the resonance ($\epsilon=E_{\rm C}$) limit we considered for our model we find

[TABLE]

Appendix F Non-Markovian character of the dynamics

Here we show that the model of Sec. IV has an explicit non-Markovian character which depends on the non zero value of the parameters $\gamma_{j}$ that gauge the coupling between $\mathcal{S}$ and the ancillas ${\cal A}_{j}$. For this task, given two input states $\hat{\rho}^{(1)}(0)$, $\hat{\rho}^{(2)}(0)$ of ${\cal S}$ and $\hat{\rho}^{(1)}(t)$, $\hat{\rho}^{(2)}(t)$ their corresponding dynamical evolutions under the action of the model, we consider the information-backflow BLP quantity (Breuer et al., 2016)

[TABLE]

with $\Theta$ being the Heaviside function that restrict the domain of integration to the one where the integrand is positive, and where $D(\hat{\rho}^{(1)}(t),\hat{\rho}^{(2)}(t)):=\frac{1}{2}||\hat{\rho}^{(1)}(t)-\hat{\rho}^{(2)}(t)||_{1}$ is the trace distance Holevo (2012). As discussed in Ref. (Breuer et al., 2016) value of $\mathcal{N}_{BLP}(\hat{\rho}^{(1)}(0),\hat{\rho}^{(2)}(0))$ greater then zero would imply non-Markovian character of the dynamics.

In our case, focusing only at $\hat{\rho}^{(1)}(0)$, $\hat{\rho}^{(2)}(0)$ having no coherence terms, we can simplify the analysis exploiting the fact that $D(\rho^{(1)}(t),\rho^{(2)}(t))=|p^{(1)}(t)-p^{(2)}(t)|$, where for $j=1,2$, $p^{(j)}(t)$ is the ground state population of the $\hat{\rho}^{(j)}(t)$. On-resonance ($\epsilon=E_{j}$) the solutions given in Sec. E.1 yields, in adimensional units ($\Gamma_{j}=1$),

[TABLE]

where $\kappa_{j}=\sqrt{1-16(\gamma_{j}/\Gamma_{j})^{2}}$ . Replacing this into (155) and performing the integration numerically we obtain the results reported in Fig.  9 as function of $y_{j}=\gamma_{j}/\Gamma_{j}$. As expected, the non-Markovianity is monotonously increasing with $\gamma_{j}$. Also, we find a threshold value under which this particular non-Markovianity witness is null.

Appendix G Symmetric Otto cycle has maximum power for $\tau_{\rm H}=\tau_{\rm C}$

In this appendix we prove that the power expressed by Eq. (72) is maximized, in case the coupling to the two baths is symmetric (i.e. $f_{\rm H}(t)=f_{\rm C}(t):=f(t)$ in Eq. (72)), by choosing the time durations $\tau_{\rm H}=\tau_{\rm C}$ equal. In fact under this assumption the power can be written as

[TABLE]

We show that when $\tau_{\rm H}\neq\tau_{\rm C}$ at least one between $C(\tau_{\rm C},\tau_{\rm C})$ and $C(\tau_{\rm H},\tau_{\rm H})$ is greater than $C(\tau_{\rm C},\tau_{\rm H})$, meaning that $\{\tau_{\rm H},\tau_{\rm C}\}$ would be outperformed by one of the two choices. To prove it we demonstrate that $C(\tau_{\rm C},\tau_{\rm H})\leq\sqrt{C(\tau_{\rm C},\tau_{\rm C})C(\tau_{\rm H},\tau_{\rm H})}$, which implies555Note that $C\geq 0$, by the definition (157) and $0\leq f\leq 1$. the thesis. This is equivalent to verify the following inequality holds

[TABLE]

which is true by noting the numerator is the same and on the denominator by direct inspection

[TABLE]

[TABLE]