Strong coupling corrections in quantum thermodynamics

TL;DR

This paper investigates how strong coupling between quantum systems and thermal baths affects thermodynamic laws, providing corrections to the second law, bounds on heat engine power, and illustrating these effects with quantum Brownian motion.

Contribution

It introduces strong-coupling corrections to the second law of thermodynamics and bounds on heat engine power, extending thermodynamics to the strong coupling regime.

Findings

Strong-coupling corrections to the second law are derived.

Bounds on power enhancement due to strong coupling are established.

Results are exemplified using non-Markovian quantum Brownian motion.

Abstract

Quantum systems strongly coupled to many-body systems equilibrate to the reduced state of a global thermal state, deviating from the local thermal state of the system as it occurs in the weak-coupling limit. Taking this insight as a starting point, we study the thermodynamics of systems strongly coupled to thermal baths. First, we provide strong-coupling corrections to the second law applicable to general systems in three of its different readings: As a statement of maximal extractable work, on heat dissipation, and bound to the Carnot efficiency. These corrections become relevant for small quantum systems and always vanish in first order in the interaction strength. We then move to the question of power of heat engines, obtaining a bound on the power enhancement due to strong coupling. Our results are exemplified on the paradigmatic situation of non-Markovian quantum Brownian motion.

Equations

$$W_{\text{on}}^{(i)}=\mathrm{Tr}({\rho }^{(i)}V)=-W_{\text{off}}^{(i)}.$$

$$W^{(i)}=\mathrm{Tr}({\rho }_S^{(i)}(H_S^{(i)}-H_S^{(i+1)}))$$

$${\rho }_S^{(i+1)}={\mathrm{Tr}}_B({\omega }_{\beta }(H^{(i+1)})),$$

$$W=W^{(\mathrm{weak})}-\Delta F^{(\mathrm{res})}-\Delta F^{(\mathrm{irr})},$$

$$\Delta F^{(\mathrm{irr})}$$

$$\Delta F^{(\mathrm{res})}$$

$${\mathrm{Tr}}_B\left({\omega }_{\beta }(H^{(1)})\right)={\rho }_S,$$

$${\mathrm{Tr}}_B\left({\omega }_{\beta }(H^{(N)})\right)=\frac{{\mathrm{Tr}}_B\left({\omega }_{\beta }(H^{(N)})(H_{H^{(N)},\beta }^{(0)}-H^{(N)})\right)}{\mathrm{Tr}\left((H_{H^{(N)},\beta }^{(0)}-H^{(N)}){\omega }_{\beta }(H^{(N)})\right)},$$

$$H_S^{(1)}={\tilde{H}}_S-g{\mathrm{Tr}}_B(\omega (H_B)V)+O(g^2),$$

$$H_S^{(N)}=H_S-g{\mathrm{Tr}}_B(\omega (H_B)V)+O(g^2),$$

$$\Delta F_{\min }^{(\mathrm{irr})/(\mathrm{res})}=\frac{\beta g^2}{2}{\mathrm{cov}}_{{\omega }_{\beta }({\tilde{H}}^{(0)}/H^{(0)})}(\tilde{V},\tilde{V})+O(g^3).$$

$$Q=T\Delta \mathrm{S}-(\Delta F_B^{(\mathrm{res})}+TI({\omega }^{(N)};S:B)+\Delta F^{(\mathrm{irr})})$$

$$Q=T\Delta \mathrm{S}-K_q{g}^2+O(g^3),$$

$${\omega }_{{\beta }_h}(H_S^{(B)/(A)})={\omega }_{{\beta }_c}(H_S^{(C)/(D)}),$$

$${\mathrm{Tr}}_{B_h}{\omega }_{{\beta }_h}(H^{(B)/(A)})={\mathrm{Tr}}_{B_c}{\omega }_{{\beta }_c}(H^{(C)/(D)}),$$

$$\eta ={\eta }^{\mathrm{C}}-g^2\frac{T_c}{T_h}\left(\frac{K_q^{(h)}}{Q_h^{(\mathrm{weak})}}+\frac{K_q^{(c)}}{Q_c^{(\mathrm{weak})}}\right)+O(g^3),$$

$$\Delta F_{\min }^{(\mathrm{res})/(\mathrm{irr})}\le 2\Vert V\Vert ,$$

$$TI({\omega }^{(N)};S:B)\le 2\Vert V\Vert ,$$

$$V=\sum _{j=1}^{\alpha n^2}{h}_j,$$

$$\Vert V\Vert \le \sum _{j=1}^{\alpha n^2}\Vert h_j\Vert \le \alpha n^2\underset{j}{\max }\Vert h_j\Vert .$$

$$P:=\frac{W}{\Delta t}\le \frac{g{r}_c\eta (g)}{1-\eta +r_c/r_h}<g{r}_h\eta (g),$$

$$P(g)\le g{r}_c{\eta }^{\mathrm{C}}/(1-{\eta }^{\mathrm{C}}+r_c/r_h)-O(g^3).$$

$$H=H_S+gV+H_B+H_L,$$

$${\omega }_{\beta }(H)=\exp (-\beta H)/\mathrm{Tr}(\exp (-\beta H)).$$

$$\Vert X\Vert =\underset{|\psi ⟩}{\sup }⟨\psi |X|\psi ⟩,$$

$$\displaystyle\frac{d}{dg}\mathrm{Tr}(f(A+gB))\bigg{|}_{g=0}=\mathrm{Tr}(Bf^{\prime}(A)),$$

$$\displaystyle\frac{d}{dg}F(\omega_{\beta}(H),H)\bigg{|}_{g=0}=\frac{-T}{\mathcal{Z}_{0}}\frac{d\mathcal{Z}}{dg}\bigg{|}_{g=0}=\mathrm{Tr}(V\omega_{\beta}(H_{0})).$$

$$\frac{\mathrm{d}}{\mathrm{d}g}{\mathrm{e}}^{-\beta H(g)}$$

$$e^{-\beta (H_0+gV)}=e^{-\beta H_0}\left(\mathbb{I}-\beta g{\int }_0^{1}d\tau e^{\beta \tau H_0}V{e}^{-\tau \beta H_0}\right)+O(g^2).$$

$$Y_{H,\beta }:={\int }_0^{1}d\tau e^{\beta \tau H}Y{e}^{-\tau \beta H}.$$

$$|Y_{H,\beta }⟩$$

$$={\int }_0^{1}d\tau e^{\tau \beta (\mathbb{I}\otimes H^T-Z\otimes \mathbb{I})}|Y⟩$$

$$=\underset{ϵ\to 0}{\lim }\frac{e^{\beta (\mathbb{I}\otimes H^T-H\otimes \mathbb{I})}-\mathbb{I}}{\beta (\mathbb{I}\otimes H^T-H\otimes \mathbb{I})+iϵ}|Y⟩$$

Full text

I Introduction

Thermodynamics is the fundamental theory concerned with heat and temperature and their relation to energy and work. In phenomenological thermodynamics, an implicit assumption is that couplings between the working systems and their heat baths are so weak so that effects of the interaction can be neglected. As a consequence, the equilibrium states of the working systems are thermal states, in fact thermal states of local Hamiltonians. For small-scale systems governed by quantum mechanical laws, however, such a weak-coupling limit can be far from being reasonable, as the surface area of such systems is often not much smaller than their volume. An impressive body of literature in a related field, namely equilibration and thermalisation of closed quantum many-body systems Gogolin and Eisert (2016), strongly suggests that a system coupled strongly to a thermal bath should be described by the local reduced state of the global Gibbs state $\rho_{S}=\mathrm{Tr}_{B}(e^{-\beta H}/{\rm tr}(e^{-\beta H}))$ – and not by a Gibbs state of the local Hamiltonian itself Subasi et al. (2012); Müller et al. (2015); Brandao and Cramer (2015); Farrelly et al. .

In this work we take this basic but profound insight seriously when studying in detail quantum thermal machines strongly coupled to heat baths. First, we prove exact and general bounds on work extraction from a non-equilibrium system that can be brought in contact with a single heat bath. These results can be captured as universal corrections to the weak coupling limit – showing that strong coupling unavoidably leads to irreversibility and is hence detrimental for work extraction. Similar corrections are obtained for heat dissipation and the Carnot efficiency, hence providing strong-coupling corrections to the different formulations of the second law of thermodynamics. For thermal machines, we also show that strong interactions lead to power enhancements. Finally, we illustrate these considerations by means of the paradigmatic quantum Brownian motion Caldeira and Leggett (1983).

These results are put into context of the field of quantum thermodynamics Vinjanampathy and Anders (2016); Goold et al. (2016); Kosloff (2013), in particular, of recent efforts to describe the thermodynamics of quantum systems with strong interactions between system and bath Gelbwaser-Klimovsky and Aspuru-Guzik (2015); Uzdin et al. (2016); Katz and Kosloff (2016); Strasberg et al. (2016); Newman et al. (2017); Xu et al. (2016); Lee et al. (2012); Hänggi et al. (2008); Gelin and Thoss (2009); Campisi et al. (2009); Kim and Mahler (2010); Campisi et al. (2010); Hilt et al. (2011); Lee et al. (2012); Gallego et al. (2014); Lee et al. (2015); Esposito et al. (2015a, b); Bruch et al. ; Perarnau-Llobet et al. (2016); Kato and Tanimura (2016); Philbin and Anders (2016); Carrega et al. (2016) (see Refs. Seifert (2016); Jarzynski (2017) for classical systems). These include considerations on heat engines Gelbwaser-Klimovsky and Aspuru-Guzik (2015); Uzdin et al. (2016); Katz and Kosloff (2016); Strasberg et al. (2016); Newman et al. (2017); Xu et al. (2016), equilibrium and non-equilibrium thermodynamics Hänggi et al. (2008); Gelin and Thoss (2009); Campisi et al. (2009); Kim and Mahler (2010); Campisi et al. (2010); Hilt et al. (2011); Lee et al. (2012); Gallego et al. (2014); Lee et al. (2015); Esposito et al. (2015a, b); Bruch et al. ; Perarnau-Llobet et al. (2016); Kato and Tanimura (2016); Philbin and Anders (2016); Carrega et al. (2016) and, in a more abstract level, generic limitations on transformations between states using thermal resources Brandao et al. (2013, 2015); Cwiklinski et al. (2015); Lostaglio et al. (2015); Wilming et al. (2016). The key contribution of the present work, compared with earlier strong-coupling analyses of heat engines Gelbwaser-Klimovsky and Aspuru-Guzik (2015); Uzdin et al. (2016); Katz and Kosloff (2016); Strasberg et al. (2016); Newman et al. (2017); Xu et al. (2016); Lee et al. (2012), is to provide bounds on work and efficiency, without having to restrict to any particular model for the systems involved. Our bounds apply to thermodynamic scenarios in which the system equilibrates to the reduced of a global Gibbs state, and for which the coupling can be switched on and off. More precisely, our results are derived within a framework applicable to general situations; after all, also phenomenological thermodynamics is widely applicable by largely abstracting from the specifics of a given setting.

II Framework

We consider a system $S$, a heat bath $B$, with internal Hamiltonians $H_{S}$ and $H_{B}$, respectively. They can interact via a a possibly strong interaction $V$. Thermodynamic protocols then consist on transformations over $H_{S}$, and equilibration processes induced by $V$. Specifically, we consider protocols of $N$ steps, and denote by $\rho^{(i)}$ and $H^{(i)}$ the state and Hamiltonian of $SB$ in the $i$th step, consisting of three elementary operations

