Work fluctuations in slow processes: quantum signatures and optimal control

TL;DR

This paper reveals that quantum coherence causes violations of classical work fluctuation-dissipation relations in slow processes, leading to non-Gaussian work distributions and a fundamental trade-off between dissipation and fluctuations, addressed through a quantum geometric framework.

Contribution

It introduces a quantum generalisation of the work fluctuation-dissipation relation and develops a geometric method to optimize the balance between dissipation and fluctuations in slow quantum processes.

Findings

Quantum coherence leads to violations of classical FDR.

Work distribution becomes non-Gaussian due to quantum effects.

A quantum geometric framework optimizes the dissipation-fluctuation trade-off.

Abstract

An important result in classical stochastic thermodynamics is the work fluctuation--dissipation relation (FDR), which states that the dissipated work done along a slow process is proportional to the resulting work fluctuations. Here we show that slowly driven quantum systems violate this FDR whenever quantum coherence is generated along the protocol, and derive a quantum generalisation of the work FDR. The additional quantum terms in the FDR are found to lead to a non-Gaussian work distribution. Fundamentally, our result shows that quantum fluctuations prohibit finding slow protocols that minimise both dissipation and fluctuations simultaneously, in contrast to classical slow processes. Instead, we develop a quantum geometric framework to find processes with an optimal trade-off between the two quantities.

Equations

$$W_{\mathrm{diss}}=\frac{1}{2}\beta {\sigma }_w^2.$$

$${\sigma }_w^2=2{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_2\text{Tr}\left({\dot{H}}_{t_1}\overleftarrow{P}(t_1,t_2)\left[{\mathbb{S}}_{{\varrho }_{t_2}}({\dot{H}}_{t_2})\right]\right),$$

$$\displaystyle\mathbb{S}_{\varrho}(A):=\frac{1}{2}\{\varrho,\Delta_{\varrho}A\},\ \ \$$

$${\sigma }_w^2\simeq -2{\int }_0^{\tau }dt\text{Tr}\left({\dot{H}}_t{\mathcal{L}}_t^{+}\left[{\mathbb{S}}_{{\pi }_t}({\dot{H}}_t)\right]\right).$$

$$\displaystyle\mathscr{L}^{+}_{t}[A]:=\int^{\infty}_{0}d\nu\ e^{\nu\mathscr{L}_{t}}\big{[}\pi_{t}\,\text{Tr}\left(A\right)-A\big{]},$$

$$W_{\mathrm{diss}}=-\beta {\int }_0^{\tau }dt\text{Tr}\left({\dot{H}}_t{\mathcal{L}}_t^{+}\left[{\mathbb{J}}_{{\pi }_t}({\dot{H}}_t)\right]\right).$$

$${\mathbb{J}}_{\varrho }(A):={\int }_0^{1}da{\varrho }^a{\Delta }_{\varrho }A{\varrho }^{1-a}.$$

$$\frac{1}{2}\beta {\sigma }_w^2=W_{\mathrm{diss}}+Q_w,$$

$$W_{\mathrm{diss}}\le \frac{1}{2}\beta {\sigma }_w^2.$$

$${\mathcal{L}}_t[{\varrho }_t]=({\pi }_t-{\varrho }_t)\Gamma ,$$

$$\Delta F=-{\beta }^{-1}\ln ⟨e^{-\beta w}⟩=\sum _{k=1}^{\infty }\frac{(-\beta {)}^{k-1}}{k!}{\kappa }_w^{(k)}.$$

$$\sum _{k=3}^{\infty }\frac{(-\beta {)}^{k-1}}{k!}{\kappa }_w^{(k)}=Q_w\ge 0.$$

$${\Lambda }_{ij}({\vec{\lambda }}_t):=-\frac{\beta }{2}\text{Tr}\left(X_i{\mathcal{L}}_t^{+}\left[{\mathbb{S}}_{{\pi }_t}(X_j)\right]+X_j{\mathcal{L}}_t^{+}\left[{\mathbb{S}}_{{\pi }_t}(X_i)\right]\right).$$

$${\xi }_{ij}({\vec{\lambda }}_t):=-\frac{\beta }{2}\text{Tr}\left(X_i{\mathcal{L}}_t^{+}\left[{\mathbb{J}}_{{\pi }_t}(X_j)\right]+X_j{\mathcal{L}}_t^{+}\left[{\mathbb{J}}_{{\pi }_t}(X_i)\right]\right).$$

$${\mathcal{C}}_{\alpha }:=\alpha {\tilde{\sigma }}_w^2+(1-\alpha )W_{\mathrm{diss}}\text{mbox}for\alpha \in [0,1],$$

$$P(w)=\sum _{n,m}\delta [w-{ϵ}_m(\tau )+{ϵ}_n(0)]p_{m|n}{p}_n,$$

$$G(\lambda )$$

$$\displaystyle\langle w\rangle=\int^{\tau}_{0}dt\ \text{Tr}_{\scriptscriptstyle S}\big{(}\dot{H}_{\scriptscriptstyle S}(t)\varrho_{\scriptscriptstyle S}(t)\big{)},$$

$$⟨w^2⟩=2\text{Re}{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_2{\text{Tr}}_{SB}\left({\dot{H}}_S^H(t_1){\dot{H}}_S^H(t_2){\pi }_{SB}(0)\right).$$

$$\displaystyle\dot{\varrho}_{\scriptscriptstyle S}(t):=-\frac{i}{\hbar}\text{Tr}_{\scriptscriptstyle B}\big{(}[H_{\scriptscriptstyle SB}(t),\varrho_{\scriptscriptstyle SB}(t)]\big{)}=\mathscr{L}_{t}[\varrho_{\scriptscriptstyle S}(t)],$$

$${\varrho }_{SB}(t)\simeq {\varrho }_S(t)\otimes {\pi }_B,$$

$$⟨A_S^H(t^{\prime })B_S^H(t)⟩$$

$$\displaystyle=\text{Tr}_{\scriptscriptstyle SB}\left(U^{\dagger}(t^{\prime},0)\,A_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t^{\prime})}\,U(t^{\prime},0)\,U^{\dagger}(t,0)\,B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,U(t,0)\,\pi_{\scriptscriptstyle SB}(0)\right),$$

$$\displaystyle=\text{Tr}_{\scriptscriptstyle SB}\left(A_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t^{\prime})}\,U(t^{\prime},t)\,B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle SB}(t)\,U^{\dagger}(t^{\prime},t)\right),$$

$$\displaystyle=\text{Tr}_{\scriptscriptstyle S}\big{(}A_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t^{\prime})}\,\text{Tr}_{\scriptscriptstyle B}\big{(}U(t^{\prime},t)\,B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle SB}(t)U^{\dagger}(t^{\prime},t)\big{)}\big{)},$$

$$\displaystyle\langle A^{H}_{\scriptscriptstyle S}(t+\tilde{t})B^{H}_{\scriptscriptstyle S}(t)\rangle=\text{Tr}_{\scriptscriptstyle S}\big{(}A_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t+\tilde{t})}\ \chi_{\scriptscriptstyle S}(\tilde{t})\big{)}$$

$$\displaystyle\left\{\begin{array}[]{l}\chi_{\scriptscriptstyle SB}(\tilde{t})=U(t+\tilde{t},t)\,B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle SB}(t)U^{\dagger}(t+\tilde{t},t),\\ \chi_{\scriptscriptstyle S}(\tilde{t})=\text{Tr}_{\scriptscriptstyle B}\big{(}\chi_{\scriptscriptstyle SB}(\tilde{t})\big{)}.\end{array}\right.$$

$$\frac{d}{d\tilde{t}}{\chi }_{SB}(\tilde{t}):=-\frac{i}{\hslash }[H_{SB}(\tilde{t}),{\chi }_{SB}(\tilde{t})],$$

$$\displaystyle\chi_{\scriptscriptstyle S}(\tilde{t})=\overleftarrow{P}(t+\tilde{t},t)[B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle S}(t)],$$

$$\displaystyle\overleftarrow{P}(t_{1},t_{2})=\overleftarrow{\mathcal{T}}\text{exp}\bigg{(}\int_{t_{2}}^{t_{1}}d\nu\ \mathscr{L}_{\nu}\bigg{)},$$

$$\displaystyle\langle A^{H}_{\scriptscriptstyle S}(t+\tilde{t})B^{H}_{\scriptscriptstyle S}(t)\rangle=\text{Tr}_{\scriptscriptstyle S}\big{(}A_{\scriptscriptstyle S}\ \overleftarrow{P}(t+\tilde{t},t)[B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle S}(t)]\big{)}.$$

$$\displaystyle\langle w^{2}\rangle=2\text{Re}\int^{\tau}_{0}dt_{1}\int^{t_{1}}_{0}dt_{2}\ \text{Tr}_{\scriptscriptstyle S}\big{(}\dot{H}_{\scriptscriptstyle S}(t_{1})\ \overleftarrow{P}(t_{1},t_{2})[\dot{H}_{\scriptscriptstyle S}(t_{2})\varrho_{\scriptscriptstyle S}(t_{2})]\big{)}.$$

Full text

Thermodynamics traditionally deals with macroscopic systems at thermal equilibrium, and its laws relate averages of quantities such as work and heat. When bringing the theory to the microscale, fluctuations become significant and can no longer be neglected with respect to average quantities. As a consequence, a stochastic description of thermodynamic processes is needed, which has triggered enormous attention to the understanding of work (and heat) fluctuations  Esposito et al. (2009); Jarzynski (2011); Seifert (2012); Campisi et al. (2011). In the regime of slow but finite-time classical processes, work fluctuations are governed by a single relation, known as the work fluctuation-dissipation-relation (FDR) Hermans (1991); Wood et al. (1991); Hendrix and Jarzynski (2001); Jarzynski (1997):