• (A)

Turning on/off interaction: With this, we model the process of bringing $S$ and $B$ into contact, so that the Hamiltonian takes the form $H^{(i)}=H_{S}^{(i)}+H_{B}+V$. Similarly, the interaction can be turned off at any step of the process. Treating such processes as quenches, the average work gain when placing/removing $V$ is

[TABLE]

• (B)

A quench on $S$: A fast transformation of $H_{S}$ is implemented, so that $H^{(i)}=H_{S}^{(i)}+H_{B}+V$ is changed to $H^{(i+1)}=H_{S}^{(i+1)}+H_{B}+V,$ whereas the state $\rho^{(i)}$ remains unchanged. The corresponding work gain reads

[TABLE]

which depends only on the state of $S$, since the interaction energy and the bath energy remain constant 111Quenches on $S$ can be similarly applied when the interaction $V$ is turned off. Note that in this case Eq. (2) is also valid..

• (C)

A thermalisation process. This operation models the closed free evolution of $SB$ when $V$ is present, i.e., under $H^{(i+1)}$. In this case, $S$ and $B$ exchange energy, while the total energy is preserved. Hence, this operation has no work cost. When they reach equilibrium, we assume that the state of $S$ is well described by

[TABLE]

where $\omega_{\beta}(H)=e^{-\beta H}/{\rm tr}(e^{-\beta H})$. Similarly, we assume that the boundary between $S$ and $B$, i.e., the support of $V$, can also be described by the reduced of a global thermal state.

Both assumptions are reasonable for locally interacting systems and are backed by a body of rigorous arguments Gogolin and Eisert (2016); Farrelly et al. 222Indeed, in the specific situation of the system $S$ and its bath $B$ being translationally invariant, Ref. Farrelly et al. shows that the reduced state of the infinite time average is close in $1$-norm to $\rho_{S}={\rm tr}_{B}(e^{-\beta H}/{\rm tr}(e^{-\beta H}))$, once one perturbs a Gibbs state by acting with a completely positive map on it supported on a constant number of sites only. Here, in contrast, we consider Hamiltonian quenches in $S$. However, in the high temperature limit Kliesch et al. (2014), the effect of such quenches can be approximated by completely positive maps acting on $S$ and a constant number of degrees of freedom of $B$. The closeness to the infinite time average is guaranteed for the overwhelming majority of times by invoking the equilibration results of Ref. Linden et al. (2009) in conjunction to the bounds on the effective dimension from Ref. Farrelly et al. .. When it is clear from the context we will use the notation $\omega^{(i)}:=\omega_{\beta}(H^{(i)})$ and $\omega_{S}^{(i)}:=\omega_{\beta}(H_{S}^{(i)})$. We also use the convention $\hbar=1$, $k_{B}=1$, and that when $SB$ decrease their global energy, then work is extracted and $W>0$.

A thermodynamic protocol then consists of an arbitrary sequence of operations of the type (A)-(C). The total expected work $W$ gained in the process is the sum of all the contributions of the form (1) and (2). Note that in this framework, the Hamiltonian terms $V$ and $H_{B}$ remain fixed throughout the protocol, reflecting the fact that an experimenter will in many realistic situations not have precise control over $B$ and the coupling between $S$ and $B$, at least not beyond the capability of turning it on and off. After every transformation of the form (C), $S$ is assumed to be brought to equilibrium after sufficiently long time. That is, possible finite-time effects are not included in this framework.

III Maximal work extraction for arbitrary coupling strengths.

We now study work extraction from an out-of equilibrium state of $S$. In order to avoid the possibility of extracting work from the energy stored in $V$, we consider that $S$ is initially isolated from $B$. The initial Hamiltonian is hence non-interacting, $H^{(0)}=H_{S}+H_{B}$, and the initial state is uncorrelated, $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B}).$ Given these initial conditions, the task is to optimize the extracted work over all cyclic Hamiltonian processes under the operations (A)-(C). Cyclicity here means that in a protocol of $N+1$ steps, we have $H^{(N+1)}=H^{(0)}$, where $N$ can be arbitrarily large.

It is instructive to first recall the optimal protocol in the weak-coupling regime Alicki et al. (2004); Strasberg and Esposito (2017); Anders and Giovannetti (2013); Gallego et al. (2014). It consists of four steps: (i) a quench from $H_{S}^{(0)}$ to $\tilde{H}_{S}$, where $\omega_{\beta}(\tilde{H}_{S})=\rho_{S}$, (ii) turning on $V$, (iii) an isothermal process from $\tilde{H}_{S}$ back to $H_{S}$, and (iv) turning off $V$. In our framework, isothermal processes correspond to a concatenation of infinitesimally small quenches followed by equilibration steps – we refer the reader to Refs. Anders and Giovannetti (2013); Gallego et al. (2014) for more details. The protocol (i)-(iv) has no dissipation, and is hence reversible, in the limit of an arbitrarily weak $V$. In the strong-coupling regime, where the energy of $V$ can no longer be neglected, we show in Appendix C that the optimal protocol also has the form (i)-(iv), but the initial and final Hamiltonians of the isothermal process need to be modified. Let $H_{S}^{(1)}$ and $H_{S}^{(N)}$ be the Hamiltonians of $S$ when $V$ is turned on and off, respectively. Then, the total work $W$ of the protocol can be expressed as

[TABLE]

where $W^{\rm(weak)}=F(\rho_{S},H_{S})-F(\omega_{\beta}(H_{S}),H_{S})$ is the maximal extractable work in the weak coupling regime, $F(\rho,H):=\mathrm{Tr}(\rho H)+T\mathrm{Tr}(\rho\ln\rho)$ is the (non-equilibrium) free energy, and

[TABLE]

with $H^{(1)/(N)}=H_{S}^{(1)/(N)}+H_{B}+V$. It is important to note that $F(\rho,H)-F(\omega_{\beta}(H),H)=T{\rm S}(\rho\|\omega_{\beta}(H))\geq 0$ with ${\rm S}(\rho\|\sigma)\coloneqq\mathrm{Tr}\left(\rho(\log\rho-\log\sigma)\right)$ the quantum relative entropy. It follows that always $\Delta F^{\mathrm{(irr)/(res)}}\geq 0$, and we can already conclude that strong coupling cannot be beneficial for work extraction as $W\leq W^{\mathrm{(weak)}}$. The correcting term $\Delta F^{\rm(irr)}$ can be interpreted as the energy dissipated when $S$ is put in contact with $B$, whereas $\Delta F^{\rm(res)}$ is the extractable work left on the final state.

The extracted work $W$ in (4) is maximised when $H_{S}^{(1)}$ and $H_{S}^{(N)}$ minimise the correcting terms $\Delta F^{\rm(irr)}$ and $\Delta F^{\rm(res)}$, respectively. Assuming that $\rho_{S}$ is a full rank state, we show in Appendix C that this happens for

[TABLE]

where we have defined $Y_{H,\beta}:=\int_{0}^{1}d\tau\hskip 2.84526pte^{\beta\tau H}Ye^{-\tau\beta H}$ for Hermitian operators $Y$ – an integral that can be solved analytically. Furthermore, $\Delta F^{\rm(irr)/(res)}$ have at least one minimum, so that Eqs. (7), (8) always provide the desired solution (see Appendix C). If more than one minimum exists, the solution corresponds to the absolute one.

It is important to stress that although a priori (7), (8), and in general the extracted work $W$, depend on the entire bath $B$, commonly its Hamiltonian is local and the correlations between its degrees of freedom decay rapidly with the distance. Therefore, only the degrees of freedom that are geometrically close to the boundary between $S$ and $B$ will contribute. This has the important consequence that we can solve  (7), (8) by considering a small buffer region in $B$ involving a few degrees of freedom only, while maintaining tight bounds for the error made in such a prescription (see Refs. Kliesch et al. (2014); Ferraro et al. (2012); Hernández-Santana et al. (2015)). This renders the solution practically and efficiently computable.

Altogether, our techniques provide a procedure to determine for any model the optimal protocol for work extraction in the strong coupling regime. Essentially, it consists of an isothermal process, where $S$ is put in contact with $B$ according to Eqs. (7) and (8). In what follows, we solve explicitly these equations at lowest order in the interaction strength.

IV Corrections at lowest order of work

Interestingly, in a perturbative treatment, the problem at hand can be essentially solved by computing covariances. We start by replacing $V$ by $gV$, where the dimensionless $g>0$ quantifies the interaction strength. Expanding Eqs. (7) and (8) in $g$, we get

[TABLE]

where we recall that $\tilde{H}_{S}$ is defined via $\rho_{S}=\omega_{\beta}(\tilde{H}_{S})$. Inserting (9) and (10) into (5) and (6) respectively, in Appendix C we obtain $\Delta F^{\rm(irr)/(res)}_{\rm min}:=\min_{H_{S}^{(1)/(N)}}\Delta F^{\rm(irr)/(res)}$ at lowest non-vanishing order in $g$

[TABLE]

Here, we have defined $\tilde{V}:=V-\mathrm{Tr}_{B}(V\omega_{\beta}(H_{B}))$, $\tilde{H}^{(0)}:=\tilde{H}_{S}+H_{B}$, and ${\rm cov}_{\omega_{\beta}(H)}(A,B)=\mathrm{Tr}(A_{H,\beta}B\omega_{\beta}(H))-\mathrm{Tr}(A\omega_{\beta}(H))\mathrm{Tr}(B\omega_{\beta}(H))$ is the generalized covariance Kliesch et al. (2014), also known as Kubo-Mori inner product in linear response theory Mori (1956); Kubo (1957); Mori (1965). Some important remarks are now in order,

• •

The first order correction to $W$ vanishes for any $H_{S}^{(1)}=\tilde{H}_{S}+O(g)$, $H_{S}^{(N)}=H_{S}+O(g)$. This follows from the penalty terms $\Delta F^{\rm(irr)/(res)}$ being differentiable functions of $g$ and having a minimum at $g=0$. The choice (9), (10) provides the minimum coefficient of $O(g^{2})$.

• •

The first order correction in (9), (10) exactly compensates for the term $g\mathrm{Tr}_{B}(\omega(H_{B})V)$ which often appears in open quantum systems as an effective action of $B$ on $S$ Weiss (1993); Breuer and Petruccione (2002).

• •

The generalised covariance ${\rm cov}_{\omega_{\beta}(H)}(A,B)$ captures the linear response of the thermal state under perturbations Mori (1956); Kubo (1957); Mori (1965).

V Heat and dissipation

Let us now turn to heat dissipation in an isothermal process in the strong coupling regime. For that, we do not consider a cyclic process, and we instead fix the intial and final Hamiltonian to be $H_{S}$ and $H_{S}^{(N)}$ respectively – specifically $H^{(N)}=H_{S}^{(N)}+V+H_{B}$ and $H^{(N+1)}=H_{S}^{(N)}+H_{B}$. We consider the same initial state as in the work-extracting protocol, i.e., $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B})$.

From the first law of thermodynamics, the total heat reads $Q=\Delta E_{S}+W$, with $\Delta E_{S}=\mathrm{Tr}(H_{S}^{(N)}\omega^{(N)})-\mathrm{Tr}(H_{S}\rho^{(0)})$. Since $\Delta E_{S}$ is fixed by $\rho_{S}$ and $H_{S}^{(N)}$, it becomes clear that the optimal protocol for maximising $W$ also minimises dissipation. Then from (4), we obtain (see Appendix D for details)

[TABLE]

where $\Delta{\rm S}={\rm S}(\mathrm{Tr}_{B}(\omega^{(N)}))-{\rm S}(\rho_{S})$ is the gain of entropy of $S$, $\Delta F^{\mathrm{(res)}}_{B}=T{\rm S}(\mathrm{Tr}_{S}(\omega^{(N)})\|\omega_{\beta}(H_{B}))$ is the increase of the free energy of $B$, and $I(\omega^{(N)};S:B)>0$ is the mutual information between $S$ and $B$. Note that in the strong coupling case, $Q<T\Delta{\rm S}$, even when the isothermal process is accomplished reversibly. Again, this is due to the penalising terms $\Delta F^{\mathrm{(res)}}_{B}$ and $\Delta F^{\mathrm{(irr)}}$, in addition to the correlations captured by the mutual information. Minimising dissipation corresponds to minimising the negative terms in (12). However, in this case $\Delta F^{\mathrm{(res)}}_{B}$ and $I(\omega^{(N)};S:B)$ are fixed through $H^{(N)}$. Hence, we only have freedom to minimise $\Delta F^{\mathrm{(irr)}}$, a problem that has been solved in (7).

Similar to the case of work, we can now expand the correcting terms over the interaction strength $g$. As before, the first order correction vanishes, so that the series expansion reads

[TABLE]

with $K_{q}>0$, and where we note that $\Delta{\rm S}$ depends on $\omega_{\beta}(H^{(N)})$ and hence indirectly also on $g$ 333If one wishes to expand $Q(g)$ with respect to $Q^{\rm(weak)}$, then the first order corrections do not vanish. This is discussed in more detail in Appendix D.. From (12), we note that a simple and useful lower bound for $K_{g}$ is given by $K_{q}\geq\Delta F^{\rm(irr)}_{\rm min}$ as given by (11). In other words, Eq. (11) also provides a strong coupling correction to dissipation, and to the Clausius formulation of the second law.

VI Heat engines

Given (4) and (12), it is straightforward to study the efficiency of a heat engine in the strong coupling regime. We consider engines made up of two baths at different temperatures which can sequentially interact strongly with $S$. More precisely, we extend our formalism to account for equilibrations, always in the form (3), with two baths $B_{c}$ or $B_{h}$, at two different (inverse) temperatures, $\beta_{c}$ and $\beta_{h}$. The task is then to maximise the efficiency of a cycle of the engine.

Not surprisingly, the optimal cycle turns out to have the same form as a Carnot engine, as we show in Section E of the Appendix. The Carnot cycle can be described with four steps

• •

an isothermal transformation in contact with $B_{h}$ from $H^{(A)}_{S}$ to $H^{(B)}_{S}$,

• •

a quench from $H^{(B)}_{S}$ to $H^{(C)}_{S}$,

• •

an isothermal transformation with $B_{h}$ from $H^{(C)}_{S}$ to $H^{(D)}_{S}$,

• •

a quench from $H^{(D)}_{S}$back to $H^{(C)}_{S}$.

In the weak coupling regime, the efficiency is maximised through the choice of Hamiltonians

[TABLE]

which guarantees no dissipation. Given our previous results, these conditions are naturally extended in the strong coupling regime to

[TABLE]

where $H^{(X)}:=H_{S}^{(X)}+V+H_{B}$. This provides a simple recipe for constructing minimal dissipation engines in the strong coupling regime. The corresponding (maximal) efficiency $\eta$, using $\eta=1-|Q_{c}|/|Q_{h}|$, (13) and expanding in $g$, reads

[TABLE]