[TABLE]

Here, $\sigma_{w}^{2}\equiv\langle w^{2}\rangle-\langle w\rangle^{2}$ is the variance of the work distribution $P(w)$ and $W_{\rm diss}\equiv\langle w\rangle-\Delta F\geq 0$ the average dissipated work along the protocol, i.e. the difference between average work done and the change of equilibrium free energy $\Delta F$, which is always non-negative due to the second law, and $\beta=1/k_{B}T$ with $T$ the inverse temperature of the environment. The work FDR (1) is one of the pillars of classical stochastic thermodynamics; it shows that near equilibrium work fluctuations are responsible for dissipation, and conversely that any optimal slow process that minimises dissipation will subsequently minimise the fluctuations Crooks (2007); Sivak and Crooks (2012). For many slow classical processes the work distribution $P(w)$ is Gaussian Speck and Seifert (2004); Chvosta et al. (2007); Šubrt and Chvosta (2007); Speck (2011); Hoppenau and Engel (2013), and if the process also fulfils Jarzynski’s equality then this immediately implies Eq. (1) Jarzynski (1997).

For quantum systems, developing a definition of work and understanding how quantum effects influence its statistics has raised much attention recently Talkner et al. (2007, 2009); Jarzynski et al. (2015); Talkner and Hänggi (2016); Miller and Anders (2017); Perarnau-Llobet et al. (2017); Solinas et al. (2017); Åberg (2018); Bäumer et al. (2018); Holmes et al. (2019); Wu et al. (2019). Previous studies on the non-classicality of work distributions have considered the emergence of quasiprobabilities due to weak measurement Allahverdyan (2014a); Solinas and Gasparinetti (2015a); Potts (2019), contextuality Lostaglio (2018), and violations of macrorealism Blattmann and Mølmer (2017); Miller and Anders (2018a). Despite the wealth of research on this topic, a quantum generalisation of (1) has not been addressed.

Based initially on the Two-Projective-Measurement (TPM) distribution $P(w)$  Talkner et al. (2007, 2009); Esposito et al. (2009), in this article we derive a quantum work FDR and find that it differs from (1) through an additional contribution arising due to quantum fluctuations generated along the protocol. This extra term is positive-definite implying that slow quantum processes are governed by the inequality $W_{\rm diss}\leq\beta\sigma_{w}^{2}/2$, with equality obtained when no coherences in energy are created during the dynamics. We further demonstrate that the extra quantum term in the FDR leads to a non-Gaussian $P(w)$, and show that the same quantum FDR is also valid for work distributions accessed from weak measurements of the system.

While quantum work fluctuations are of fundamental interest, understanding their behaviour also provides a method for minimising them in practical implementations. Indeed, the design of reliable and minimally-dissipative thermodynamic engines is of utmost importance in quantum thermodynamics. In the regime of slow processes, the minimisation of dissipation can be obtained using techniques from differential geometry: one can equip the thermodynamic state space with a Riemannian metric Weinhold (1975); Ruppeiner (1979), and optimal protocols can be found by calculating the associated geodesics Nulton et al. (1985); Salamon and Berry (1983); Crooks (2007); Sivak and Crooks (2012); Zulkowski et al. (2013); Zulkowski and DeWeese (2015); Bonança and Deffner (2018); Scandi and Perarnau-Llobet (2019). Here, we show that the quantum work fluctuations can also be related to a Riemannian metric. However, due to quantum modifications this new metric only coincides with the metric responsible for minimising dissipation in the classical commutative regime. While this result rules out protocols that simultaneously minimise both $W_{\rm diss}$ and $\sigma_{w}$ for quantum coherent processes, our framework can be used to find optimal trade-offs between dissipation and fluctuations.

These results are derived under three main assumptions: (i) the coupling between system and bath is weak, (ii) the system reaches thermal equilibrium when interacting with the bath, (iii) the driving is slow, so that we can expand the magnitudes of interest in the driving velocity and keep only leading terms. Under these assumptions, we now derive a quantum version of the FDR in Eq. (1).

The quantum work FDR. We study the thermodynamics of an open quantum system S coupled to a thermal bath B with total Hamiltonian $H_{\scriptscriptstyle SB}(t)=H_{\scriptscriptstyle S}(t)+H_{\scriptscriptstyle B}+V_{\scriptscriptstyle SB}$, where $H_{\scriptscriptstyle S}(t)=H_{\scriptscriptstyle S}(t)\otimes\mathbb{I}_{\scriptscriptstyle B}$ is the driven system Hamiltonian and $V_{\scriptscriptstyle SB}$ a small but finite coupling Hamiltonian. We take a finite-time interval $t\in[0,\tau]$ and consider processes where the two system Hamiltonian endpoints are fixed, $H_{\scriptscriptstyle S}(0)=H_{0}$ and $H_{\scriptscriptstyle S}(\tau)=H_{\tau}$. We assume that the initial density matrix of S and B is a product $\varrho_{\scriptscriptstyle SB}(0)=\pi_{\scriptscriptstyle S}(0)\otimes\pi_{\scriptscriptstyle B}(0)$ where $\pi_{\scriptscriptstyle S}(0)=e^{-\beta H_{\scriptscriptstyle S}(0)}/Z_{\scriptscriptstyle S}(0)$ and $\pi_{\scriptscriptstyle B}=e^{-\beta H_{\scriptscriptstyle B}}/Z_{\scriptscriptstyle B}$ are the respective Gibbs states for the bare system and bath. The compound system evolves as $\varrho_{\scriptscriptstyle SB}(t)=U(t)\,\varrho_{\scriptscriptstyle SB}(0)\,U^{\dagger}(t)$ with the time-ordered exponential U(t)=\overleftarrow{\mathcal{T}}{\text{exp}}\big{(}-(i/\hbar)\int^{t}_{0}dt^{\prime}\ H_{\scriptscriptstyle SB}(t^{\prime})\big{)}. Work is required to perform $U(t)$, and because only the system Hamiltonian changes in time while coupling is weak, this work can be associated with work on the system alone Talkner et al. (2009). The work statistics can be defined via the TPM scheme, where a projective energy measurement of the total Hamiltonian is performed at the beginning, $H_{\scriptscriptstyle SB}(0)$, and the end, $H_{\scriptscriptstyle SB}(\tau)$, of the process, with the energy differences measured identified as the work values $w$. From the statistics, the work distribution can then be constructed and becomes $P(w)={1\over 2\pi}\int d\lambda\ e^{-i\lambda w}\,G(\lambda)$ with a moment generating function $G(\lambda)=\text{Tr}_{\scriptscriptstyle SB}\left(U^{\dagger}({\tau})e^{i\lambda H_{\scriptscriptstyle SB}({\tau})}U({\tau})e^{-i\lambda H_{\scriptscriptstyle SB}({0})}\varrho_{\scriptscriptstyle SB}(0)\right)$ Talkner et al. (2007, 2009); Esposito et al. (2009), which directly gives the work moments via \langle w^{k}\rangle=(-i)^{k}(d^{k}/d\lambda^{k})G(\lambda)\big{|}_{\lambda=0}. While in the following we will use the TPM work distribution to establish the quantum FDR, we show in Appendix C that our results are also valid for alternative work distributions based on weak measurements Allahverdyan (2014a); Solinas and Gasparinetti (2015a); Potts (2019).

From now on we shall use the more compact notation $X_{t}\equiv X_{\scriptscriptstyle S}(t)$, with $X=\varrho,H,\pi$ and denote $\text{Tr}\left(.\right)$ as the trace over the system degrees of freedom. In general, the reduced dynamics of the system can be written as \dot{\varrho}_{t}=-{i\over\hbar}\,\text{Tr}_{\scriptscriptstyle B}\big{(}[H_{\scriptscriptstyle SB}(t),\varrho_{\scriptscriptstyle SB}(t)]\big{)}=\mathscr{L}_{t}[\varrho_{t}]. Here, we will assume that the system follows an adiabatic Markovian master equation with a unique instantaneous steady state given by the thermal state at each $t\in[0,\tau]$: $\mathscr{L}_{t}[\pi_{t}]=0,$ with $\pi_{t}=e^{-\beta H_{t}}/Z_{t}$ (a precise form of $\mathscr{L}_{t}[\varrho_{t}]$ is presented in Appendix D). This is well–justified whenever the bath dynamics are fast compared to the driving rate of the system Hamiltonian Albash et al. (2012); Dann et al. (2018), and the coupling between S and B is weak enough to satisfy the Born-Markov approximation and the rotating wave approximation Breuer et al. (2002). Importantly, under these assumptions the TPM statistics can be determined by unravelling the master equation in terms of quantum jump trajectories Silaev and Heikkil (2014); Liu (2016); Liu and Xi (2016). These trajectories can then be accessed via local measurements of a quantum detector Pekola et al. (2013), circumventing the need to perform global energy measurements. Under these assumptions, we show in Appendix A that the work fluctuations $\sigma^{2}_{w}\equiv\langle w^{2}\rangle-\langle w\rangle^{2}$ are given by

[TABLE]

where $\overleftarrow{P}(t_{1},t_{2})=\overleftarrow{\mathcal{T}}\text{exp}\left(\int_{t_{2}}^{t_{1}}d\nu\ \mathscr{L}_{\nu}\right)$ is the propagator for the Lindbladian, and we have introduced the linear mapping

[TABLE]