where $\eta^{C}$ is the Carnot efficiency, $Q_{h/c}^{\rm(weak)}=T_{h/c}\Delta S^{\rm(weak)}$ and the entropy change in the weak-coupling regime is defined as $\Delta S^{\rm(weak)}=S(\omega_{\beta}(H_{S}^{(B)})-S(\omega_{\beta}(H_{S}^{(D)})$ and $K_{q}^{(h/c)}$ are coefficients obtained from (13) for $B_{h/c}$ (see details in Appendix E). By recalling the bound $K_{q}\geq\Delta F^{\rm(irr)}_{\rm min}$ (at order $O(g^{2})$), through (11) we obtain strong coupling corrections to the Carnot efficiency.

VII Macroscopic Limit

Let us briefly discuss the macroscopic limit, in which $S$ becomes large. The correcting terms to work and heat in (4) and (12) can be bounded by the interaction strength as

[TABLE]

where $\left\|V\right\|$ is the operator norm of $V$. The first bound is derived in Appendix C, whereas the second one follows from Ref. Wolf et al. (2008). Now we use these bounds to provide a simple argument to show that the limiting terms disappear when dealing with large systems. Let $S$ be a locally interacting system made up of $n^{3}$ particles. Let it be coupled also locally with $B$, so that the number of particles interacting with $B$ is $\alpha n^{2}$, $\alpha$ being a parameter that depends on the specific geometry of the boundary between $S$ and $B$. Let us write the interaction as

[TABLE]

where $h_{j}$ contains all interactions with the $j$-th particle of $S$ in the boundary. Now we have that,

[TABLE]

On the other hand, the extractable work from $S$, given by $\Delta F_{S}=F(\rho_{S},H_{S})-F(\omega_{\beta}(H_{S}),H_{S})$, will in general scale with the size of $n^{3}$, as the free energy is an extensive quantity. Hence, in the limit of large $n$ we have that both $\Delta F_{S}/\Delta F^{\rm(irr)}$ and $\Delta F_{S}/\Delta F^{\rm(res)}$ scale as $O(1/n)$, and hence disappear in the macroscopic limit. Thus, the above corrections become negligible in the limit of large systems. In other words, macroscopic phenomenological thermodynamics is insensitive to the strength of the underlying interactions; making these effects only relevant for small systems.

VIII Power

Although non-zero interactions between $S$ and $B$ tend to increase dissipation, they can help to enhance power of an engine, as they can decrease the time scale of thermalisation $\tau$. Following, we provide an upper bound for the power enhancement of a Carnot-like engine due to the interaction strength $g$. In order to do so, we need some considerations on how $\tau$ is related to $g$. A dimensional analysis argument rapidly suggests that $\tau\propto g^{-1}$, and through a more careful analysis implemented in Appendix F we obtain $\tau\geq\delta Q/gr$, with $\delta Q$ being the energy change of $B$ during the equilibration and $r\coloneqq\left\|[H_{B},V]\right\|$ the maximum rate in which $S$ and $B$ can exchange energy. These considerations allow us to obtain the following bound

[TABLE]

where $W$ is the work produced in a cycle, $\Delta t$ the cycle time, $\eta$ the efficiency of the machine, $r_{c/h}\coloneqq\left\|[H_{B_{c/h}},V_{c/h}]\right\|$ is the maximum rate at which the cold/hot bath with Hamiltonian $H_{B_{c/h}}$ can lose/gain energy and $V_{c/v}$ is the interaction that couples $S$ to the cold/hot bath (see for details Appendix F). Using (16) we can expand (20) in $g$, obtaining

[TABLE]

This result suggests that the power of an engine first increases with $g$, reaches a maximum, to then decrease for larger values of $g$ (see also discussion of the limit $g\rightarrow\infty$ in Appendix C). This behaviour is also observed in other treatments of heat engines in strong coupling for various figures of merit Gelbwaser-Klimovsky and Aspuru-Guzik (2015); Uzdin et al. (2016); Katz and Kosloff (2016); Strasberg et al. (2016); Newman et al. (2017); Xu et al. (2016).

Let us point out that while our bound for $P$ is valid for arbitrary systems, it is expected to be very crude in general, as it depends on the norm of the interaction, $\left\|[H_{B_{c/h}},V_{c/h}]\right\|$. Further progress could be made by either considering particular models or by obtaining better bounds on the time scales of equilibration of generic systems, a notoriously hard and diverse problem Gogolin and Eisert (2016). Finally, note that relations between power and efficiency of Carnot engines have also been obtained in Refs. Allahverdyan et al. (2013); Shiraishi et al. (2016); Shiraishi and Tajima , yielding complementary results.

IX Thermodynamic protocols within the Caldeira-Leggett model

We now illustrate our findings by applying them to the model of quantum Brownian motion, captured in the standard Ullersma or Caldeira-Leggett (CL) model Ullersma (1966); Caldeira and Leggett (1983). In this model,

[TABLE]

where $S$ is an harmonic oscillator, $H_{S}=(m\omega^{2}x^{2}+{p^{2}}/{m})/2,$ $B$ a bosonic bath , $H_{B}=\sum_{k}(m_{k}\omega_{k}^{2}x_{k}^{2}+p_{k}^{2}/m_{k})/2$, which is linearly coupled to $S$ through $gV$ with $V=x\sum_{k}g_{k}x_{k}$ and where $g$ quantifies the strength of the interaction, and finally $H_{L}=x^{2}g^{2}\sum_{k}{g_{k}^{2}}/({m_{k}\omega_{k}^{2}})$ is a renormalization term. Indeed, for Ohmic spectral densities, $g$ at the same time quantifies the deviation from Markovian dynamics Groeblacher et al. (2015). The couplings $g_{k}$ are determined by the spectral density $J(w)=({\pi}/{2})\sum_{k}({g_{k}^{2}}/{\omega_{k}^{2}})\delta(\omega-\omega_{k}),$ which in the continuum limit is often taken to be well-approximated by an Ohmic function, $J(\omega)\propto\omega$, for low frequencies until some cut-off $\Omega>0$. This model plays a crucial role in the study of open quantum systems Weiss (1993) and finds numerous applications in thermodynamics Hänggi et al. (2008); Eisert and Plenio (2002); Allahverdyan and Nieuwenhuizen (2000); Nieuwenhuizen and Allahverdyan (2002); Kim and Mahler (2010); Hilt et al. (2011); Carrega et al. (2015); Philbin and Anders (2016); Strasberg et al. (2016); Newman et al. (2017).

The equilibrium state of $S$ in the CL-model is given by $\mathrm{Tr}_{B}(\omega_{\beta}(H))$, irrespectively of $g$, thus satisfying (3) Ambegaokar (2007); Subasi et al. (2012) 444Numerical evidence also shows that the energy contained in the interaction $V$ also equilibrates to the reduced of a global thermal state.. Furthermore, the total Hamiltonian is quadratic, and hence can be solved exactly with matrices of order $O(n^{2})$, where $n$ is the number of oscillators in the bath (see, e.g., Refs. Rivas et al. (2010); Eisert et al. (2010)). This allows us to numerically simulate thermodynamic protocols exactly for arbitrary strong coupling and large (but finite) baths. Details on the discretisation of the CL model and its simulation are provided in Appendix G. There, the equilibration time is also discussed, finding $\tau\propto 1/g^{2}$, for $g\leq 1$, using techniques from Refs. Oliveira et al. (2017); Wilming et al. (2017), a result which agrees with standard perturbative approaches in the macroscopic limit Weiss (1993).

We now illustrate our results for work extraction using the CL model. Crucially, the first order corrections in (9) and (10) vanish, as the thermal state of $H_{B}$ is symmetric under $x_{k}\leftrightarrow-x_{k}$. This implies that the optimal protocol in the weak coupling regime is in fact also optimal for small but non-zero $g$. This is perfectly illustrated in Fig. 1, where we plot the work extracted using the weak coupling protocol and the optimal one, which is obtained by numerically minimising $\Delta F^{(\rm res)/(irr)}$. It is clear that differences between the two start appearing only at higher orders than $O(g^{2})$. Note also that Fig. 1 shows an excellent agreement between the exact unitary dynamics and our framework, in which (3) is assumed, even when many quenches are performed.

Now we turn to the question of power. Here we keep the number of quenches $N$ fixed and vary the coupling strength $g$. Since we deal with isothermal processes, for which $N\rightarrow\infty$, we take $N$ large but finite. As a result of the equilibration time $\tau\propto 1/g^{2}$ for $g\leq 1$, the power $P(g)=W(g)/\tau(g)$, scales as $P(g)\propto g^{2}W^{\text{(weak)}}-O(g^{3})$. This relation is shown in Fig. 1, where we see that $P(g)$ first increases as $g^{2}$, reaches a maximum, and then decays to zero for large $g$.

X Conclusion

Bringing together arguments from quantum thermodynamics and the theory of equilibration in closed many-body systems, we have derived general strong coupling corrections to the second law of thermodynamics. These corrections are applicable to any model of interest, and have been obtained by designing optimal thermodynamic protocols in the strong coupling regime. The corrections become relevant if the working body is a small system, and vanish in first order with the interaction strength. An upper bound on the power enhancement due to the interaction strength has also been derived. A particularly relevant open problem is to extend these considerations to scenarios where the system is simultaneously strongly coupled to more than one thermal bath. In this case reaction coordinate mappings Strasberg et al. (2016); Newman et al. (2017) appear as a promising technique to extend results in the weak coupling limit Brunner et al. (2012). It is the hope that this work further stimulates the emerging field of strong-coupling quantum thermodynamics, aiming at identifying the potential and burden coming along with such interactions.

XI Acknowledgements.

We thank Ignacio Cirac, Andras Molnar, and Tao Shi for insightful discussions. We also thank P. Strasberg, G. Schaller and D. Gelbwaser-Klimovsky for useful comments on the manuscript. HW acknowledges funding from the ERC (TAQ) and the Studienstiftung des Deutschen Volkes, MPL from the Alexander von Humboldt foundation, RG from the DFG (GA 2184/2-1). JE thanks for support by the DFG (CRC 183, EI 519/7-1), the ERC (TAQ), the EU (AQuS), and the Templeton Foundation. AR is supported by the Spanish MINECO (QIBEQI FIS2016-80773-P, FISICATEAMO FIS2016-79508-P and Severo Ochoa Grant SEV-2015-0522), Fundació Privada Cellex, Generalitat de Catalunya (Grant SGR 874, 875, and CERCA Programme), the EC [FP7-ICT-2011-7, Grant No. 288263; OSYRIS (ERC-2013-AdG Grant 339106)] and the Beatriu de Pinós fellowship (BP-DGR 2013). All authors are grateful for support from the EU COST Action MP1209 on “Thermodynamics in the quantum regime”.

Appendix

This Appendix contains eight sections. In Section A and B, notation and the basic mathematical tools are introduced. Section C contains all results concerning maximal work extraction in the strong coupling regime. Section D deals with heat dissipation, Section E with heat engines, and Section F with power. Finally, in Section G the simulation and equilibration times of the Caldeira-Legget model are discussed.

Appendix A Notation

For the sake of clarity, here we recall the notation introduced in the main text:

• •

$T,\beta$: Temperature and inverse temperature, respectively.

• •

$H_{S}$, $H_{B}$, and $V$: Hamiltonian of $S$, $B$ and the interaction, respectively.

• •

$H^{(i)}$: Hamiltonian of $SB$ in the $i$-th step of the protocol. It may contain $V$, so that $H^{(i)}=H_{S}^{(i)}+V+H_{B}$.

• •

$\omega_{\beta}(H)$: Thermal state,

[TABLE]

We will often use the shorthand notation $\omega^{(i)}:=\omega_{\beta}(H^{(i)})$.

• •

$\mathcal{Z}$: Partition function, $\mathcal{Z}=\mathrm{Tr}(e^{-\beta H})$.

• •

$\rho_{S}^{(i)},\rho^{(i)}_{B}$: Reduced states of $S$ and $B$ at the $i$-th step of the protocol, $\rho_{S/B}^{(i)}=\mathrm{Tr}_{B/S}\omega^{(i)}$.

• •

$\rho_{S},\tilde{H}_{S},\tilde{H}^{(0)}$: $\rho_{S}$ initial state of $S$, $\tilde{H}_{S}$ is defined through $\rho_{S}=\omega_{\beta}(\tilde{H}_{S})$, and finally $\tilde{H}^{(0)}=\tilde{H}_{S}+H_{B}$.

• •

$F(\rho,H)$: Non-equilibrium free energy, $F(\rho,H)=\mathrm{Tr}(H\rho)-TS(\rho)$, with $S(\rho)=-\mathrm{Tr}(\rho\ln\rho)$.

• •

The operator norm $\left\|X\right\|$ of an operator $X$ is given by

[TABLE]

for normalized state vectors $\left.\left|{\psi}\right.\right\rangle$. In the finite-dimensional spaces considered here, it is simply the largest spectral value of $X$. If $X$ is Hermitian, then we have $\left\|X\right\|=\max_{j}|\lambda_{j}|$, where $\lambda_{j}$ are the eigenvalues of $X$.

Appendix B Preliminaries

Let us also introduce a few technical tools.

B.1 Differentiation of matrix functions

Let $A,B$ be two matrices. Then we can differentiate the trace of some matrix function $f(A)$ as

[TABLE]

where $f^{\prime}$ is the derivative of $f$. For example, given $H=H_{0}+gV$, and with $\mathcal{Z}=\mathrm{Tr}(e^{-\beta H})$, $\mathcal{Z}_{0}=\mathrm{Tr}(e^{-\beta H_{0}})$ we obtain

[TABLE]

B.2 Derivative and expansion of exponential operators

We will also use the derivative of exponentials of operators

[TABLE]

From (26), we can easily obtain the first order of the Dyson series (in imaginary time)

[TABLE]

In what follows we will use the short hand notation for the first term of the expansion,

[TABLE]

This integral can in fact be carried out. For that, we transform the operators of the Hilbert space into vectors, $\left.\left|{Y}\right.\right\rangle=\sum_{i,j}Y_{i,j}\left.\left|{i,j}\right.\right\rangle$ where $Y=\sum_{i,j}Y_{i,j}\left.\left|{i}\right.\right\rangle\left.\left\langle{j}\right.\right|$, obtaining

[TABLE]

where $\epsilon>0$ is introduced to ensure that the result is well defined when the denominator has a zero eigenvalue. Of course, the solution is a vector that should be transformed back to the original operator basis.

B.3 Perturbations of thermal states and generalised covariance

In the following we will make use of the generalised covariance Kliesch et al. (2014), also known as Kubo-Mori inner product Mori (1956); Kubo (1957); Mori (1965) or Bogoliubov inner product Petz and Toth (1993). For any state $\rho$ and two observables $A,B$ it is defined as

[TABLE]

For thermal states $\omega_{\beta}(H)$ we have

[TABLE]

From a physical point of view, the generalised covariance measures how thermal perturbation values respond to a change of the Hamiltonian in a thermal state Mori (1956); Kubo (1957); Mori (1965). Let $H_{t}$ be a smooth family of Hamiltonians. Then it follows from (26) that

[TABLE]

where $H^{\prime}_{0}=\mathrm{d}H_{t}/\mathrm{d}t|_{t=0}$. For finite perturbations of the Hamiltonian $H_{0}$ one can integrate the above equation (see Kliesch et al. (2014)).

B.4 Symmetry of relative entropy

In general, the quantum relative entropy is an asymmetric function, i.e., $S(\rho\|\sigma)\neq S(\sigma\|\rho)$. Nevertheless, it is known that it is symmetric up to second order in the difference between $\rho$ and $\sigma$ (see, e.g., Ref.  Amari and Nagaoka (2007)). The following Lemma, which might be of independent interest, shows this result when $\rho$ and $\sigma$ are thermal states, and their distance is modified by changing their Hamiltonians. This result will be key for the derivation of the exact form of the corrections to work and efficiency for small but finite coupling strengths.

Lemma 1** (Perturbative symmetry of relative entropy).**

Let $H_{t}$ be a one-parameter family of operators. Then for small $t$, we have

[TABLE]

Proof.

In the following we will use $f^{\prime}_{t}$ as shorthand for the derivative of the function $f_{t}$ with respect of $t$. Let us also define the partition function $\mathcal{Z}_{t}\coloneqq\mathrm{Tr}(\mathrm{e}^{-\beta H_{t}})$. Using (26), we obtain $\log(\mathcal{Z}_{t})^{\prime}=-\beta\mathrm{Tr}(\omega_{t}H^{\prime}_{t})$. We can then calculate the derivative of $S(\omega_{0}\|\omega_{t})$ as

[TABLE]

To compute the derivative of $S(\omega_{t}\|\omega_{0})$, let us first compute the derivative of the entropy of $\omega_{t}$, to get

[TABLE]

We then obtain for the derivative of the relative entropy

[TABLE]

where (36) follows from $\mathrm{Tr}(\omega^{\prime}_{t})=0$. Then, we have for the first derivative of $\Delta_{t}$

[TABLE]

From this expression we can easily compute the second derivative as

[TABLE]

In particular, we get $\Delta_{0}^{\prime}=0$ and $\Delta_{0}^{\prime\prime}=0$, which proves the claim. ∎

Appendix C Maximizing work extraction in the strong coupling regime

In this section, we optimize work extraction protocols over all cyclic protocols that can be constructed with the operations (A)-(C) described in the main text, and given the initial non-interacting Hamiltonian $H^{(0)}=H_{S}+H_{B}$ and the initial state $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B})$. Here cyclicity is understood in terms of the Hamiltonian, so that at the end of the protocol consisting of $N+1$ steps we have that $H^{(N+1)}=H^{(0)}$. Let us stress that these $N+1$ steps can in principle consist of any of the operations (A), (B) and (C) of the main text, and $N$ can be arbitrarily large.

We now focus on protocols where the interaction is turned on and off only once in the protocol (i.e., operation (A) is only implemented twice). Later we will show that this is in fact optimal. Any protocol can then always be described as,

A series of quenches are applied to the local Hamiltonian of $S$, $H_{S}\rightarrow...\rightarrow H_{S}^{(1)}$, and the interaction between $S$ and $B$ is turned on. The total Hamiltonian becomes $H^{(1)}=H_{S}^{(1)}+V+H_{B}$. The expected work gain of this process reads

[TABLE] 2. 2.

The Hamiltonian of $S$ is modified while $S$ is in contact with $B$, amounting to a sequence of quenches $H_{S}^{(1)}\mapsto H_{S}^{(2)}\mapsto\dots\mapsto H_{S}^{(N)}$, each followed by an equilibration to the corresponding thermal state, until $H^{(N)}=H_{S}^{(N)}+V+H_{B}$ is reached. Note that we assume that $H_{S}^{(N)}$ is independent of $N$, so that increasing $N$ means doing the same Hamiltonian transformation in more steps (i.e., slower). During this process the state of $S$ reads $\rho_{S}^{(j)}=\mathrm{Tr}_{B}(\omega^{(j)})$, $j=1,\dots,N$. The total work gain is then

[TABLE] 3. 3.

The interaction between $S$ and $B$ is turned off, and the Hamiltonian of $S$ is brought back to the initial form by an arbitrary series of quenches, $H^{(N)}_{S}\rightarrow...\rightarrow H_{S}$. The expected work gain is simply

[TABLE]

These three steps conclude a cyclic Hamiltonian process. Two remarks are now in order. Firstly, the last term of (40) is positive and tends to zero in the limit $N\rightarrow\infty$, i.e., for isothermal processes. That implies

[TABLE]

where $W^{\mathrm{(isoth)}}_{\rm sc}$ stands for the work gain of an isothermal transformation in the strong coupling regime Gallego et al. (2014). Hence, in the optimal protocol contacts between $S$ and $B$ must be in form of isothermal transformations.

The second remark concerns points 1 and 3. If the interaction is weak, the energy of turning on and off the interaction can be neglected. Exactly in this case, one can obtain the usual expression of optimal work given by the free energy. This is indeed obtained with the choice of $H_{S}^{(N)}=H_{S}$ and $H_{S}^{(1)}$ such that $\omega_{\beta}(H_{S}^{(1)})=\rho_{S}$. Then one obtains

[TABLE]

with $\omega^{(0)}:=\omega_{\beta}(H^{(0)})$. However, if the interaction is non-negligible, we now show that steps 1 and 3 become a source of irreversibility. As a consequence, the optimal protocol must consist of only one thermal contact between $S$ and $B$ – hence justifying the form of the protocol considered – further contacts can only decrease the extractable work.

Now we perform some simple algebra to express the total work of the optimal protocol in a convenient form. Adding up the three work contributions, and adding and subtracting $F(\omega^{(0)},H^{(0)})$, we obtain

[TABLE]

Using that $F(\rho^{(0)},H^{(1)})+\mathrm{Tr}((H_{S}-H_{S}^{(1)}-V)\rho^{(0)})=F(\rho^{(0)},H^{(0)})$ and $F(\omega^{(N)},H^{(N)})-F(\omega^{(N)},H^{(0)})=\mathrm{Tr}((V+H_{S}^{(N)}-H_{S})\omega^{(N)})$ we further obtain

[TABLE]

where we have defined

[TABLE]

By noting that $F(\rho,H)-F(\omega_{\beta}(H),H)=T{\rm S}(\rho\|\omega_{\beta}(H))$, we already obtain that always $W\leq W^{\rm(weak)}$, so that interactions are detrimental for work extraction. However, at the moment, the terms (45) depend on $H_{S}^{(1)}$ and $H_{S}^{(N)}$, i.e., on the particular points where $S$ is connected and disconnected from $B$. In order to maximise $W$, we hence need to minimize $\Delta F^{\rm(res)}$ and $\Delta F^{\rm(irr)}$ over $H_{S}^{(1)}$ and $H_{S}^{(N)}$, i.e.,

[TABLE]

with

[TABLE]

We now proceed to solve these minimizations, which can be carried out independently. For that, we can use that in a (local) minimum of a function, its derivative must vanish. Since here we minimise over a matrix, the relation (25) becomes useful.

C.1 Minimization of $\Delta F^{\rm(irr)}$

Let $X_{S}$ be the choice of $H_{S}^{(1)}$ yielding $\Delta F^{\rm(irr)}_{\rm min}$, and define

[TABLE]

Then, for any $Y_{S}$, it must hold that

[TABLE]

That is, if we perturb the solution $X_{S}$ by $tY_{S}$, then the derivative w.r.t $t$ must be zero for all $Y_{S}$. In other words, we are standing in a minimum. Conversely, we can use (49) to find $X_{S}$. Indeed, computing (49) using (25), we obtain

[TABLE]

Since this holds for any $Y_{S}$, and setting $X:=X(0)$, it follows that

[TABLE]

This matrix equation (of the size of $S$) implicitly provides $X_{S}$, and hence also $\Delta F^{\rm(irr)}_{\rm min}$.

We now show that there is at least one minimum of $\Delta F^{\rm(irr)}$, so that (51) always provides the desired solution, provided that $\rho_{S}$ is a full rank state. For that, we first show that for any full rank state $\rho_{S}$ and non-trivial Hamiltonian $H_{S}^{(1)}$, the term $\Delta F^{\rm(irr)}$ diverges with the operator norm $H_{S}^{(1)}$. More explicitly, let us parametrize the Hamiltonians as $H_{S}^{(1)}=\lambda\hat{H}_{S}^{(1)}$ with $\hat{H}_{S}^{(1)}$ controlling its direction, i.e. $\|\hat{H}_{S}^{(1)}\|=1$, and $\lambda$ its operator norm. Let $P_{S}^{(1)}$ denote the ground-state space of $\hat{H}_{S}^{(1)}$. Then the state

[TABLE]

is supported only with the subspace $P_{S}^{(1)}\otimes\mathcal{H}_{B}$, where $\mathcal{H}_{B}$ denotes the Hilbert-space of the bath. But then we have

[TABLE]

since the relative entropy diverges if the support of the first argument is not contained in that of the second. Now, the function $\Delta F^{\mathrm{(irr)}}$ is (i) positive, (ii) continuous and differentiable, (iii) and tends to $+\infty$ for large Hamiltonians $H^{(1)}_{S}$. It follows from (i)-(iii) that $\Delta F^{\rm(irr)}$ has at least one minimum, which can be obtained through (51). In case that the solution of (51) is not unique and there is more than a minimum, the global one is chosen.

C.2 Minimization of $\Delta F^{\rm(res)}$

Let us now proceed similarly for $\Delta F^{\rm(res)}_{\rm min}$, a calculation which turns out to be bit more involved. Similarly as in the last section, assume that $Z_{S}$ is the solution to $\min_{H_{S}^{(N)}}\Delta F^{\rm(res)}$, and $tY_{S}$ is some perturbation over this solution. Then, we define

[TABLE]

and $\mathcal{Z}(t)=\mathrm{Tr}\left(e^{-\beta Z(t)}\right)$ and $\omega_{\beta}(t)\coloneqq\omega_{\beta}(Z(t))$. The condition of minimum implies

[TABLE]

From the calculation in the proof of Lemma 1, we find that

[TABLE]

We can now use the perturbation-formula (32) and obtain

[TABLE]

Setting $Z\coloneqq Z(0)$ and writing out the definition of the generalised covariance, we then obtain

[TABLE]

Since this condition has to hold for all $Y_{S}$, it is, per definition of the partial trace, equivalent to

[TABLE]

This is the desired implicit solution for $Z_{S}$.

To finish this section, we show that the residual free energy $\Delta F^{({\rm res})}$ has at least one local minimum, and hence (58) always provides a solution. Unlike $\Delta F^{{\rm(irr)}}$, $\Delta F^{({\rm res})}$ does not diverge to $+\infty$ with the operator norm of $H_{S}^{(N)}$, but tends to a constant. As above, the Hamiltonian is parametrised as $H_{S}^{(N)}=\lambda\hat{H}_{S}^{(N)}$, where $\lambda$ controls its norm and $\hat{H}_{S}^{(N)}$ controls its direction, i. e. $\|\hat{H}_{S}^{(N)}\|=1$. As in the last section, let $P_{S}^{(N)}$ denote the ground-state subspace of $\hat{H}_{S}^{(N)}$. In the limit $\lambda\to\infty$, the thermal state $\omega_{\beta}^{(N)}(\lambda)=\omega_{\beta}(\lambda\hat{H}^{(N)}_{S}+V+H_{B})$ tends to the state

[TABLE]

which is again supported within the subspace $P_{S}^{(N)}\otimes\mathcal{H}_{B}$. Since $\Delta F^{\rm(res)}=F_{\beta}(\omega_{\beta}^{(N)}(\lambda),H^{(0)})-F_{\beta}(\omega_{\beta}(H^{(0)}),H^{(0)})$, we then have

[TABLE]

In order to see that the minimum of the residual free energy takes place for a finite $\lambda$ we show that the asymptotic value above is approached from below, and equivalently that the derivative for large $\lambda$ tends to 0 from the positive side. The first derivative of the residual free energy reads

[TABLE]

where $\omega^{(N)}(\lambda)\coloneqq\omega_{\beta}(\lambda\hat{H}^{(N)}_{S}+V+H_{B})$. We can use (32) to obtain

[TABLE]

In the limit $\lambda\to\infty$, we have

[TABLE]

since $\sigma_{SB}^{(N)}$ is supported within the ground-state space of $H^{(N)}_{S}$. Furthermore, in the limit $\lambda\rightarrow\infty$, we have $\exp(-\tau\lambda\hat{H}^{(N)}_{S}+V+H_{B})\hat{H}^{(N)}_{S}=0$ for any $\tau>0$. This then also implies

[TABLE]

We can then expand an expression of the form

[TABLE]

in $y=\exp(-\beta\lambda)$ around $y=0$. The above considerations show that the first order vanishes, hence for very large $\lambda$ we have

[TABLE]

Since

[TABLE]

for all $\lambda$, we have $f(\hat{H}^{(N)}_{S})\geq 0$. We can now use these results in Eq. 61 and find for large $\lambda$

[TABLE]

up to terms $O(\mathrm{e}^{-2\beta\lambda})$ and higher. We hence see that the derivative of $\Delta F^{\rm(res)}$ becomes positive for large enough $\lambda$ and $\Delta F^{\rm(res)}$ approaches its limiting value from below for large $\lambda$. Since it is a positive, smooth function, this implies the existence of at least one minimum and (58) always provides a solution to the problem of maximising work extraction. Again, if several solutions exist, the global optimum must be chosen.

In sum, so far we have reduced the problem of finding the optimal protocol for work extraction to the solution of two matrix equations, (51) and (58). We now provide the solution of these equations, and hence an expression for $W_{\rm max}$ in (46), at lowest non-vanishing order in the interaction strength. For that, let us replace $V$ by $gV$, and expand relevant quantities over $g$.

C.3 Expansion in orders of $g$

Recall that $X_{S}$ and $Z_{S}$ are the local Hamiltonians on $S$ minimizing $\Delta F^{(\rm irr)}$ and $\Delta F^{(\rm res)}$ respectively, and also

[TABLE]

The Hamiltonians $X_{S}$ and $Z_{S}$ are functions of $g$ since they depend on the interaction given by $gV$. This can be expressed by defining the functions $X_{S}(g)$ and $Z_{S}(g)$, which we can formally expand in powers of $g$ as

[TABLE]

where $X^{\prime}_{S}$ is the derivative with respect to $g$ of $X_{S}(g)$ and equivalently for $Z_{S}(g)$. Using (68) and (69) we can write an expansion in $g$ for the Hamiltonians

[TABLE]

Consider also the expansions of the corresponding thermal states,

[TABLE]

Note that, from (51) and (58), it is clear that $Z_{0}=H^{(0)}$, $X_{0}=\tilde{H}^{(0)}$, and hence $\omega_{0}^{(Z)}=\omega_{\beta}(H^{(0)})$, $\omega_{0}^{(X)}=\omega_{\beta}(\tilde{H}^{(0)})$. Moreover, we have that,

[TABLE]

Let us now expand these expressions. The calculation follows from the proof of Lemma 1 and yields

[TABLE]

and

[TABLE]

Proceeding analogously for $\Delta F^{\rm(res)}$, we obtain

[TABLE]

Importantly, from these calculations we learn that (i) the first non-vanishing terms are of order $O(g^{2})$ and only depend on first order corrections of $Z$ and $X$, and that (ii) the two correcting terms have the same behavior at second order in $g$, which is a consequence of Lemma 1.

We now compute explicitly the expansions (72) and (73). We start by computing the expansion of $X$ as determined by the solution (51), which, expanded using (72) and (73), can be simply written as

[TABLE]

From (26) and (28), and using the cyclic property of the trace, we obtain