with $\Delta_{\varrho}A=A-\text{Tr}\left(A\varrho\right)$ and $\{,\}$ denoting the anticommutator. We now assume that the total time $\tau$ of the process is large with respect to the time scale(s) of thermalisation, which are encoded in $\mathscr{L}_{t}$. Since the two endpoints of the trajectory are fixed at $H_{0}$ and $H_{\tau}$, one has $\dot{H}_{t}\propto\tau^{-1}$. In this case, we can expand the relevant expressions in terms of $\tau^{-1}$ and keep the leading orders, which we refer to as the slow driving regime. This assumption allows us to further simplify Eq. (2) in Appendix B, using techniques similar to the ones developed in Mandal and Jarzynski (2016) for classical systems. To first order in $\tau^{-1}$ the work fluctuations are

[TABLE]

Note that the integrand is proportional to $\tau^{-2}$, and so for the whole integral $\sigma^{2}_{w}\propto\tau^{-1}$ as desired. In Eq. (4), we have introduced the so-called Drazin inverse $\mathscr{L}^{+}_{t}$ of the Lindblad operator $\mathscr{L}_{t}$ Boullion and Odell (1971); Scandi and Perarnau-Llobet (2019). This inverse is defined as

[TABLE]

and satisfies three conditions Scandi and Perarnau-Llobet (2019):  (i) commutation with the Lindbladian, i.e. $\mathscr{L}_{t}\mathscr{L}^{+}_{t}[A]=\mathscr{L}^{+}_{t}\mathscr{L}_{t}[A]=A-\pi_{t}\text{Tr}\left(A\right)$,  (ii) invariance of the thermal state, i.e. $\mathscr{L}^{+}_{t}[\pi_{t}]=0$,  and (iii) tracelessness, i.e. $\text{Tr}\left(\mathscr{L}^{+}_{t}[A]\right)=0$.

An expression similar to Eq. (4) describes the dissipated work, $W_{\rm diss}$, in slow quantum processes Cavina et al. (2017); Scandi and Perarnau-Llobet (2019)

[TABLE]

Note, that in place of $\mathbb{S}_{\pi_{t}}$ in Eq. (4) the map $\mathbb{J}_{\pi_{t}}$ appears, with

[TABLE]

In the special case that $A$ commutes with $\varrho$ the maps $\mathbb{S}_{\varrho}(A)$ and $\mathbb{J}_{\varrho}(A)$ both reduce to $\mathbb{S}_{\varrho}(A)=\varrho\,\,\Delta_{\varrho}A=\mathbb{J}_{\varrho}(A)$.

Taking the expressions for work fluctuations, $\sigma^{2}_{w}$, and the dissipated work, $W_{\rm diss}$, together, we obtain the quantum work FDR:

[TABLE]

where $Q_{w}=\beta\int^{\tau}_{0}dt\ \text{Tr}\left(\dot{H}_{t}\,\,\mathscr{L}^{+}_{t}\left[(\mathbb{J}_{\pi_{t}}-\mathbb{S}_{\pi_{t}})(\dot{H}_{t})\right]\right)$ is a quantum correction coming from the difference between the maps $\mathbb{S}_{\varrho}(A)$ and $\mathbb{J}_{\varrho}(A)$.

In Appendix D we show that $Q_{w}\geq 0$, with equality if and only if $[\dot{H}_{t},H_{t}]=0$ for $\beta>0$ and $\forall t\in[0,\tau]$. This implies that for slow quantum processes with $[H_{t},\dot{H}_{t}]\not=0$ the classical FDR in Eq. (1) breaks down and the work fluctuations are in fact greater than dissipation. In general, one has an inequality:

[TABLE]

We can then interpret the quantum work FDR (8) as follows: during a slow process where the state remains close to a thermal state $\pi_{t}$, the work fluctuations $\beta\sigma^{2}_{w}/2$ can be divided into two positive contributions: a thermal contribution $W_{\rm diss}$, which arises from the thermal fluctuations in $\pi_{t}$, and a purely quantum contribution $Q_{w}$, which appears whenever quantum fluctuations are created in the dynamics as $[\pi_{t},\dot{H}_{t}]\neq 0$.

Let us rewrite $Q_{w}=\beta\int^{\tau}_{0}dt\mathcal{I}_{t}(\pi_{t},\dot{H}_{t})$ where we have introduced the dynamical skew information

$\mathcal{I}_{t}(\pi_{t},A):=\text{Tr}\left(A\,\,\mathscr{L}^{+}_{t}\left[(\mathbb{J}_{\pi_{t}}-\mathbb{S}_{\pi_{t}})(A)\right]\right)$ for an arbitrary observable $A$. To further elaborate the idea that $Q_{w}$ measures the quantum work fluctuations, for now suppose S evolves under a perfectly thermalising map with a single time-scale $1/\Gamma$, i.e. the Lindbladian satisfies

[TABLE]

which has the Drazin inverse $\mathscr{L}_{t}^{+}(.)=(\text{Tr}\left(.\right)\pi_{t}-\mathbb{I}(.))/\Gamma$. In this case, $\mathcal{I}_{t}(\pi_{t},A)$ becomes proportional to the average Wigner-Yanase-Dyson skew information Hansen (2008); Frérot and Roscilde (2016); Miller and Anders (2018b): $\mathcal{I}_{t}(\pi_{t},A)=-\frac{1}{2\Gamma}\int^{1}_{0}da\,\,\text{Tr}\left([A,\pi^{a}]\,[A,\pi^{1-a}]\right)$ which can be understood as a measure of quantum uncertainty in the observable $A$ Luo (2005): it is positive and vanishes iff $[A,\pi_{t}]=0$, reduces to the usual variance for pure $\pi_{t}=\ket{\psi}\bra{\psi}$, and decreases under classical mixing. For more general Lindbladians, $\mathcal{I}_{t}(\pi_{t},A)$ also takes into account the presence of different timescales of thermalisation through the additional dependence on $\mathscr{L}_{t}^{+}$. Summarising, in Eq. (8) we can interpret $Q_{w}$ as a measure of the time-integrated quantum fluctuations in the power $\dot{H}_{t}$.

Non-Gaussianity of the work distribution. Here we show that these quantum coherences necessarily lead to a non-Gaussian shape of the TPM work distribution $P(w)$. For this $P(w)$ the Jarzynski equality holds Talkner et al. (2009), which relates the change in equilibrium free energy to the cumulants of work done on the system that are computed from $P(w)$:

[TABLE]

Here $\kappa^{(k)}_{w}$ are the cumulants of work, with $\kappa^{(1)}_{w}=\langle w\rangle$ and $\kappa^{(2)}_{w}=\sigma^{2}_{w}$. After rearranging terms in (11) and combining this with the quantum FDR (8), we find

[TABLE]

In fact, as we have seen, $Q_{w}$ vanishes iff $[\dot{H}_{t},H_{t}]=0$ for

$\forall t\in[0,\tau]$. Since a Gaussian work distribution has zero cumulants for $k\geq 3$, we conclude that $P(w)$ necessarily becomes non-Gaussian whenever the process generates coherences of the power operator with respect to the instantaneous Hamiltonian. This contrasts with the classical expectation that slow processes lead to Gaussian work distributions Hendrix and Jarzynski (2001); Speck and Seifert (2004). The equality (12) further demonstrates that measuring the work cumulants can provide a direct witness of quantum fluctuations in power.

Thermodynamic geometry and optimal paths. Now that we have established a relationship between work dissipation and fluctuations, we are in a position to determine optimal protocols. In order to find protocols with minimal fluctuations, one can take a geometric approach similar to Scandi and Perarnau-Llobet (2019); Crooks (2007); Sivak and Crooks (2012).

Considering a decomposition of the system Hamiltonian of the form $H_{t}=X_{0}+\vec{\lambda}_{t}\cdot\vec{X}$, where $\vec{\lambda}_{t}=(\lambda_{1}(t),\lambda_{2}(t),...)$ is the vector of scalar controllable parameters and

$\vec{X}=\partial H_{t}/\partial\vec{\lambda}_{t}=(X_{1},X_{2},...)$ are the corresponding generalised conjugate forces. Then, Eq. (4) can be recast in the form

$\sigma_{w}^{2}={2\over\beta}\int_{0}^{\tau}dt\,\ \left[\frac{d\vec{\lambda}_{t}}{dt}\right]^{T}\bm{\Lambda}(\vec{\lambda}_{t})\left[\frac{d\vec{\lambda}_{t}}{dt}\right]$, where $\bm{\Lambda}(\vec{\lambda}_{t})$ has elements

[TABLE]

It follows that since the rate of dissipated work and dynamical skew information are both positive, $\bm{\Lambda}(\vec{\lambda}_{t})$ is a positive-definite matrix. Since $\bm{\Lambda}(\vec{\lambda}_{t})$ is also symmetric and depends smoothly on $\pi_{t}$, it induces a Riemannian metric on the space of quantum thermal states Petz (1994). Differential geometry then provides an efficient and systematic approach to find optimal protocols by solving Euler-Lagrange equations for the functional $\sigma_{w}^{2}$ of the curve $\vec{\lambda}_{t}$. Curves of minimal fluctuations are identified as geodesics of constant velocity.

The work–fluctuation metric $\bm{\Lambda}(\vec{\lambda}_{t})$ given in Eq. (13) should be compared to the work–dissipation metric $\bm{\xi}(\vec{\lambda}_{t})$, for which W_{\rm diss}=\int_{0}^{\tau}dt\,\ \big{[}\frac{d\vec{\lambda}_{t}}{dt}\big{]}^{T}\bm{\xi}(\vec{\lambda}_{t})\big{[}\frac{d\vec{\lambda}_{t}}{dt}\big{]}, with elements Scandi and Perarnau-Llobet (2019)

[TABLE]

The two metrics $\bm{\Lambda}(\vec{\lambda}_{t})$ and $\bm{\xi}(\vec{\lambda}_{t})$ coincide whenever the conjugate forces commute ie. $[X_{i},X_{0}]=[X_{i},X_{j}]=0\ \forall i,j$. In this special case both metrics reduce to the classical Fisher-Rao metric over the space of thermal states, multiplied with $k_{B}T$ and an integral relaxation time related to the open system dynamics Sivak and Crooks (2012).