[TABLE]

Taking the partial trace over $B$, one can express (79) as

[TABLE]

Recalling that $X_{1}=X^{\prime}_{S}(0)+V$ and defining $V^{\prime}:=-\mathrm{Tr}_{B}(V\omega_{\beta}(H_{B}))$ one can re-express the equation above as

[TABLE]

where $L$ is a linear super-operator defined as $L(\cdot):=\mathrm{Tr}(\omega_{\beta}(\tilde{H}_{S})\>\>\cdot\>\>)\mathbf{1}-(\>\cdot\>)_{\tilde{H}_{S}}$. The solution is then unique, given by $X_{S}^{\prime}(0)=V^{\prime}=-\mathrm{Tr}_{B}(V\omega_{\beta}(H_{B}))$ up to addition of any operator $M$ fulfilling $L(M)=0$. We will later see that adding any operator $M$ does not alter the correction to work, hence all the possible solutions perform equally good regarding work extraction. Altogether we find that

[TABLE]

for any $M$ fulfilling $L(M)=0$. In practice, since $M$ does not contribute to the extracted work, we can just set $M=0$.

We can proceed analogously for $\Delta F^{\rm(res)}$, expanding (58) at first order in $g$. One obtains an equivalent condition as for the case of $X_{S}$, that is, one obtains

[TABLE]

which can be again re-expressed as

[TABLE]

where $R(\cdot):=\mathrm{Tr}(\omega_{\beta}(H_{S}^{(0)})\>\>\cdot\>\>)\mathbf{1}-(\>\cdot\>)_{H^{(0)}_{S}}$. This in turn yields the solution

[TABLE]

for any $M^{\prime}$ fulfilling $R(M^{\prime})=0$.

We now use the solutions for the optimal Hamiltonians $X_{1}$ and $Z_{1}$ at first order, given by (82) and (85), to obtain the lowest order corrections to $\Delta F^{\rm(irr)}_{\rm min}$. Introducing the dressed interaction

[TABLE]

we obtain families of solutions given by

[TABLE]

with any $M$ and $M^{\prime}$ fulfilling $L(M)=0$ and $R(M^{\prime})=0$. Using the perturbation formula (32), we then have

[TABLE]

Eq. (88) follows since $L(M)=0$ implies that $\mathrm{Tr}(\omega_{\beta}(\tilde{H}^{(0)})M)\mathbf{1}=M_{\tilde{H}^{(0)}}$ and from the the definition of the generalised covariance (B.3). An equivalent argument can be used to show (89). Finally, using (77) and (78), we reach our main result

[TABLE]

Let us now summarize our results. The problem we have addressed is to maximise work extraction in a cyclic Hamiltonian process under the set of operations $(A)-(C)$ of the main text, and for an initial non-interacting Hamiltonian $H^{(0)}=H_{S}+H_{B}$ and the initial state $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B})$. Firstly, by expressing the maximal work as the difference between the maximal extractable work in the weak coupling regime and two positive terms (46), we have shown that strong interactions can only be detrimental for work extraction. More importantly, we have reduced the initial optimization over all protocols to two minimizations over local Hamiltonians as given by (C). These two minimizations have been solved for arbitrary Hamiltonians, giving rise to the two implicit solutions (51) and (58). Finally, by expanding the interaction $V$ over its strength $g$, we have obtained explicit solutions for the minimizations, c.f. Eqs. (82) and (85), which have allowed us to compute the correcting terms due to strong coupling up to order $O(g^{3})$, see (91).

C.4 Upper bounds and the macroscopic limit

In this section, we derive simple upper bounds to the strong-coupling correction terms, $\Delta F^{\rm(irr)}_{\rm min}$ and $\Delta F^{\rm(res)}_{\rm min}$, which in particular imply that they are bounded by $2\left\|V\right\|$ In this section we make no assumptions on the strength of the interaction. We begin with the residual free energy $\Delta F^{\rm(res)}$, obtaining

[TABLE]

In this expression, we have used again the relation $F(\rho,H)-F(\omega_{\beta}(H),H)=T{\rm S}(\rho\|\omega_{\beta}(H))$. Since the relative entropy is positive, we obtain

[TABLE]

which is the desired result. One can proceed in a similar manner for $\Delta F^{\rm(irr)}_{\rm min}$. First, one realizes that

[TABLE]

where recall that $\tilde{H}_{S}$ satisfies $\rho_{S}=\omega_{\beta}(\tilde{H}_{S})$. Then,

[TABLE]

This completes the derivation of the upper bounds for $\Delta F^{\rm(irr)}_{\rm min}$ and $\Delta F^{\rm(res)}_{\rm min}$.

C.5 Limit of arbitrarily strong interactions

Let us here briefly discuss the limit $g\rightarrow\infty$ of arbitrarily strong interactions. To do this, denote by $P_{V}$ the (projector onto the) ground-state subspace of the interaction $V$ and denote by $H^{(0)}|_{V}$ the operator $H^{(0)}=H^{(0)}_{S}+H_{B}$ restricted to the ground-state subspace of $V$ and similarly for $\tilde{H}^{(0)}$. Note that if $V$ is a local interaction only acting on some finite region, then its ground-state subspace is highly degenerate. In the case of very strong coupling we then obtain

[TABLE]

where the direct sum is over $P_{V}$. Thus, an arbitrarily strong interaction effectively restricts the Hilbert-space to a sub-space. This has important consequences for the correction terms. Consider the irreversible free energy $\Delta F^{\mathrm{(irr)}}$, defined as

[TABLE]

In the limit $g\rightarrow\infty$, the second argument is only supported within the ground-state subspace of $V$, while the support of the first argument is not contained in this subspace. But as is well known, the relative entropy $S(\rho\|\sigma)$ diverges whenever the support of $\rho$ is not contained in the support of $\sigma$. We thus find that

[TABLE]

Thus in the limit of arbitrarily strong interactions, the optimal protocol is to never couple the system to the bath. In this case, work can only be extracted if the initial state $\rho(0)$ is not passive Pusz and Woronowicz (1978); Lenard (1978). In particular, no work can be extracted if $\rho(0)=\omega_{\beta(0)}(H^{(0)}_{S})$ for some $\beta(0)>0$ as is the case in cyclically working thermal machines.

Appendix D Heat dissipation

In thermodynamics, it is often the case that the optimal protocols for one task turns out to be also optimal for others. Here we apply this logic to show that the previous results can be readily applied in order to minimise heat dissipation and maximise efficiency of heat engines strongly coupled to baths.

In order to study heat dissipation, we focus on (non-cyclic) processes where an initial Hamiltonian $H^{(0)}_{S}$ is modified to $H^{(N)}_{S}$ (with $H^{(N)}_{S}$ fixed and independent of $N$) by putting it in contact with a strongly interacting bath. More precisely, we consider an initially non-interacting Hamiltonian, $H^{(0)}=H^{(0)}_{S}+H_{B}$, an initial state $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B})$, and the following family of protocols,

A series of quenches are applied to the local Hamiltonian of $S$, $H_{S}^{(0)}\rightarrow...\rightarrow H_{S}^{(1)}$. After that, the interaction between $S$ and $B$ is turned on. 2. 2.

The Hamiltonian of $S$ is modified while $S$ is in contact with $B$ until $H_{S}^{(N)}$ is reached, the total Hamiltonian reading $H^{(N)}=H_{S}^{(N)}+V+H_{B}$. This amounts to a series of quenches $H_{S}^{(1)}\mapsto\cdots\mapsto H_{S}^{(N)}$, each followed by an equilibration to the corresponding thermal state. 3. 3.

The interaction between $S$ and $B$ is turned off, so that the final Hamiltonian is $H^{(N+1)}=H_{S}^{(N)}+H_{B}$.

Note that this protocol is essentially the same as the optimal protocol for work extraction but without the last step step, which ensured cyclicity. Instead, now $H_{S}^{(N)}$ is fixed, and hence so is the final state of $S$,

[TABLE]

As a consequence, the entropy change is also fixed

[TABLE]

Then, in the spirit of the second law, our aim is to relate $\Delta S$ with the dissipated heat $Q$, and to find the protocol that minimises $Q$.

Firstly, in order to define the average dissipated heat, we invoke the first law of thermodynamics

[TABLE]

Since at the beginning and at the end of the process the Hamiltonian is non-interacting, $\Delta E_{S}$ is simply the change of local energy of the system. Hence $Q=\Delta E_{B}$, which corresponds to the energy dissipated to the bath.

Secondly, since $\Delta E_{S}$ is fixed, it naturally follows that minimizing heat dissipation corresponds to maximizing extractable work. This implies that in the optimal protocol step 2. corresponds to an isothermal process Gallego et al. (2014). In this case, we have that the total work of steps 1.-3. is

[TABLE]

Hence, using (100),

[TABLE]

where we note $\rho_{S}^{(N)}=\mathrm{Tr}_{B}(\omega^{(N)})$ and $\rho_{B}^{(N)}=\mathrm{Tr}_{S}(\omega^{(N)})$. We now make use of the mutual information to write

[TABLE]

and using again $F(\rho,H)-F(\omega_{\beta}(H),H)=T{\rm S}(\rho\|\omega_{\beta}(H))$, we finally obtain

[TABLE]

with $\Delta S_{S}=S(\rho_{S}^{(N)})-S(\rho_{S}^{(0)})$. As expected, the second law is satisfied in the form $-T\Delta S\geq Q$. The protocol minimizing heat dissipation can now be found by minimizing the positive terms in (104). Note however that the second term in (104) is now fixed through $H_{S}^{(N)}$. The only task left is hence to minimise $\Delta F^{(\rm irr)}$ over the Hamiltonian $H_{S}^{(1)}$, a problem which was already solved yielding (51).

Summarizing, the protocol minimizing heat dissipation in the strong coupling regime –given the initial conditions $H^{(0)}=H^{(0)}_{S}+H_{B}$, $\rho^{(0)}=\rho_{S}\otimes\omega_{\beta}(H_{B})$– can be described as: $S$ and $B$ are put in contact through the choice (51), which ensures minimal heat dissipation, and afterwards an isothermal transformation is implemented until the desired Hamiltonian $H_{S}^{(N)}$ is reached.

As we did for work extraction, we now replace $V$ by $gV$ and find the first non-vanishing corrections of the penalizing terms in (104). First of all, to understand the scaling of the term involving the mutual information, let us write the mutual information as

[TABLE]

Since we have $\omega^{(N)}=\omega_{S}^{(N)}\otimes\omega_{B}^{(N)}$ if and only if $g=0$, the function obtains its minimal value [math] only at $g=0$. We thus conclude that the corrections in $g$ will be of the order $g^{2}$ for small $g$. Similarly, as argued in the previous section, the corrections of the relative distances $S(\cdot||\cdot)$ are also of $O(g^{2})$. Hence, without doing any explicit calculation, we can already conclude that,

[TABLE]

where $K_{q}$ is given by

[TABLE]

The computation of $K_{q}$ can be carried out for each particular model of interest. Here, by using our previous considerations, we can provide a general and simple lower bound to it,

[TABLE]

Hence we see that the correction for work extraction in the strong coupling regime also provides a correction to heat dissipation.

Appendix E Carnot engine

Let us now optimize a Carnot engine in the strong coupling regime, by considering the presence of two baths, $B_{c}$ and $B_{h}$, at temperatures $T_{c},T_{h}>0$, respectively. We consider a Carnot cycle consisting of four steps