In general, the fluctuation and dissipation metrics differ and hence their corresponding geodesics will no longer coincide, in contrast to slow processes in classical thermodynamics. In other words, for quantum processes, any slow protocol $\vec{\lambda}_{t}^{opt}$ that minimises dissipation will have non-minimal fluctuations, and vice versa. To interpolate between these two extremes, one can resort to minimising the objective function

[TABLE]

where $\alpha$ weights the relative importance of the fluctuations versus dissipation and $\tilde{\sigma}^{2}_{w}=\frac{1}{2}\beta\sigma^{2}_{w}$. The family of metrics minimising $\mathcal{C}_{\alpha}$ for weights $\alpha$ is just the convex sum $\bm{g}_{\alpha}(\vec{\lambda}_{t})=\alpha\,\bm{\Lambda}(\vec{\lambda}_{t})+(1-\alpha)\,\bm{\xi}(\vec{\lambda}_{t})$. In Appendix E we use Euler-Lagrange methods to find the optimal protocol $\lambda_{t}^{\text{opt}}(\alpha)$ that minimises $\mathcal{C}_{\alpha}$ when $\lambda_{t}$ is a one-dimensional control parameter with $H_{t}=X_{0}+\lambda_{t}X$. The optimal velocity takes the form $\dot{\lambda}_{t}^{\text{opt}}(\alpha)\propto\sqrt{\xi(\lambda_{t})+\alpha\,\,\beta\mathcal{I}_{t}(\pi_{t},X)}$ which clearly depends on $\alpha$ due to the presence of quantum coherence. This contrasts with the classical case $[X_{0},X]=0$ where the optimal protocol can be obtained for any $\alpha$ by driving the system at a constant dissipation rate Sivak and Crooks (2012).

Example. Let us illustrate our results with a slowly driven harmonic oscillator, $H_{t}=\hbar\omega_{t}\left(a^{\dagger}_{\omega_{t}}a_{\omega_{t}}+1/2\right)$, connected to a perfectly-thermalising bath described by the master equation Eq. (10). Here $\omega_{t}$ is the time-dependent frequency of the oscillator, and $a_{\omega_{t}}$ and $a^{\dagger}_{\omega_{t}}$ are the frequency-dependent creation and annihilation operators. Taking the time-derivative yields the power operator $\dot{H}_{t}=\hbar\dot{\omega}_{t}(H_{t}/\hbar\omega_{t}+((a^{\dagger}_{\omega_{t}})^{2}+a^{2}_{\omega_{t}})/2)$, which does not commute with the instantaneous Hamiltonian $H_{t}$, i.e. $[H_{t},\dot{H}_{t}]\neq 0$. In Fig. 1(a), we compare the expressions for $W_{\rm diss}$ and $\beta\sigma_{w}^{2}/2$ for a slow linear ramp of $\omega_{t}$, and it can be seen that the curves differ substantially at low temperatures (i.e. high $\beta$), where quantum fluctuations become dominant, and become closer for higher temperatures, where thermal fluctuations dominate and classical behaviour is recovered. The corresponding metrics $\Lambda(\omega_{t})$ and $\xi(\omega_{t})$ along with their difference, $\Lambda(\omega_{t})-\xi(\omega_{t})=\beta\mathcal{I}_{t}(\pi_{t},X)$, are shown in Fig. 1(b) as a function of inverse temperature. As expected, this difference vanishes in the high temperature limit ($\beta\to 0$). In the low temperature regime thermal fluctuations given by the dissipation metric $\xi(\omega_{t})$ decay, while quantum coherences contribute more significantly to the total fluctuations in power that are given by $\Lambda(\omega_{t})$. The details of all these calculations are provided in Appendix F.

Turning to optimisation, we now use the metric $g_{\alpha}(\omega_{t})=\alpha\Lambda(\omega_{t})+(1-\alpha)\xi(\omega_{t})$ associated with Eq. (15) to construct geodesics that interpolate between minimally dissipating and minimally fluctuating protocols (see Appendix F). So-called Pareto fronts can be used to bound the region of allowed protocols Solon and Horowitz (2018). This is illustrated in Fig. 2, where Pareto front curves indicate the trade-off between minimal fluctuation $\left(\beta\sigma^{2}_{w}/2\right)$ and minimal dissipation ($W_{\rm diss}$) for various values of $\beta$. Each curve is obtained by evaluating $\beta\sigma^{2}_{w}/2$ and $W_{\rm diss}$ for the geodesics minimising $\mathcal{C}_{\alpha}$ for all values $\alpha\in[0,1]$. If the classical FDR would hold, each curve would collapse into a single point along the diagonal line $\beta\sigma^{2}_{w}/2=W_{\rm diss}$. The quantum correction moves each Pareto front above this line and expands it from a single point to a curve, parametrised by $\alpha$. As expected, this effect is most significant at low temperatures where quantum fluctuations dominate.

Conclusions: In this article, we have studied the statistics of work in slowly driven open quantum systems interacting with a thermal environment. We have derived a quantum FDR for work as shown in Eq. (8), which generalises the well-known classical FDR given by Eq. (1). This result implies that whenever quantum coherence is generated during the dynamics of a slow protocol, then $W_{\rm diss}<\frac{1}{2}\,\beta\,\sigma_{w}^{2}$, which is a genuinely quantum effect. Let us briefly comment on the generality of our results. While (8) has been derived using the TPM approach with thermal initial conditions, we prove in Appendix C that (8) holds more generally for arbitrary initial states using alternative definitions of work based on weak measurements Clerk (2011); Bednorz et al. (2012); Hofer and Clerk (2016); Hofer (2017); Allahverdyan (2014b); Solinas and Gasparinetti (2015b); Solinas et al. (2017). This follows directly because these measurement schemes give rise to the same work average and variance. The validity of the quantum FDR for various work definitions highlights that the quantum effects captured by $Q_{w}$ stem from the coherent dynamics of the protocol, rather than arise as a result of measurement disturbance or a particular choice of work definition (see discussion in Appendix C).

It is also interesting to discuss how breaking any of the three main assumptions used to derive the quantum FDR –namely (i) slow driving, (ii) thermalisation, and (iii) weak coupling– can affect it. Both (i) and (ii) appear crucial: in Appendix H we compare $W_{\rm diss}$ and $\sigma_{w}^{2}$ for a spin in contact with a bosonic bath and, while we verify the validity of Eq. (8) for sufficiently slow driving, we do find violations of the FDR for faster driving. Regarding assumption (ii), one can demonstrate that the quantum FDR can break down if the system is not close to thermal equilibrium even if the dynamics are slow, as shown in Zawadzki et al. (2019) for closed unitary evolutions. On the other hand, we believe that the quantum FDR can remain valid away from the weak coupling regime (i.e. if (iii) is broken): a step towards proving this hypothesis is done in Appendix G. By using a a discrete model of quasi-isothermal processes Gallego et al. (2014); Perarnau-Llobet et al. (2018), we derive an analogous quantum FDR for a system strongly coupled to a thermal bath.

The quantum FDR also implies that it is fundamentally impossible to simultaneously minimise dissipation and fluctuations in slow coherent quantum processes. In the second part of the paper we have derived a family of metrics whose geodesics interpolate between minimally-dissipative and minimally-fluctuating thermodynamic protocols, and our results unveil a new geometric structure within quantum thermodynamics. A promising platform to observe these effects experimentally are quantum dots Pekola (2015); Ludovico et al. (2016); Alonso et al. (2018) and superconducting qubits Pekola et al. (2016); Cottet and Huard (2018), where slowly driven non-commuting protocols appear as a realistic possibility Hacohen-Gourgy et al. (2016) and proposals for observing TPM work statistics using a calorimeter have been made Pekola et al. (2013). An interesting future direction is to extend the FDR to many-body closed systems Fusco et al. (2014); Zawadzki et al. (2019); Arrais et al. (2019), and to investigate how these genuinely quantum effects can modify the thermodynamic uncertainty relations in non-equilibrium steady states Barato and Seifert (2015); Ptaszyński (2018); Agarwalla and Segal (2018); Guarnieri et al. and FDR’s in other contexts such as quantum transport Averin and Pekola (2010).

Acknowledgements. We thank Sreekanth Manikandan, Irene D’Amico and Augusto Roncaglia for interesting discussions and comments on the manuscript. This research was supported in part by the COST network MP1209 “Thermodynamics in the quantum regime” and by the National Science Foundation under Grant No. NSF PHY-1748958. J.A. acknowledges support from EPSRC (grant EP/R045577/1) and the Royal Society.

Appendix A Work moments for open system dynamics

In the main text we define the fluctuating work done on the system using the Two-Projective-Measurement (TPM) scheme. Let us denote the spectral decomposition of the total system-bath Hamiltonian $H_{\scriptscriptstyle SB}(t)=\sum_{n}\epsilon_{n}(t)\ket{\epsilon_{n}(t)}\bra{\epsilon_{n}(t)}$ and define the time-ordered unitary U(t_{f},t_{i})=\overleftarrow{\mathcal{T}}{\text{exp}}\big{(}-(i/\hbar)\int^{t_{f}}_{t_{i}}dt^{\prime}\ H_{\scriptscriptstyle SB}(t^{\prime})\big{)} generated by variations in the global Hamiltonian. The TPM work distribution $P(w)$ is constructed in terms of observed transitions between energy states Talkner et al. (2007), resulting in

[TABLE]

where $p_{n}=\bra{\epsilon_{n}(0)}\pi_{\scriptscriptstyle SB}(0)\ket{\epsilon_{n}(0)}$ is the initial energy occupation probability while p_{m|n}=\big{|}\bra{\epsilon_{m}(\tau)}U(\tau,0)\ket{\epsilon_{n}(0)}\big{|}^{2} denotes the conditional energy transition probability. After taking a Fourier transform, the expression for the corresponding moment generating function of the work is given by:

[TABLE]

Note that this function implicitly assumes that $\pi_{\scriptscriptstyle SB}(0)$ commutes with the initial Hamiltonian, which is guaranteed since $\pi_{\scriptscriptstyle SB}(0)=e^{-\beta H_{\scriptscriptstyle SB}(0)}/\text{Tr}_{\scriptscriptstyle SB}\left(e^{-\beta H_{\scriptscriptstyle SB}(0)}\right)$ is a global thermal state with negligible coupling between system and bath. From (17) one can compute the first two work moments using \langle w^{k}\rangle=(-i)^{k}(d^{k}/d\lambda^{k})G(\lambda)\big{|}_{\lambda=0}, and a lengthy but straightforward calculation yields the following Suomela et al. (2014):

[TABLE]

where we denote $X^{H}_{\scriptscriptstyle S}(t)=U^{\dagger}(t,0)\,X_{\scriptscriptstyle S}(t)\,U(t,0)$ as operator $X_{\scriptscriptstyle S}(t)$ in the Heisenberg picture. The expressions (18) and  (19) are valid for any global unitary evolution driven by a total Hamiltonian $H_{\scriptscriptstyle SB}(t)=H_{\scriptscriptstyle S}(t)+H_{\scriptscriptstyle B}+V_{\scriptscriptstyle SB}$ with time-dependence only on the system Hamiltonian, and thus $\dot{H}_{\scriptscriptstyle SB}(t)=\dot{H}_{\scriptscriptstyle S}(t)$. Note that for $\langle w^{2}\rangle$ in (19) the Heisenberg picture $\dot{H}_{\scriptscriptstyle S}^{H}(t)$ is used instead of $\dot{H}_{\scriptscriptstyle S}(t)$, and thus the trace is taken over the whole Hilbert space due to the bath-dependence of the unitary $U(t,0)$.

Our goal will be to express the second moment in terms of the reduced density operator \varrho_{\scriptscriptstyle S}(t)=\text{Tr}_{\scriptscriptstyle B}\big{(}\varrho_{\scriptscriptstyle SB}(t)\big{)}, where we define the evolved density operator for the composite state by $\varrho_{\scriptscriptstyle SB}(t)=U(t,0)\pi_{\scriptscriptstyle SB}(0)U^{\dagger}(t,0)$. To do this we now assume that the evolution of the system is of Lindblad form as defined in the main text:

[TABLE]

with $\mathscr{L}_{t}[(.)]$ a time-dependent Lindbladian. Implicit within our assumption for (20) is the Born-Markov approximation, which assumes that the global state remains factorised at all times during the evolution Breuer et al. (2002):

[TABLE]

This assumption is justified only in the weak-coupling regime.

Our goal will now be to use (20) to rewrite (19) in terms of the system degrees of freedom. Let us now consider two hermitian time-dependent operators $A_{\scriptscriptstyle S}(t),B_{\scriptscriptstyle S}(t)$ acting on the system Hilbert space alone. We are concerned with evaluating the two-time correlation function $\langle A^{H}_{\scriptscriptstyle S}(t^{\prime})\,B^{H}_{\scriptscriptstyle S}(t)\rangle$ in Heisenberg picture for $t^{\prime}\geq t$, which can be expressed as follows:

[TABLE]

where in the third line we used $\pi_{\scriptscriptstyle SB}(0)=U^{\dagger}(t,0)\varrho_{\scriptscriptstyle SB}(t)U(t,0)$ and the cyclicity of the trace, while in the final line we used the fact that B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}=B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\otimes\mathbb{I}_{\scriptscriptstyle B}. Setting $\tilde{t}=t^{\prime}-t\geq 0$, a simple change in variables gives

[TABLE]

where

[TABLE]

Now observe that $\chi_{\scriptscriptstyle SB}(\tilde{t})$ is the solution to the following equation of motion:

[TABLE]

with initial condition \chi_{\scriptscriptstyle SB}(0)=B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\,\varrho_{\scriptscriptstyle SB}(t). We now use the Born-Markov approximation (21), which implies that initial condition to (27) factorises according to \chi_{\scriptscriptstyle SB}(0)=B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t)}\rho_{\scriptscriptstyle S}(t)\otimes\pi_{\scriptscriptstyle B}. Given that the initial operator $\chi_{\scriptscriptstyle SB}(0)$ here factorises and obeys the same global equation of motion given by (27) with respect to $\tilde{t}$ as the state $\varrho_{\scriptscriptstyle SB}(t)$, we obtain the following solution after tracing out the bath degrees of freedom:

[TABLE]

where

[TABLE]

is the propagator for the Lindbladian in (20). Combining this with (23) we have

[TABLE]

Setting A_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t_{1})}=\dot{H}_{\scriptscriptstyle S}(t_{1}) and B_{\scriptscriptstyle S}{\color[rgb]{0,0,0}(t_{2})}=\dot{H}_{\scriptscriptstyle S}(t_{2}), and combining this all together gives us an expression for $\langle w^{2}\rangle$ from (19) in terms of the system degrees of freedom:

[TABLE]

Furthermore, the squared average work (18) can be written as follows:

[TABLE]

We now define $\Delta_{\varrho}A=A-\text{Tr}\left(A\varrho\right)$ and combine (31) and (A) to get

[TABLE]

where in the third line we used the fact that $\text{Re}\ \text{Tr}\left(A\overleftarrow{P}(t_{1},t_{2})[B]\right)=(1/2)\text{Tr}\left(A\overleftarrow{P}(t_{1},t_{2})[B+B^{\dagger}]\right)$ for $A=A^{\dagger}$, and in the fourth we introduced the definition for $\mathbb{S}_{\varrho}$ from Eq. (3). This concludes the derivation of Eq. (2) in the main text.

Appendix B Derivation of Eq. (4)

We want to take the slow driving limit of the expression:

[TABLE]

Recalling the definition of \overleftarrow{P}(t_{1},t_{2})=\overleftarrow{\mathcal{T}}{\text{exp}}\big{(}\int_{t_{2}}^{t_{1}}d\nu\mathscr{L}_{\nu}\big{)}, we notice that the trace will decay to zero exponentially fast in

$|t_{1}-t_{2}|\sim\tau$, since $\mathbb{S}_{\varrho_{t_{2}}}(\dot{H}_{t_{2}})$ is traceless. For this reason, we can substitute at first order in $1/\tau$ the varying Liouvillian with the initial one:

[TABLE]

where in the second line we made the substitution $s=t_{1}-t_{2}$. Again, since $s$ will be typically much bigger than the thermalisation timescales, not only we can approximate $t_{1}-s$ with $t_{1}$ in all the expression (since the correction for finite $s$ will be exponentially suppressed), but also we can send the limit of the integration to infinity. Then, equation (35) becomes:

[TABLE]

where in the last step we used the integral expression of the Drazin inverse $\mathscr{L}_{t}^{+}$ in Eq. (5) and the fact that $\mathbb{S}_{\varrho_{t_{1}}}(\dot{H}_{t_{1}})$ is traceless. Finally, at first order in $1/\tau$, we can substitute $\varrho_{t}\simeq\pi_{t}$. This concludes the derivation of Eq. (4) in the main text.

Appendix C Weak measurements of fluctuating work

Our analysis throughout the paper defines work via the standard two-projective measurement (TPM) scheme. This definition of work is typically adopted for states that are initially thermal, as one recovers the usual Jarzynski equality and laws of thermodynamics at the ensemble level Talkner et al. (2009) while remaining consistent with the usual definition of work in the classical regime Jarzynski et al. (2015). Furthermore, since the initial state is diagonal in the energy basis one can neglect the effect of measurement backaction caused by the initial projective energy measurement. However, for states that are non-diagonal the TPM scheme can remove initial coherences due to disturbances caused by the first projective measurement. As a result, one can no longer associate the average TPM work to the change in total energy of the system and bath Bäumer et al. (2018).

In order to characterise work when initial coherences are present, alternative definitions of work based around weak measurement have been proposed that preserve the coherent evolution of the system and bath Allahverdyan (2014b); Solinas and Gasparinetti (2015b); Solinas et al. (2017). In these measurement schemes one may obtain negative quasi-probabilities in the work distribution, signifying uniquely quantum behaviour such as contextuality Lostaglio (2018). In this Appendix we will demonstrate that definitions of work based on weak measurement gives rise to the same fluctuation-dissipation relation Eq. (8), but are now applicable to arbitrary initial states of the system that may be non-diagonal in the energy basis. While different choices of measurement scheme typically lead to different work statistics, any discrepancies between different definitions of quantum work only apply to the moments of third order and higher, which are not relevant to the quantum work FDR Eq. (8). We will therefore demonstrate that our results hold more generally beyond the TPM definition of work.

We first describe a continuous weak measurement scheme that can be implemented by coupling only to the system degrees of freedom. Recall that in the absence of any measurement, the system and bath evolve according to state $\varrho_{\scriptscriptstyle SB}(t)=U(t,0)\varrho_{\scriptscriptstyle SB}(0)U^{\dagger}(t,0)$, with unitary U(t_{f},t_{i})=\overleftarrow{\mathcal{T}}{\text{exp}}\big{(}-i\int^{t_{f}}_{t_{i}}dt^{\prime}\ H_{\scriptscriptstyle SB}(t^{\prime})\big{)} and setting $\hbar=1$. In contrast to the main text, we now assume arbitrary initial conditions for the system state, such that $\varrho_{\scriptscriptstyle SB}(0)=\rho_{\scriptscriptstyle S}(0)\otimes\pi_{\scriptscriptstyle B}$. Here $\rho_{\scriptscriptstyle S}(0)$ may not be thermal or even commute with its Hamiltonian, ie. $[\rho_{\scriptscriptstyle S}(0),H_{\scriptscriptstyle S}(0)]\neq 0$.