Coupling to $B_{h}$ + isothermal process+ decoupling from $B_{h}$. We consider an isothermal transformation –following the three steps described in Section D – from $H^{(A)}_{S}$ to $H^{(B)}_{S}$. Note that the state after the transformation is,

[TABLE]

where $H^{(B)}=H^{(B)}_{S}+V_{h}+H_{Bh}$, and $V_{h}$ corresponds to the interaction $SB_{h}$ and $H_{Bh}$ to the local Hamiltonian of $B_{h}$. 2. 2.

Adiabatic expansion. This is represented by a quench from $H^{(B)}_{S}$ to $H^{(C)}_{S}$. 3. 3.

Coupling to $B_{c}$ + isothermal process+ decoupling from $B_{c}$. An isothermal transformation is applied from $H^{(A)}_{S}$ to $H^{(B)}_{S}$. The state after the transformation reads

[TABLE]

with $H^{(D)}=H^{(D)}_{S}+V_{c}+H_{Bc}$. 4. 4.

Adiabatic compression. A quench back to $H^{(A)}_{S}$ is implemented.

With regard to efficiency, such protocols, which only contain one isothermal process with each bath per cycle, are optimal. Intuitively, this follows since every time the working system is coupled to a heat bath, there is unavoidable dissipation, which decreases the efficiency. It is thus optimal to couple to each bath the minimal number of times so that work extraction remains possible. In section E.1 we prove that this intuition is indeed correct.

In the weak coupling regime, the efficiency of the engine is maximised (obtaining the well-known Carnot bound $T=1-\beta_{h}/\beta_{c}$) when

[TABLE]

which ensures no dissipation at any point. Hence, given two Hamiltonians, e.g. $H_{S}^{(B)}$ and $H_{S}^{(D)}$, the engine becomes maximally efficient when the other two Hamiltonians, $H_{S}^{(A)}$ and $H_{S}^{(C)}$, satisfy condition (111). We now look for the corresponding condition in the strong coupling regime. We find convenient to fix $H_{S}^{(B)}$ and $H_{S}^{(D)}$ – so that states $\rho_{1}$ and $\rho_{2}$ are fixed – and optimize over $H_{S}^{(A)}$ and $H_{S}^{(C)}$.

The efficiency of the cycle is defined as

[TABLE]

the minus sign appearing due to the sign convention. By reasonably assuming that the interaction of $S$ with the cold (hot) bath destroys the correlations of $S$ with the other bath (or equivalently, use different baths for each round of the cycle), the initial state before interacting with the hot (cold) bath is $\rho_{1}\otimes\omega_{B_{h}}$ ($\rho_{2}\otimes\omega_{B_{c}}$), where $\omega_{B_{h/c}}=\omega_{\beta}(H_{B_{h/c}})$. Then we can readily apply Eq. (104) obtaining

[TABLE]

where $\Delta{\rm S}={\rm S}(\rho_{1})-{\rm S}(\rho_{2})$ is the entropy loss of the system $S$, $\rho^{(B)}_{B_{h}}=\mathrm{Tr}_{S}\omega_{\beta_{h}}(H^{(B)})$, $\rho^{(D)}_{B_{c}}=\mathrm{Tr}_{S}\omega_{\beta_{c}}(H^{(D)})$ and

[TABLE]

The first law of thermodynamics implies that for a cyclic process $W=-Q_{h}-Q_{c}$ and the efficiency becomes

[TABLE]

where $x_{c/h}$ is the fraction of free energy irreversibly lost during the cold/hot part of the cycle,

[TABLE]

Now, for an engine producing work, we have that $\Delta E_{B_{h}}=-Q_{h}>0$, and hence also $\Delta{\rm S}>0$, which implies that $x_{c/h}>0$. It then follows that the efficiency is lower than the Carnot efficiency, i.e., $\eta<\eta^{{\rm C}}\coloneqq 1-T_{c}/T_{h}$, as expected.

In order to maximise the efficiency $\eta$, one needs to minimise both $x_{c/h}$. Crucially, since $H_{S}^{(B)}$ and $H_{S}^{(D)}$ are fixed, we only need to minimise the term $\Delta F_{c/h}^{\mathrm{(irr)}}$, which was already minimise leading to (51). This implies that the condition (111) naturally generalises in the strong coupling to

[TABLE]

These conditions define the optimal choice of $H_{S}^{(A)}$ and $H_{S}^{(C)}$ in order to maximise $\eta$ in the strong coupling regime.

In order to study the limit of small interactions for the Carnot engine, we first note that $\Delta{\rm S}$ does depend on $g$, as the initial/final state of S before/after being coupled with the heat bath is the reduced of a thermal state. We can then expand it over $g$ obtaining

[TABLE]

where $\Delta{\rm S}^{\rm(weak)}=S(\omega_{\beta}(H_{S}^{(B)}))-S(\omega_{\beta}(H_{S}^{(D)}))$. Then it follows

[TABLE]

and

[TABLE]

where $Q_{c/h}^{\rm(weak)}=T_{c/h}\Delta S^{\rm(weak)}$. Using (115) and the lower bound (108) for $K_{q}^{c/h}$, we obtain the desired corrections to Carnot due to strong coupling.

E.1 Carnot-like protocols are optimal

In this section we prove that the Carnot-like protocols considered above are optimal from the point of view of efficiency. To do this, let us first discuss some preliminaries. Consider an arbitrary cyclic protocol. In every such protocol, there are parts of the protocol where the system remains coupled to one of the heat baths, while multiple quenches from some Hamiltonian $H$ to some other Hamiltonian $H^{\prime}$ are done on the system. The first observation to make is that the protocol can only become more efficient if we replace this part of the protocol by an isothermal reversible process from $H$ to $H^{\prime}$ since such processes are reversible. This shows that optimal protocols will consist only of two kinds of operations:

isothermal reversible processes in contact with one of the baths, 2. 2.

quenches while not being coupled to the baths (adiabatic compression).

Any such protocol is thus composed of $n_{h}$ isothermal process with the hot bath and $n_{c}$ isothermal process with the cold bath, with adiabatic quenches in-between. We can then describe any such protocol by $n_{h}$ pairs of Hamiltonians $(H_{h,c}^{(i)},H_{h,c}^{(f)})$, denoting the initial and final Hamiltonian of the $j$-th IRP with the hot bath, and analogously $n_{c}$ pairs of Hamiltonians $(H_{c,j}^{(i)},H_{c,j}^{(f)})$ for the cold bath. Here, we take the convention that

[TABLE]

Suppose now that a given protocol would have two isothermal processes with one of the baths after each other, only separated by an adiabatic quench. Then this part of the protocol would in general be irreversible and have fixed initial and final states. For concreteness, suppose this would happen with the cold bath and suppose the two initial and final Hamiltonians would be given by $(H_{c,j}^{(i)},H_{c_{j}}^{(f)})$ with $j=1,2$. Then we could replace this part of the protocol by an isothermal reversible process from $H_{c,1}^{(i)}$ to $H_{c,2}^{(f)}$ and the efficiency could only increase. This shows that in optimal protocols, the isothermal reversible process at the two different baths alternate, so that an isothermal process at the hot bath is necessarily followed by an adiabatic quench and an isothermal process at the cold bath (and vice-versa).

Since the protocols have to be cyclic, we already know that $n_{h}=n_{c}$. What remains to be shown is that optimal protocols have $n_{h}=n_{c}=1$. To see this, suppose, for concreteness, that $n_{h}=n_{c}=2$. We will now show that we can always describe such a protocol by two sub-cycles with $n_{h}=n_{c}=1$ that are run sequentially. To do this, let us start the description of the total protocol at the end of the second isothermal at the cold bath, thus starting with state $\omega_{\beta_{c}}(H_{c,2}^{(f)})$. Then the original protocol proceeds as follows:

Adiabatic quench to $H_{h,1}^{(i)}$, 2. 2.

Isothermal process to $H_{h,1}^{(f)}$, 3. 3.

Adiabatic quench to $H_{c,1}^{(i)}$, 4. 4.

Isothermal process to $H_{c,1}^{(f)}$, 5. 5.

Adiabatic quench to $H_{h,2}^{(i)}$, 6. 6.

Isothermal process to $H_{h,2}^{(f)}$, 7. 7.

Adiabatich quench to $H_{c,2}^{(i)}$, 8. 8.

Isothermal process back to $H_{c,2}^{(f)}$.

We have here omitted turning on and off the interaction between the system and bath. Let us now replace this protocol by the following protocol, which has exactly the same efficiency:

Adiabatic quench to $H_{h,1}^{(i)}$, 2. 2.

Isothermal process to $H_{h,2}^{(f)}$ (change here), 3. 3.

Isothermal process from $H_{h,2}^{(f)}$ to $H_{h,1}^{(f)}$, 4. 4.

Adiabatic quench to $H_{c,1}^{(i)}$, 5. 5.

Isothermal process to $H_{c,1}^{(f)}$, 6. 6.

Adiabatic quench to $H_{h,2}^{(i)}$, 7. 7.

Isothermal process to $H_{h,1}^{(f)}$ (change here), 8. 8.

Isothermal process from $H_{h,1}^{(f)}$ to $H_{h,2}^{(f)}$. 9. 9.

Adiabatic quench to $H_{c,2}^{(i)}$, 10. 10.

Isothermal process back to $H_{c,2}^{(f)}$.

The only change that occurred is a splitting off of isothermal processes into two pieces. Due to reversibility of isothermal processes, this does not change the work and heat flows and hence this protocol has the same efficiency. However, we can now see that we have turned the protocol into two simple cycles: One composed of one isothermal in contact with the hot bath from Hamiltonian $H_{h,1}^{(i)}$ to $H_{h,2}^{(f)}$ and an isothermal with the cold bath connecting the Hamiltonians $H_{c,2}^{(i/f)}$. The second cycles consists of the isothermal process in contact with the hot bath connecting the Hamiltonians $H_{h,2}^{(i)}$ and $H_{h,1}^{(f)}$ and an isothermal with the cold bath connecting the Hamiltonians $H_{c,1}^{(i/f)}$ as before. The two cycles are connected by one isothermal from $H_{h,1}^{(f)}$ to $H_{h,2}^{(i)}$ and the same isothermal run backwards. Due two reversibility, these two isothermals ”cancel out” when calculating heat and work. We thus conclude that the total work of the protocol and the total heat absorbed from the hot baths are given by

[TABLE]

where $W_{1}$ and $W_{2}$ denote the work in the individual cycles and similarly for the heat. The total efficiency then yields

[TABLE]

A similar construction can be made for any protocol with $n_{h}=n_{c}>1$. Hence, the efficiency of any such protocol can be bounded as

[TABLE]

where $\eta_{j}$ dentoes the efficiency of the $j$-th sub-cycle. We conclude that optimal protocols are simple Carnot-cycles as analyzed in the previous sections.

Appendix F Power and lower bound on the equilibration time

We have argued in the main text that a dimensional analysis suggests that the equilibration time should scales as $1/g$. In this section, we show that the equilibration time satisfies $\tau\geq C/g$ for some constant $C>0$. By equilibration we mean the process in which the initial expectation value of any operator $A(0)$ evolves in time towards a certain value $\bar{A}$ in which it (approximately) remains.

F.1 Lower bound on the equilibration time for a single equilibration

In order to give a lower bound for the equilibration time let us consider the fastest rate of change of $\left\langle A(t)\right\rangle$, i.e., the quantity

[TABLE]

which is trivially upper bounded by

[TABLE]

This upper bound implies a lower bound on the equilibration time by means of