In this weak measurement approach, the fluctuations in work can be characterised as the time-integrated fluctuations in the power observable $\dot{H}_{\scriptscriptstyle S}(t)$ Hofer (2017). To determine these fluctuations we couple the system and bath to a detector modelled as a two-level system, initially uncorrelated such that $\varrho_{\scriptscriptstyle SBD}(t)=\varrho_{\scriptscriptstyle SB}(0)\otimes\varrho_{\scriptscriptstyle D}(0)$. The modified Hamiltonian including the interaction with the detector is now given by

[TABLE]

where $\lambda$ is a coupling constant and $\sigma_{z}$ is the Pauli spin-z operator. This generates a new evolution

\tilde{U}(t,0)=\overleftarrow{\mathcal{T}}{\text{exp}}\big{(}-i\int^{t}_{0}dt^{\prime}\ H_{\scriptscriptstyle SBD}(t)\big{)}, and we denote the reduced state of the detector as a function of $\lambda$ by $\varrho_{\scriptscriptstyle D}(t;\lambda)$. By measuring the relative change in phase of the detector one gets the so-called Keldysh-ordered full-counting statistics Clerk (2011); Bednorz et al. (2012); Hofer and Clerk (2016); Hofer (2017):

[TABLE]

where

[TABLE]

where $\dot{H}_{\scriptscriptstyle S}^{H}(t^{\prime})$ is in the Heisenberg picture with respect to the isolated system and bath unitary dynamics. In this approach one interprets (38) as a moment generating function for the fluctuating work, with moments

[TABLE]

One can see that these statistics generally differ from those obtained via the TPM scheme via (17), as highlighted in Hofer (2017). However, we will now prove that the first two moments of (38) are in fact equivalent to those obtained from (17) used in the main text, but valid for any initial system state.

We first use the Magnus expansion for $V_{\lambda}(t)=\exp{(\Omega_{\lambda}(t))}$, where

[TABLE]

The generating function (38) can then be expressed in powers of $\lambda$:

[TABLE]

It should be noted here that the second order term in (41) is skew hermitian and hence does not contribute to the second order term in (42). Finally using the definition of the work moments, we get

[TABLE]

Comparison with (18) and (19) confirms that these first two moments of work (43) and (44) take the same form as the TPM definition used in the main text, but now hold for any initial conditions chosen for the system. With these expressions the proof of the quantum work FDR follows exactly as before for any initial state.

As an alternative to continuous weak measurement, full-counting statistics for work can be obtained through discrete coupling to the Hamiltonian. This gives rise to a different moment generating function Solinas and Gasparinetti (2015b):

[TABLE]

It is straightforward to show that the first and second moments here become

[TABLE]

As shown in Suomela et al. (2014), these expressions are exactly equivalent to the time-integrated expressions (43) and (44) since we assume only the system Hamiltonian to be time-dependent. We remark that another definition of work based on the Margenau-Hill quasi-probability distribution has been proposed by Allahverdyan in Allahverdyan (2014b), which is obtained from an alternative two-point measurement scheme combining both strong and weak energy measurements. Finally, a quasi-probability for work constructed within the quantum histories approach is proposed in Miller and Anders (2017). Again, in these cases one finds the same first two moments of work (46), which equate to (43) and (44). We reiterate that all of these quasi-probabilistic generalisations of the TPM work distribution for arbitrary initial states differ only for third moments and higher, and therefore the particular choice is not relevant for deriving the quantum work FDR, which is a statement only about the first and second work moments.

In conclusion, we have shown that the quantum work FDR Eq. (8) continues to hold for arbitrary initial states using alternative definitions of quantum work based on weak measurement. Therefore Eq. (8) is not only restricted to the TPM definition of work, or a particular choice of initial conditions. This implies that the quantum modifications to the work FDR are a manifestation of the coherent dynamics generated during the slow driving protocol, and are not a signature of quantum measurement effects stemming from a particular choice of measurement scheme. Furthermore, we have presented a continuous weak measurement scheme that can verify Eq. (8) using only local interactions with the system degrees of freedom, as opposed to measurements of the full system-bath Hilbert space.

Appendix D Proof of Eq. (9)

In this Appendix we prove the positivity of the dynamical skew information $\mathcal{I}_{t}(\pi_{t},\dot{H}_{t})$, which under time integration gives the quantum correction $Q_{w}=\beta\int^{\tau}_{0}dt\ \mathcal{I}_{t}(\pi_{t},\dot{H}_{t})$ appearing in Eq. (8). Consider the Hilbert space $M_{d}$ of $d\times d$ complex matrices with Hilbert-Schmidt inner product $\langle A,B\rangle=\text{Tr}\left(B^{\dagger}A\right)$. Then any superoperator $\mathcal{M}(.)$ acting on the elements of this Hilbert space can be expressed as a $d^{2}\times d^{2}$ matrix. The matrix describing $\mathcal{M}(.)$ is positive if $\text{Tr}\left(A^{\dagger}\mathcal{M}(A)\right)\geq 0$ for any $A\in M_{d}$, and we define the adjoint $\mathcal{M}^{\dagger}(.)$ as the superoperator satisfying $\text{Tr}\left(\mathcal{M}^{\dagger}(B^{\dagger})A\right)=\text{Tr}\left(B^{\dagger}\mathcal{M}(A)\right)$. We begin by assuming a generic interaction between system and bath formed by a sum of hermitian operators

[TABLE]

As stated in the main text, we will work in the slow driving regime and assume that the bath dynamics are much faster than the driving rate of the system Hamiltonian. This means that one can neglect any non-adiabatic contributions to the reduced system dynamics Yamaguchi et al. (2017); Dann et al. (2018). In addition the system dynamics are assumed to satisfy detailed balance along with the Born-Markov and rotating-wave approximations Breuer et al. (2002); Alicki and Lendi (2007). When taken together these assumptions result in a time-dependent Markovian master equation describing the system dynamics that can be expressed in a Lindblad form $\mathscr{L}_{t}(.)=-i[H_{t},(.)]+\mathcal{D}_{t}(.)$, and a precise derivation of its structure and regime of validity can be found in Albash et al. (2012). Throughout this derivation we will only be concerned with the structure of the Lindbladian at some fixed point in time. At any time $t$ the time-dependent Lindbladian takes the following form Albash et al. (2012); Yamaguchi et al. (2017); Dann et al. (2018):

[TABLE]

where $\gamma_{\alpha\beta}(\omega_{t})$ is a hermitian matrix representing the Fourier transform of the bath correlation function, and we have defined the eigenoperators

[TABLE]

with $H_{t}=\sum_{j}\epsilon_{j}(t)\ket{\epsilon_{j}(t)}\bra{\epsilon_{j}(t)}$ the spectral decomposition of the system Hamiltonian at some fixed point in time. The eigenoperators satisfy

[TABLE]

It then follows that such a Lindbladian has a unique zero eigenvalue corresponding to a thermal fixed point $\mathscr{L}_{t}(\pi_{t})=0$, while all other eigenvalues have a strictly negative real part Rivas and Huelga (2012). This ensures that at each fixed configuration $\mathscr{L}_{t}$, any initial state $\rho$ will converge to the instantaneous fixed point:

[TABLE]

We now observe some important properties of $\mathscr{L}_{t}(.)$ Breuer et al. (2002). Firstly, due to the rotating-wave approximation the unitary and dissipative parts commute:

[TABLE]

Secondly, $\mathscr{L}_{t}(.)$ satisfies the condition of detailed balance, which implies

[TABLE]

Finally, the bath correlation function satisfies the KMS condition and hence

[TABLE]

Now note that the dynamical skew information is a real-valued trace functional, thus it is sufficient to prove positivity of the quantity

[TABLE]

where $A=A^{\dagger}$ is an arbitrary hermitian operator, $\mathscr{L}_{t}^{+}$ the Drazin inverse of the Lindbladian defined in Eq. (5) and

[TABLE]

Here $\mathcal{M}(.)$ represents the difference between the arithmetic and logarithmic matrix means, and is hence a positive superoperator due to the Kubo-Ando inequality Kubo and Ando (1980). Alternatively, we can see this by looking at the spectrum of $\mathcal{M}_{t}(.)$. The eigenvectors are given by the energy state elements $\ket{\epsilon_{i}(t)}\bra{\epsilon_{j}(t)}$, and one finds

[TABLE]

where

[TABLE]

and $p_{i}(t)$ represent the eigenvalues of $\pi_{t}$. In addition, since $\pi_{t}$ commutes with Hamiltonian $H_{t}$ one can verify the commutation relation

[TABLE]

Let us now consider the relation between $\mathcal{M}_{t}(.)$ and the dissipator $\mathcal{D}_{t}(.)$. It is first useful to introduce the following integral representation for the matrix power $\pi^{s}$ for positive $\pi$ Briët and Harremoës (2009):

[TABLE]

with $\mu_{s}(x)$ a positive measure on $(0,\infty)$ that we leave unspecified for convenience. Using detailed balance (D) we get the following:

[TABLE]

Similarly one finds

[TABLE]

This then implies

[TABLE]

Using this one obtains the following:

[TABLE]

where in the second line we used (64), in the third line we used (51) and (55), in the fourth line swapped indices $-\omega_{t}\to\omega_{t}$ and $\alpha\to\beta$, and in the final line used the definition of the adjoint superoperator $\text{Tr}\left(\mathcal{D}_{t}^{\dagger}(B^{\dagger})A\right)=\text{Tr}\left(B^{\dagger}\mathcal{D}_{t}(A)\right)$ and the fact that the second term is self-adjoint. We next introduce a complementary Lindbladian of form

[TABLE]

Notably the real part of the spectrum of $\tilde{\mathscr{L}}_{t}(.)$ coincides with that of the original Lindbladian $\mathcal{L}_{t}(.)$ due to the fact that $\mathcal{U}_{t}(.)$ is skew hermitian. Recalling that by assumption $\pi_{t}$ is a unique fixed point of $\mathcal{L}_{t}(.)$, this implies that $\tilde{\mathscr{L}}_{t}(.)$ also has a unique fixed point $\tilde{\pi}_{t}$. We can thus define a pair of Drazin inverses given by

[TABLE]

for any $A\in M_{d}$. These inverses act according to

[TABLE]

By using (60) and (D) we have

[TABLE]

For any traceless matrix $\{B\ |\ B\in M_{d},\ \text{Tr}\left(B\right)=0\}$, we can combine (D) and (70) to get

[TABLE]

We also define the following superoperator:

[TABLE]

Using (71) one can also see that

[TABLE]

Furthermore, since by assumption the real part of the eigenvalues of the Lindbladian $\mathscr{L}_{t}$ are negative, the same holds true for both Drazin inverses $\mathscr{L}_{t}^{+}$ and $\tilde{\mathscr{L}}_{t}^{+}$. To see this, let us consider any non-zero eigenvalue of $\mathscr{L}_{t}$ such as $z=x+iy$, with $x,y\in\Re e$. For the Drazin inverse $\mathscr{L}_{t}^{+}$, which shares the same eigenvectors as $\mathscr{L}_{t}$, the corresponding eigenvalue is given by $z^{-1}=(x+iy)^{-1}=(x-iy)/(x^{2}+y^{2})$ Boullion and Odell (1971). By assumption $x<0$, and thus all eigenvalues of $\mathscr{L}_{t}^{+}$ must also have a negative real part $\Re e(z^{-1})=x/(x^{2}+y^{2})<0$. We also note that the non-zero eigenvalues of $\tilde{\mathscr{L}}_{t}$ must also have a negative real part since the real part of the spectrum coincides with that of $\mathscr{L}_{t}$. By the same argument as above, this means that the non-zero eigenvalues of $\tilde{\mathscr{L}}^{+}_{t}$ have a negative real part.

As a result, the eigenvalues of $\mathcal{V}_{t}(.)$ must have a positive real part. By Corollary 4.2 of Hassi (2005), a matrix product $XY$ with $Y\geq 0$ is positive if the eigenvalues of $X$ have no negative real part and $XY=YX^{\dagger}$. Since $\mathcal{M}_{t}(.)$ is positive, (73) implies that

[TABLE]

Finally, we return to the trace functional (56). Let us introduce the projection onto the traceless subspace $\mathcal{P}_{T}(A)=A-\text{Tr}\left(A\right)\mathbb{I}/d$. Taking all results together one gets

[TABLE]

where in the second line we used the fact that only traceless elements contribute to the functional due to (59), in the third line we set $\mathcal{P}_{T}(A)=B$ and $\mathcal{M}_{t}^{\dagger}(.)=\mathcal{M}_{t}(.)$, in the fourth line we used (71) and in the final line we used the matrix positivity (74). Since the above holds for any hermitian matrix $A$, we conclude that the dynamical skew information is positive. Under time integration we therefore have $Q_{w}\geq 0$ and inequality Eq. (9) To conclude, we now prove $\mathscr{I}(\pi_{t},A)=0$ if and only if $[A,\pi_{t}]=0$. Without loss of generality we may assume $\text{Tr}\left(A\right)=0$. If $[A,\pi_{t}]=0$ then $\mathcal{M}_{t}(A)=0$ and thus $\mathscr{I}(\pi_{t},A)=0$. On the other hand, let us instead suppose $[A,\pi_{t}]\neq 0$, in which case one necessarily has $\mathcal{M}_{t}(A)\neq 0$ which follows from (59). Now note that the unique fixed point of $\mathcal{V}_{t}$ is $\pi_{t}$. Clearly $\mathcal{M}_{t}(A)\not\propto\pi_{t}$ if $[A,\pi_{t}]\neq 0$, and hence $\mathcal{V}_{t}\mathcal{M}_{t}(A)\neq 0$. Using the positivity of $\mathcal{V}_{t}\mathcal{M}_{t}$ we then rewrite (D) as $\mathscr{I}(\pi_{t},A)=\langle\sqrt{\mathcal{V}_{t}\mathcal{M}_{t}(A)},\sqrt{\mathcal{V}_{t}\mathcal{M}_{t}(A)}\rangle>0$. Therefore $\mathscr{I}(\pi_{t},A)=0\Leftrightarrow[A,\pi_{t}]=0$. This further implies that the total quantum correction $Q_{w}$ vanishes if and only if $[\pi_{t},\dot{H}_{t}]=[H_{t},\dot{H}_{t}]=0$ at all times.

Appendix E Finding geodesics for a single parameter

Consider a single parameter change $\lambda_{0}\to\lambda_{\tau}$ with Hamiltonian $H_{t}=X_{0}+\lambda_{t}X$, where in general $[X_{0},X]\neq 0$. Furthermore, denote the rescaled work fluctuations by $\tilde{\sigma}^{2}_{w}=\frac{1}{2}\beta\sigma^{2}_{w}$. The aim is to minimise the linear objective function

[TABLE]

with respect to protocol $\lambda_{t}$. Using the metrics $\Lambda(\lambda_{t})$ and $\xi(\lambda_{t})$ from the main text in Eq. (13) and Eq. (14), we have

[TABLE]

where $\beta\mathcal{I}_{t}(\pi_{t},X)={\Lambda}({\lambda}_{t})-{\xi}({\lambda}_{t})$ is the dynamical skew information for a single parameter change (which depends on $\lambda_{t}$ through $\pi_{t}$), following from the general definition of $\mathcal{I}_{t}(\pi_{t},A)$ in the main text before Eq. (10). The functional (77) is minimised by the solution to the Euler-Lagrange equation for the cost function C_{\alpha}(\lambda_{t},\dot{\lambda}_{t}):=\dot{\lambda}_{t}^{2}\big{(}\xi(\lambda_{t})+\alpha\beta\mathcal{I}_{t}(\pi_{t},X)\big{)}, i.e.

[TABLE]

which gives

[TABLE]

Solving (78) yields an equation for the optimal velocity of the control parameter for a given $\alpha$:

[TABLE]

One concludes that the optimal velocity is proportional to the following:

[TABLE]

Appendix F Thermodynamic metrics for the harmonic oscillator

We wish to evaluate the following metrics for the single parameter $\omega_{t}$:

[TABLE]

For the harmonic oscillator the Hamiltonian and power operator are given by

[TABLE]

where $n_{t}=a_{\omega_{t}}^{\dagger}a_{\omega_{t}}$ with $a_{\omega_{t}}=\sqrt{m\omega_{t}\over 2\hbar}(x+i{p\over m\omega_{t}})$. The metrics (82) can then be simplified as

[TABLE]

where we introduce $\delta H_{t}=\hbar\omega_{t}((a^{\dagger}_{\omega_{t}})^{2}+a^{2}_{\omega_{t}})/2$ for shorthand. By working in the number basis $H_{t}=\hbar\omega_{t}\sum_{n=0}^{\infty}(n_{t}+\frac{1}{2})\ket{n}\bra{n}$ and using the standard relations $a_{\omega_{t}}^{\dagger}\ket{n_{t}}=\sqrt{n_{t}+1}\ket{n_{t}+1}$ and $a_{\omega_{t}}\ket{n_{t}}=\sqrt{n_{t}}\ket{n_{t}-1}$, a textbook calculation reveals the following expressions:

[TABLE]

where $x=\beta\hbar\omega_{t}$. After plugging these expressions into (85) we get

[TABLE]

These quantities are used to plot Fig. 1(a)-(b) in the main text.

From this we can analyse the limiting behaviour of the metrics in different temperature regimes. In the limit $\beta\to 0$ these metrics reduce to

[TABLE]

so that

[TABLE]

In the limit $\beta\to\infty$ the above metrics reduce to

[TABLE]

Therefore in the high temperature limit ($\beta\rightarrow 0$), the metrics $\xi(\omega_{t})$ and $\Lambda(\omega_{t})$ become equal, indicating vanishing of the quantum fluctuations $Q_{w}$. On the other hand, when the temperature is low ($\beta\rightarrow\infty$), the dissipation metric $\xi(\omega_{t})$ converges to a constant in $\beta$, $\lim_{\beta\to\infty}\xi(\omega_{t})=\hbar/2\Gamma\omega_{t}$, while the fluctuation metric $\Lambda(\omega_{t})$ grows linearly $\Lambda(\omega_{t})\approx\beta\,\hbar^{2}/2\Gamma$ for $\beta\to\infty$, see Fig. 1(b) in the main text.

Given a fixed initial and final frequency $(\omega_{0},\omega_{\tau})$, we can now use (91) and (92) to minimise the objective function

[TABLE]

for any $\alpha\in[0,1]$. Let $\omega^{opt}_{t}=\omega^{opt}_{t}(\alpha)$ denote the optimal solution minimising (F) for a given $\alpha$. Using the solution (80), we find an implicit equation for $\omega^{opt}_{t}(\alpha)$:

[TABLE]

The numerical solutions for (97) at each value of $\alpha$ are finally used to compute the points on the Pareto front presented in

Fig. 2 from the main text.

Appendix G Work fluctuation-dissipation relation for discrete processes beyond weak-coupling