[TABLE]

In the particular case that $A=H_{S}$, the rate at which the energy of the system changes during the equilibration is bounded by

[TABLE]

with $c=\left\|[H_{S},V]\right\|$. This leads to an equilibration time lower bounded by

[TABLE]

F.2 Lower bound on the equilibration time for the entire cycle of the heat engine

In order to lower bound the equilibration time of the entire cycle, it will be useful to use the bound of the equilibration time of a single equilibration by means of the energy change of the bath

[TABLE]

where now $r\coloneqq\left\|[H_{B},V]\right\|$. Although in general $r$ could scale as $\left\|H_{B}\right\|$, in practice the Hamiltonian of the bath has a locality structure and the commutator is only non-trivial on the degrees of freedom close to the boundary and $r$ is independent of the bath’s size. Let us decompose the cycle in a heat engine described above into two main parts.

• •

Coupling to the cold bath + isothermal reversible process with it + decoupling from it. The system is initially already decoupled from the hot bath and thus the global Hamiltonian is non interacting. After performing the isothermal process with the cold bath, the system is also decoupled from the bath. Hence, the accumulated energy variation of the bath during all the protocol steps is given by

[TABLE]

A lower bound on the time that such part of the protocol requires is given by

[TABLE]

where $r_{c}=\left\|[H_{B_{c}},V_{c}]\right\|$ is the maximum rate at which the bath loses or gains energy, $H_{B_{c}}$ is the Hamiltonian of the cold bath, and $V_{c}$ is the interaction that couples the system to the cold bath.

• •

Coupling to the hot bath + isothermal reversible process with it + decoupling from it. By means of exactly the same argument, the time required to run the second part of the cycle can be lower bounded by

[TABLE]

where $r_{h}=\left\|[H_{B_{h}},V_{h}]\right\|$, $H_{B_{h}}$ is the Hamiltonian of the hot bath, and $V_{h}$ is the interaction that couples the system to the hot bath.

The total run-time of the cycle is then bounded by

[TABLE]

where we have considered that for an engine $|Q_{c}|=|Q_{h}|-|W|$.

F.3 Upper bound on the power of the heat engine

From equation (F.2) we obtain a limit in the power of a heat engine in terms of its efficiency and coupling strength, that is,

[TABLE]

where we have used the definition of efficiency $\eta=W/Q_{h}$ (with $W,Q_{h}>0$ for an engine) and the fact that $1-\eta>0$.

Appendix G Caldeira-Leggett model

In this section, we provide more details about the Caldeira-Leggett (CL) or Ullersma model and the exact spectral density that we are using. The CL-model describes a central harmonic oscillator (the system) coupled to $N$ peripheral modes (constituting the bath), so that the Hamiltonian takes the form

[TABLE]

for which

[TABLE]

here, the coordinates $\{x,p\}$ refer to the system $S$ and $\{x_{\mu},p_{\mu}\}$ to the bath oscillators. $H_{L}$ is the frequency-shift (Lamb-shift) needed to compensate for the distortion induced by the coupling term on the effective potential of the central oscillator. We also note that we introduced the parameter $g$ at hand (see (135) and (139)), in order to quantify the strength of the system-bath interaction. The dynamics of $S$ depends only on the spectral density of the bath $J$, which is defined as

[TABLE]

In the continuum limit, the spectral density is often assumed to be well approximated by a continuous function. A common choice of $J$ is the so called Ohmic spectral density that takes the form

[TABLE]

for frequencies $\omega$ significantly smaller than some cut-off $\Omega>0$. In our work we are interested in large but finite $n$. In order to discretise the above considerations, we assume that the bath frequencies are uniformly distributed,

[TABLE]

$\mu=1,\dots,n$, where $\Omega$ is the highest frequency, and then,

[TABLE]

where we note that $dx\sim\Omega/n$. From (141), we obtain,

[TABLE]

Going back to the discrete indices, and with $\eta=1$, we finally obtain

[TABLE]

Together with the assumption $m_{\mu}=1$, the relations (142) and (145) fully determine the Hamiltonian (135). In fact, the only relevant variable that we consider is the interaction strength $g$ in (135).

Let us now briefly outline how to solve quadratic bosonic systems exactly, which we apply to the Caldeira Leggett model. We refer the reader to Refs. Rivas et al. (2010); Eisert et al. (2010) for more detailed and extensive derivations and explanations. First, we define the vector of canonical coordinates as ${\bf r}=(x,x_{1},x_{2},\dots,x_{n},p,p_{1},\dots,p_{n})^{T}$. Then we can express the total Hamiltonian (135) as

[TABLE]

Now we invoke Williamson’s theorem to write

[TABLE]

where $D={\rm diag}(d_{1},\dots,d_{n+1})$ and the main diagonal elements are given by the strictly positive square roots of the spectrum of $(i\sigma H)^{2}$. The total Hamiltonian can then be written as

[TABLE]

with ${\bf q}=S{\bf r}$. This expression gives rise to

[TABLE]

where ${\bf b}=\Omega^{-1}S{\bf r}$, and

[TABLE]

Let us now define the first and second moments of $\rho:=\rho$ as ${\bf m}_{i}(\rho)=\mathrm{Tr}(\rho r_{i})$, and

[TABLE]

Given these definitions, it can be shown that the time evolution of ${\bf m}(\rho(t))$ and $\gamma(\rho(t))$ under $H$, with $\rho(t)=e^{-iHt}\rho e^{iHt}$,

[TABLE]

where

[TABLE]

Through the exact evolution (152), one can compute, e.g., the energy as a function of time as in Fig. 2.

G.1 Equilibration in the Caldeira-Leggett model

In this section, we aim at understanding under what conditions the expectation value of the operator $A$ equilibrates and how long such equilibration process takes. To do so, we make use of the arguments put forward in Refs. Oliveira et al. (2017); Wilming et al. (2017) in the Hilbert space, but here in the space of modes. The time evolution of an observable quadratic in the canonical coordinates in the Caldeira-Leggett model,

[TABLE]

for an initial state the covariance matrix of which has entries $\gamma_{i,j}(0)=\mathrm{Tr}\left((r_{i}r_{j}+r_{j}r_{i})\rho(0)\right)$ reads

[TABLE]

where $\tilde{A}=\Omega^{T}Q\Omega=\Omega^{T}\sigma S\sigma A\sigma S^{T}\sigma\Omega$ is a matrix associated to the observable, $C_{k,l}=\Omega^{-1}(S\gamma^{\bf r}S^{T}-\textrm{i}\sigma)(\Omega^{-1})^{T}/2$ is the covariance matrix of the initial state, $\tilde{d}_{k}=d_{k}$ if $k\leq L$ and $\tilde{d}_{k}=-d_{k-L}$ for $k>L$, and

[TABLE]

If the observable equilibrates, its equilibrium value is the infinite time average

[TABLE]

Let us restrict us for simplicity to the generic case in which the spectrum of $D$ has $\tilde{d}_{k}=-\tilde{d}_{l}$ if and only if $l=k+L\mathrm{\,mod\,}L$. The equilibrium value of $A$ in such a situation reads

[TABLE]

where the sum in the subindices is taken modulo $L$ and we have identified a Kronecker delta.

Let us now introduce the time signal of an observable $A$ for the initial state $\rho(0)$ as the distance from the equilibrium value of the instantaneous expectation value of $A$ at time $t$

[TABLE]

and for in absence of degeneracies

[TABLE]

It is useful to write the time signal as

[TABLE]

where $\omega_{\alpha}=\omega_{(k,l)}:=\tilde{d}_{k}+\tilde{d}_{l}$ with $k\leq l$ and $l\neq k+L$, which in general forms a set of the $2L^{2}$ different frequencies, and

[TABLE]

is the relevance of each one. The restrictions of the sums over $k$ and $l$ are due to the fact $\tilde{d}_{k}+\tilde{d}_{l}=\tilde{d}_{l}+\tilde{d}_{k}$ and the subtraction of the equilibrium value $\bar{A}$. In this new form (161), the time signal can be seen as the sum of a cloud of points in the complex plane which are initially in the position $v_{\alpha}$ and rotate at angular velocity $\omega_{\alpha}$.

In the same spirit of the works of equilibration in closed quantum systems and in order to define a notion of equilibration we compute the average distance from equilibrium

[TABLE]

where to simplify the calculations we have assumed that the spectrum of $D$ is generic such that $\omega_{\alpha}=\omega_{\alpha^{\prime}}$ if and only if $\alpha=\alpha^{\prime}$.

The average in time of the signal $f(t)$ gives as a notion of to which extent the observable $A$ equilibrates. If $\langle|f(t)|^{2}\rangle_{t}\ll 1$, then the observable takes for most of times the equilibrium value. In contrast, if $\langle|f(t)|^{2}\rangle_{t}\simeq O(f(0))$ then the system is most of times out of equilibrium. Here we assume that the system equilibrates and hence

[TABLE]

The above condition (164) implies some type of synchronization of the initial phases of the complex numbers $v_{\alpha}$. In particular, if the phases of $v_{\alpha}$ were isotropically distributed, then the value of $g(0)\simeq\langle g(t)\rangle_{t}$ instead of (164). To see this, let $v_{\alpha}=|v_{\alpha}|\mathrm{e}^{\textrm{i}\theta_{\alpha}}$ with $\alpha=1,\ldots,d$ be a set of $d$ independent random complex variables with an isotropic probability distribution $p_{\alpha}(r,\theta)=p_{\alpha}(r)=\delta(r-r_{\alpha})$, i. e. the random variable $v_{\alpha}$ has fixed modulus $r_{\alpha}$ and a random phase $\theta_{\alpha}$. Then, the variance of the random variable is given by

[TABLE]

where we have used the fact that the variance of a sum of independent random variables is the sum of variances and the first moments $<v_{\alpha}>=0$. In other words, if the phases are random, the typical value of $f(0)=\sum_{\alpha}v_{\alpha}$ will be of the order

[TABLE]

Thus, the relaxation to equilibrium has to be understood as the dephasing process of the set of points $v_{\alpha}$ in the complex plane. Initially, the points $v_{\alpha}$ are “more or less” synchronized in phase, as time runs, they separate each other due to their different angular velocities $\omega_{\alpha}$. Once they have completely dephased and have formed an isotropic cloud, the system is at equilibrium. As argued in Ref. Oliveira et al. (2017), the estimate of the equilibration time $\tau$ is the inverse of the dispersion of the relevant angular velocities $\omega_{k}$, that is,

[TABLE]

with

[TABLE]

where the relevance $p_{\alpha}=|v_{\alpha}|^{2}/\sum_{\alpha^{\prime}}/|v_{\alpha}|^{2}$ is the normalized relevance of the frequency $\omega_{\alpha}$. In order to understand the behavior of the equilibration time with the strength of the interaction $g$, we need to study how the $\Delta\omega$, and specifically the matrix-elements $|\tilde{A}_{k,l}|$ and $|C_{k,l}|$, change with $g$. In particular, we study the scaling of their dispersion in $\omega$ of $|\tilde{A}_{k,l}|$ and $|C_{k,l}|$ for different $g$’s in the Caldeira-Leggett model taking $A=H_{S}$ and find that both scale as $g^{2}$. This together with Eq. (167) sets a time-scale which behaves as

[TABLE]

This is numerically confirmed in Fig. 3, where the time of equilibration is plot with respect to $1/g^{2}$ in the Caldeira-Leggett model. This supports the idea that the underlying mechanism of equilibration in integrable models is also dephasing.