Our analysis in the main text assumes that at all times the system undergoes dissipative Markovian evolution, which relies on the assumption that the coupling to the bath is sufficiently weak. In this Appendix we present an alternative picture in which the system is subject to a series a fast changes in its Hamiltonian, with each quench proceeded by thermalisation with respect to the bath. However, here no restriction will be placed on the strength of coupling, and we will prove a quantum work FDR analogous to our main result Eq. (8).

Let $\{H_{\scriptscriptstyle S}^{(1)},H_{\scriptscriptstyle S}^{(2)},...H_{\scriptscriptstyle S}^{(N)}\}$ represent a particular sequence of $N$ quenches in the system Hamiltonian, where at each step we take the total Hamiltonian to be of form

[TABLE]

with $V_{\scriptscriptstyle SB}$ a time independent interaction with arbitrary coupling strength $\gamma$. The spectral decomposition of the total Hamiltonian is denoted by $H_{\scriptscriptstyle SB}^{(i)}=\sum_{n}\epsilon_{n}^{(i)}\ket{\epsilon_{n}^{(i)}}\bra{\epsilon_{n}^{(i)}}$. In addition to being initially thermal, at the end of the $(i-1)$’th quench we assume that:

• •

the system and bath equilibrates to the time-averaged state, i.e., it decoheres (see the review Gogolin and Eisert (2016) for details)

• •

the system relaxes to the reduced of a global thermal state, $\tilde{\pi}_{\scriptscriptstyle S}^{(i)}=\mbox{tr}_{B}(\pi_{\scriptscriptstyle SB}^{(i)})=\mbox{tr}_{B}(e^{-\beta H_{\scriptscriptstyle SB}^{(i)}}/Z_{\scriptscriptstyle SB}^{(i)})$. Note that $\tilde{\pi}_{\scriptscriptstyle S}^{(i)}$ deviates from a local Gibbs state unless the weak-coupling assumption is taken, namely if $\gamma^{2}\ll 1$.

It is important to realise that we make weak assumptions on the state of B, and that we only assume that thermalisation takes place at the level of S (i.e. one can imagine that SB evolve unitarily leading to thermalisation at the local level). Alternatively, one may assume that SB is put in contact via weak coupling to an external super-bath at inverse temperature $\beta$, leading to thermalisation of the full SB, but this stronger requirement is not needed for our derivation.

Crucially, each change $H_{\scriptscriptstyle S}^{(i)}\to H_{\scriptscriptstyle S}^{(i+1)}$ is sufficiently fast such that the system-bath state does not change, while the thermalisation stages occur with no work done on the system due to the fixed Hamiltonian. Since work is performed only during the quench stages, with evolution effectively unitary, the resulting work distribution can be obtained by successive iterations of the two-projective measurement protocol applied to system and bath at each stage. By treating $w=\sum^{N-1}_{i=1}w^{(i)}$ as the total work done along the N steps, which is composed of a sum of independent random variables, the work distribution is given by

[TABLE]

The first and second cumulants of work are subsequently given by

[TABLE]

Note that the above terms depend only on the system degrees of freedom since each quench applies only locally to the system Hamiltonian.

Our focus will be on slow processes, which in this context implies that the number of steps $N$ is sufficiently large. By re-expressing $\langle w\rangle$ in terms of SB as: $\langle w\rangle=\sum^{N-1}_{i=1}\text{Tr}\left((H_{\scriptscriptstyle SB}^{(i+1)}-H_{\scriptscriptstyle SB}^{(i)})\pi_{SB}^{(i)}\right)$, it is straightforward to show that the dissipated work $W_{\rm diss}=\langle w\rangle-\Delta F$ can be rewritten as Gallego et al. (2014); Perarnau-Llobet et al. (2018)

[TABLE]

where $S(\rho||\sigma)=\text{Tr}\left(\rho(\ln\rho-\ln\sigma)\right)$ is the quantum relative entropy. Note that for a density operator $\rho(t)$ that depends smoothly on some parameter $t$, the relative entropy between close states $\rho(t+\delta t)$ and $\rho(t)$ is approximated up to second order in $\delta t$ as follows Hayashi (2002):

[TABLE]

where $\mathbb{J}_{\rho}(A)$ is defined in Eq. (9). By defining $\Delta H_{\scriptscriptstyle S}^{(i)}/N=(H_{\scriptscriptstyle S}^{(i+1)}-H_{\scriptscriptstyle S}^{(i)})$ and identifying $1/N$ as a small parameter, the expansion (103) yields the following approximation for $S(\pi_{\scriptscriptstyle SB}^{(i)}||\pi_{\scriptscriptstyle SB}^{(i+1)})$ after taking the partial trace over the bath degrees of freedom:

[TABLE]

At this stage we introduce a one-parameter family of Hamiltonians $\{H_{\scriptscriptstyle S}(t)\ |\ t\in[0,1],\ H_{\scriptscriptstyle S}(t=i/N)=H_{\scriptscriptstyle S}^{(i)}\}$ and denote the corresponding system state by

[TABLE]

Then one has $\Delta H_{\scriptscriptstyle S}^{(i)}=\dot{H}_{\scriptscriptstyle S}(i/N)+\mathcal{O}(1/N)$ where $\dot{H}_{\scriptscriptstyle S}(t)=\frac{\partial}{\partial t}H_{\scriptscriptstyle S}(t)$ represents the system’s power operator. Combining (102) and (104) and taking the continuum limit gives

[TABLE]

This gives the dissipated work in the limit of many discrete steps up to first order.

We now turn to evaluating the work fluctuations $\sigma_{w}^{2}=\langle w^{2}\rangle-\langle w\rangle^{2}$ in this limit. From (101) we find in the continuum limit

[TABLE]

We are now ready to obtain the work fluctuation-dissipation relation for discrete processes. For a large number of steps $N^{2}\gg 1$, subtracting (106) from (G) gives

[TABLE]

where

[TABLE]

is the time-integrated average Wigner-Yanase-Dyson skew information, defined as in the main text by

[TABLE]

It follows from the positivity of the skew information that

[TABLE]

with equality iff $[\dot{H}_{\scriptscriptstyle S}(t),\tilde{\pi}_{\scriptscriptstyle S}(t)]=0\ \ \forall t\in[0,1]$.

Comparing (108) with Eq. (8) obtained in the Lindblad approach, we see structural similarities. Here one finds a quantum correction term stemming from the fact that the reduced state of the system does not necessarily commute with the power operator $\dot{H_{\scriptscriptstyle S}}(t)$. In the weak coupling limit where $\tilde{\pi_{\scriptscriptstyle S}}\simeq\pi_{\scriptscriptstyle S}$ we recover the same quantum correction term for Lindblad equation Eq. (10), with $N=\tau\Gamma/2$ playing the role of the ratio between the evolution and equilibration timescales. Clearly we see that the resulting work distribution has an increased spread due to additional quantum fluctuations in power, quantified by the skew information. This suggests that the assumption of weak coupling used to obtain Eq. (8) in the main text is not crucial, and that the result is more fundamentally linked to the linear expansion of the system state close to the quasi-static limit.

Appendix H Numerical verification of Eq. (8)

In order to show how the slow driving approximation affects the FDR we numerically simulate here the exact dynamics of a two level system in contact with a bosonic thermal bath with flat spectrum. The Lindblad equation we use in the eigenbasis of the system Hamiltonian is given by Breuer et al. (2002):

[TABLE]

where $r$ is the time dependent energy spacing and $P_{r}$ is simply the Planck distribution $P_{r}:=\frac{1}{e^{2\beta r}-1}$. For simplicity, we assume the proper equilibration timescale $\gamma$ to be one. Since we have to account for the explicit change of Hamiltonian we define the unitary transformation on $SU(2)$:

[TABLE]

Notice that $SU(2)$ can be parametrised by spherical angles; for this reason we write the qubit Hamiltonian in spherical coordinates as:

[TABLE]

Then in this diagonal representation the Lindbladian equation is given as:

[TABLE]

The formal solution is given by

[TABLE]

After determining the solution (116) numerically, we can compute the first two moments of work using the following exact expressions from Appendix A:

[TABLE]

These quantities are shown in Figure 3 for two different protocols as a function of the total duration $\tau$. We clearly see that our slow driving approximations (ie. $1/\tau^{2}\ll 1$) becomes valid at large values of $\tau$, thereby verifying our work FDR Eq. (8). In the first case we consider a commuting protocol where only the energy levels are changed in time, which in spherical coordinates is simply given by

[TABLE]

The Hamiltonian stays in the same energy basis throughout the protocol, and thus the quantum correction $Q_{w}=0$. In the large $\tau$ limit we see that work dissipation and fluctuations converge to the classical FDR $W_{\rm diss}=\frac{1}{2}\beta\sigma_{w}^{2}$ as expected.

For the second case we rotate the Hamiltonian energy basis using the following protocol:

[TABLE]

As a result we generate coherence during the process and hence obtain a non-zero quantum correction $Q_{w}$. In Figure 3 one sees that again the slow driving approximations for $W_{\rm diss}$ and $\frac{1}{2}\beta\sigma_{w}^{2}$ are valid in the large $\tau$ limit. In contrast to the commuting protocol, at large $\tau$ the work fluctuations exceed the dissipation, confirming the inequality Eq. (9) in the main text. The difference between the curves for $\frac{1}{2}\beta\sigma_{w}^{2}$ and $W_{\rm diss}$ represents the time-integrated skew information $Q_{w}$. Interestingly, we see that $W_{\rm diss}\leq\frac{1}{2}\beta\sigma_{w}^{2}$ remains valid at all times. It is not clear whether or not this is a general feature of the Lindblad master equation, or if this is simply a property of this particular spin-boson model.