*Removing Coherences*

To this point, we have discussed subjecting the system *S* to a cyclic process under isothermal dynamics. Now, following Ref. [9], we impose an additional condition:

$$\text{diag}\,\rho\_f = \text{diag}\,\rho\_i\tag{13}$$

where

$$\text{diag}\,\hat{\rho} = \sum\_{n} |n\rangle\langle n|\hat{\rho}|n\rangle\langle n| \tag{14}$$

is the density matrix obtained from *ρ*ˆ by setting to zero its off-diagonal elements, in the reference energy basis. In other words, we now restrict ourselves to processes that alter the system's energy coherences *m*|*ρ*<sup>ˆ</sup>|*n*-, *m* = *n*, while leaving the probabilities *n*|*ρ*<sup>ˆ</sup>|*n*- unchanged. We will refer to Equation (13), and to its classical counterpart, Equation (42), as the *isoenergetic constraint*. As in Ref. [9], our motivation for imposing this condition is to isolate and accentuate the thermodynamic implications of quantum energy coherences. From Equation (13), it follows that Tr[*H*<sup>ˆ</sup> 0*ρ*<sup>ˆ</sup>*i*] = Tr[*H*<sup>ˆ</sup> 0*ρ*<sup>ˆ</sup>*f* ], i.e.,

$$\mathcal{U}\_i^{\emptyset} = \mathcal{U}\_f^{\emptyset} \tag{15}$$

which in turn implies that the generalized second law, Equation (5), becomes

$$\mathcal{W}^q \le \mathcal{B}^{-1} \left( \mathcal{S}\_f^q - \mathcal{S}\_i^q \right). \tag{16}$$

This bound relates the maximum extractable work to the change in the system's entropy. The thermodynamic interpretation is clear: since the system's energy undergoes no net change (Equation (15)), the only way to extract work is to withdraw energy from the bath, causing the entropy of the bath to decrease by an amount *β*W*q*. This decrease in the bath's entropy must be compensated, or over-compensated, by an increase in the entropy of the system, as reflected by Equation (16).

We are now in a position to investigate the maximum amount of work that can be extracted from energy coherences. For a given reference Hamiltonian *H* ˆ 0 and initial state *ρ* ˆ *i*, let W*q* denote the maximum extracted work, over all protocols *H*ˆ (*t*) that begin and end in *H* ˆ 0, subject to the isoenergetic constraint (13). Since the right side of Equation (16) is a function of *ρ*ˆ*i* and *ρ*ˆ*f* , we can place a bound on W*q* by maximizing that function with respect to *ρ*ˆ*f* :

$$\mathcal{W}^{q\star} \le \beta^{-1} \max\_{\hat{\rho}\_f \mid \text{diag}\,\hat{\rho}\_f = \text{diag}\,\hat{\rho}\_i} \left[ \mathcal{S}^q(\hat{\rho}\_f) - \mathcal{S}^q(\hat{\rho}\_i) \right] \tag{17}$$

For fixed diagonal elements of a density matrix *ρ*ˆ, the value of S*q* = −Tr*ρ*<sup>ˆ</sup> ln *ρ*ˆ is maximized when the off-diagonal elements are all zero. We therefore obtain

$$\mathcal{W}^{q\*} \le \mathcal{S}^{-1}[\mathcal{S}^q(\text{diag}\,\rho\_i) - \mathcal{S}^q(\rho\_i)].\tag{18}$$

This result does not ye<sup>t</sup> tell us whether the bound can be saturated, that is, whether there exist protocols for extracting this amount of work. Rather, it states that under no circumstances can we extract more than this much work, in a cyclic, isothermal process satisfying Equation (13). Moreover, if a protocol for saturating this bound exists, then that protocol will result in the system ending in the state diag *ρ*ˆ*i* at *t* = *τ*. In other words, the saturating protocol (if it exists) removes all energy coherences from the system's initial state, and effectively converts these coherences into extracted work.

In fact, protocols for saturating the bound given by Equation (18) do exist [9,28]. A simple example is given by:

$$\hat{H}(t) = \begin{cases} \hat{H}\_0 & t \le 0 \\ -\beta^{-1} \ln[(1-\lambda)\beta\_i + \lambda(\text{diag}\,\hat{\rho}\_i)] & 0 < t < \tau \\ \hat{H}\_0 & \tau \le t \end{cases} \tag{19}$$

where *λ* ≡ *t*/*τ* varies from 0 to 1 during the process, and *τ* is taken to be sufficiently large that the process is quasistatic. This protocol can be understood as follows. At the start of the process, there is a sudden change, or *quench*, in the system's Hamiltonian, from *H*ˆ 0 at *t* = 0 to −*β*−<sup>1</sup> ln *ρ*ˆ*i* at *t* = 0+. Thus, at *t* = 0<sup>+</sup> the system's state *ρ*ˆ*i* is in equilibrium with respect to the immediate post-quench Hamiltonian. (The term "quench" is often used in situations in which the system is in equilibrium before the quench, and out of equilibrium after it. Thus, the first step of this protocol (19) might be viewed as an *anti*-quench.) From *t* = 0<sup>+</sup> to *<sup>τ</sup>*<sup>−</sup>, the Hamiltonian is varied quasistatically from −*β*−<sup>1</sup> ln *ρ*ˆ*i* to −*β*−<sup>1</sup> ln(diag *ρ*<sup>ˆ</sup>*i*), and the system is dragged through the corresponding sequences of equilibrium states, from *ρ* ˆ *i* to diag *ρ*<sup>ˆ</sup>*i*—see comments after Equation (4). At *t* = *τ*, a second quench abruptly returns the Hamiltonian to *H* ˆ 0, completing the cycle. The evolution of the system's state is thus given by

$$\rho(t) = \begin{cases} \rho\_i & t = 0 \\ (1 - \lambda)\rho\_i + \lambda (\text{diag}\,\rho\_i) & 0 < t < \tau \\ \text{diag}\,\rho\_i & \tau \le t. \end{cases} \tag{20}$$

We show in Appendix A that the work extracted during this process is given by the right side of Equation (18), that is the bound is saturated. Hence, under the isoenergetic constraint (13), work extraction is optimized by removing all coherences from the system's state, and the value of this optimized work is:

$$\mathcal{W}^{q\*} = \mathcal{S}^{-1} |\mathcal{S}^q(\text{diag}\,\rho\_i) - \mathcal{S}^q(\rho\_i)|. \tag{21}$$

This result is equivalent to Equation (1) of Kammerlander and Anders [9].

#### **3. Classical Setup and Notation**

Now, imagine a classical system with *N* degrees of freedom and phase space variables

$$
\Gamma = (\mathbf{x}\_1, \dots, \mathbf{x}\_N, p\_1, \dots, p\_N). \tag{22}
$$

Adopting (as in the quantum case) an ensemble perspective, let the system's state at time *t* be described by a phase space density *ρ*(<sup>Γ</sup>, *t*). We will consider a thermodynamic process in which the system begins in a state *ρ*(<sup>Γ</sup>, 0) = *ρi*(Γ), then evolves from *t* = 0 to *τ* as its Hamiltonian is varied according to a cyclic protocol *<sup>H</sup>*(<sup>Γ</sup>, *t*), with

$$H(\Gamma, 0) = H(\Gamma, \tau) = Ho(\Gamma) \tag{23}$$

where *<sup>H</sup>*0(Γ) specifies a reference Hamiltonian; see Figure 2. We assume that no observables commute with *<sup>H</sup>*0(Γ) under the Poisson bracket, except those that are functions of *H*0:

$$\{K, H\_0\} = 0 \text{ iff } K(\Gamma) = k(H\_0(\Gamma)) \tag{24}$$

for some function *k*(·) (compare with Equation (2)). This assumption implies that energy is the only non-trivially conserved quantity along all trajectories Γ(*t*) obeying Hamiltonian dynamics *d*Γ/*dt* = {<sup>Γ</sup>, *<sup>H</sup>*0}. (This conclusion follows from the identity (*d*/*dt*)*A*(Γ(*t*)) = {*<sup>A</sup>*, *<sup>H</sup>*0}, which applies to any observable *A*(Γ) and any trajectory Γ(*t*) obeying *d*Γ/*dt* = {<sup>Γ</sup>, *<sup>H</sup>*0}). This is a necessary but not sufficient condition for the dynamics to be *ergodic* on constant-energy surfaces in phase space—an assumption often made in statistical physics. (Roughly speaking, ergodicity means that a generic Hamiltonian trajectory of energy *E* visits all regions of the surface *H*0 = *E*, given sufficient time.) For our purposes, we do not need the assumption of ergodicity, only the weaker assumption given by Equation (24).

**Figure 2.** Schematic illustration of the classical process. The system begins in state *ρi*(Γ), then evolves in contact with a thermal bath to a final state *<sup>ρ</sup>f*(Γ) as the Hamiltonian is driven through a cycle from *<sup>H</sup>*(<sup>Γ</sup>, 0) = *<sup>H</sup>*0(Γ) to *<sup>H</sup>*(<sup>Γ</sup>, *τ*) = *<sup>H</sup>*0(Γ). The constraint diag *ρi*(Γ) = diag *<sup>ρ</sup>f*(Γ) indicates that the initial and final energy distributions are identical, while inhomogeneities may differ.

If the system were thermally isolated, then its state *ρ*(<sup>Γ</sup>, *t*) would evolve under the Liouville equation, *∂ρ*/*∂t* = {*<sup>H</sup>*, *ρ*}. However, we assume that the system is in contact with a thermal bath as it undergoes the cyclic process, hence its evolution follows classical isothermal dynamics, rather than Hamiltonian dynamics. As in the quantum case, we will not specify the equations of motion that describe the isothermal dynamics, but we will make the following assumptions.

(1) If the system's Hamiltonian is held fixed, then the isothermal dynamics drive the system to the equilibrium state

$$
\pi(\Gamma) = \frac{1}{Z^{\varepsilon}} e^{-\beta H(\Gamma)} \tag{25}
$$

with partition function and free energy

$$Z^{\mathbb{C}}[H(\Gamma)] = \int d\Gamma \, e^{-\beta H(\Gamma)} \quad , \quad \mathcal{F}^{\text{c.eq}}[H(\Gamma)] = -\beta^{-1} \ln \left( Z^{\mathbb{C}} / h^N \right) . \tag{26}$$

Here, *h* is a constant with dimensions of action that ensures the argumen<sup>t</sup> of the logarithm is dimensionless. We choose *h* to coincide with Planck's constant as this will facilitate comparisons of quantum and classical work extraction in Section 4. We assume this relaxation takes place over a finite timescale *τ*rel. As a consequence, if *<sup>H</sup>*(<sup>Γ</sup>, *t*) is varied quasistatically, then the system's state follows the instantaneous equilibrium state, *ρ*(<sup>Γ</sup>, *t*) = *<sup>π</sup>*(<sup>Γ</sup>, *t*).

(2) When the system evolves over a time interval 0 ≤ *t* ≤ *τ* under isothermal dynamics and a time-dependent Hamiltonian *<sup>H</sup>*(<sup>Γ</sup>, *t*), it obeys a generalized second law

$$
\mathcal{W}^{\mathbb{C}} \le -\Delta \mathcal{F}^{\mathbb{C}} \tag{27}
$$

with

$$\mathcal{W}^{\mathbb{C}}[\rho(\Gamma, t), H(\Gamma, t)] = -\int\_0^{\tau} \left( \int \frac{\partial H}{\partial t} \rho d\Gamma \right) dt \tag{28}$$

$$\left[\mathcal{F}^{\mathbb{C}}[\rho(\Gamma), H(\Gamma)]\right] = \mathcal{U}^{\mathbb{C}} - \mathcal{S}^{\mathbb{C}}/\beta \tag{29}$$

$$\mathcal{U}^c[\rho(\Gamma), H(\Gamma)] = \int H(\Gamma)\rho(\Gamma)d\Gamma \tag{30}$$

$$\mathcal{S}^{\mathbb{C}}[\rho(\Gamma)] = -\int \rho(\Gamma) \ln[h^N \rho(\Gamma)] d\Gamma. \tag{31}$$

Unlike the quantum von Neumann entropy (9) which is always non-negative, the classical Shannon differential (or continuous) entropy (31) can become arbitrarily negative for probability distributions that are highly concentrated in phase space, as we will see in Section 4.1. (For brevity, we will henceforth refer to the Shannon differential entropy simply as the Shannon entropy.) As with the quantum bound (Equation (5)), Equation (27) is not restricted to transitions between equilibrium states, and has been derived under a variety of modeling approaches, see, e.g., Refs. [26–30,42–44].

The classical non-equilibrium free energy (29) can be rewritten as

$$\mathcal{F}^{\mathbb{C}}[\rho, H] = \mathcal{F}^{\mathbb{C}, \text{eq}}[H] + \beta^{-1} D[\rho|\pi] \tag{32}$$

where

$$D[\rho\_1|\rho\_2] = \int d\Gamma \,\rho\_1(\Gamma) \ln \frac{\rho\_1(\Gamma)}{\rho\_2(\Gamma)} \ge 0 \tag{33}$$

is the classical relative entropy or Kullback–Leibler divergence. Thus, for a cyclic process, Equation (27) becomes

$$\mathcal{W}^c \le \beta^{-1} \left( D[\rho\_i | \pi\_0] - D[\rho\_f | \pi\_0] \right) \tag{34}$$

where *ρi*(Γ) = *ρ*(<sup>Γ</sup>, 0) and *<sup>ρ</sup>f*(Γ) = *ρ*(<sup>Γ</sup>, *τ*) are the system's initial and final states, and *<sup>π</sup>*0(Γ) is the equilibrium state for the reference Hamiltonian. As in the quantum case (Equation (12)), the right side of Equation (34) provides a measure of the degree to which the process brings the system closer to equilibrium.

#### *3.1. Energy-Shell Inhomogeneities*

The evident similarity between the quantum framework for cyclic isothermal processes described by Equations (1)–(12) and the classical framework of Equations (23)–(34) motivates us to seek a classical analogue of the statement that quantum energy coherences represent a thermodynamic resource. As a step in this direction, we note that in the quantum case, density matrices that are stationary under the unitary evolution generated by *H*ˆ 0 are exactly those that lack energy coherences in the eigenbasis of *H* ˆ 0:

$$\frac{d\boldsymbol{\hat{\rho}}}{dt} = \frac{1}{i\hbar} [\hat{H}\_0, \boldsymbol{\hat{\rho}}] = 0 \qquad \text{iff} \qquad \boldsymbol{\hat{\rho}} = \text{diag}\, \boldsymbol{\hat{\rho}}.\tag{35}$$

In the classical case, phase space densities that are stationary under the Hamiltonian dynamics generated by *<sup>H</sup>*0(Γ) are exactly those that are functions of *<sup>H</sup>*0(Γ):

$$\frac{\partial \rho}{\partial t} = \{H\_0, \rho\} = 0 \qquad \text{iff} \qquad \rho(\Gamma) = k(H\_0(\mathbf{x})) \tag{36}$$

for some function *k*(·). (Equations (2) and (24) are needed for the "only if" parts of Equations (35) and (36).) These observations sugges<sup>t</sup> that we ought to view phase space distributions of the form *ρ* = *k*(*<sup>H</sup>*0) as analogues of density matrices that are diagonal in the eigenbasis of *H* ˆ 0.

To pursue this idea, let *η*(*E*) denote the distribution of energies associated with a phase space density *ρ*(Γ):

$$\eta(E) = \int d\Gamma \,\rho(\Gamma) \,\delta[E - H\_0(\Gamma)].\tag{37}$$

In addition, let *<sup>ω</sup>E*(Γ) denote the classical microcanonical density of energy *E*:

$$\omega\_E(\Gamma) = \lim\_{\Delta E \to 0} \frac{I\_{[E, E + \Delta E]}(H\_0(\Gamma))}{\int d\Gamma' I\_{[E, E + \Delta E]}(H\_0(\Gamma'))} = \frac{\delta[E - H\_0(\Gamma)]}{\Omega(E)}\tag{38}$$

where *<sup>I</sup>*[*<sup>E</sup>*,*E*+Δ*E*](·) is the indicator function over the interval [*<sup>E</sup>*, *E* + Δ*E*] (that is, *<sup>I</sup>*[*<sup>a</sup>*,*<sup>b</sup>*](*x*) = 1 when *x* ∈ [*a*, *b*], otherwise *<sup>I</sup>*[*<sup>a</sup>*,*<sup>b</sup>*](*x*) = 0), and

$$
\Omega(E) = \int d\Gamma \,\delta[E - H\_0(\Gamma)]\tag{39}
$$

is the classical density of states. The microcanonical density *<sup>ω</sup>E*(Γ) is singular, uniformly distributed over the *energy shell E* (the level set *H*0 = *E*), and zero elsewhere. Here, "uniformly distributed" is defined by Equation (38): as Δ*E* approaches zero, the phase space density remains uniform, with respect to the Liouville measure *dNx d<sup>N</sup> p*, in the region between shells *E* and *E* + Δ*E*, and zero elsewhere.

Using Equations (37)–(39), a phase space density of the form *ρ* = *k*(*<sup>H</sup>*0) can be written as 

$$
\rho(\Gamma) = \int dE \,\eta(E) \,\omega\_E(\Gamma)\tag{40}
$$

with *η*(*E*) = *k*(*E*)Ω(*E*). Such a density is a statistical mixture of microcanonical ensembles (just as a diagonal density matrix is a mixture of energy eigenstates: *ρ*ˆ = diag *ρ*ˆ = ∑*n pn*|*nn*|), hence *ρ*(Γ) is uniform, or *homogeneous*, over any specific energy shell *E*, while its value differs from one shell to another. By contrast, a phase space density that is not of the form *ρ* = *k*(*<sup>H</sup>*0) is *inhomogeneous* on energy shells: there exist points Γ and Γ such that *<sup>H</sup>*0(Γ) = *<sup>H</sup>*0(Γ) but *ρ*(Γ) = *ρ*(Γ).

We will henceforth use the terms homogeneous/inhomogeneous to distinguish between phase space densities that can/cannot be written as *ρ* = *k*(*<sup>H</sup>*0). For instance, the equilibrium distribution *<sup>π</sup>*0(Γ) ∝ exp −*βH*0(Γ) is a homogeneous density. By the stationarity argumen<sup>t</sup> given above (Equations (35) and (36)), homogeneous phase space densities

will be viewed as classical counterparts of diagonal density matrices, and inhomogeneous densities as counterparts of quantum states with energy coherences. In other words, for our purposes *the counterparts of quantum energy coherences are classical energy-shell inhomogeneities*. Weintroducethenotation

> diag *ρ*(Γ) = *dE η*(*E*) *<sup>ω</sup>E*(Γ) (41)

with *η*(*E*) given by Equation (37), to denote the phase space density obtained by "homogenizing" *ρ*(Γ). That is, diag *ρ* is the homogeneous density that has the same energy distribution as *ρ*.

#### *3.2. Removing Inhomogeneities*

Let us now focus our attention on classical, cyclic isothermal processes that satisfy the isoenergetic constraint (compare with Equation (13)):

$$\text{diag}\,\rho\_f(\Gamma) = \text{diag}\,\rho\_l(\Gamma) \tag{42}$$

where *ρi*, *f* denote the system's initial and final states. Such processes leave the energy distribution undisturbed, *ηi*(*E*) = *<sup>η</sup>f*(*E*), while allowing energy-shell inhomogeneities to change. Equation (42) implies

$$\mathcal{U}\_i^\varepsilon = \mathcal{U}\_f^\varepsilon \tag{43}$$

hence, Equation (27) becomes

$$\mathcal{W}^{\varepsilon} \le \beta^{-1} \left( \mathcal{S}\_f^{\varepsilon} - \mathcal{S}\_i^{\varepsilon} \right). \tag{44}$$

Let W*c* denote the maximum amount of work that can be extracted, over all conceivable cyclic protocols, for a given reference Hamiltonian *<sup>H</sup>*0(Γ) and initial state *ρi*(Γ). Equation (44) implies

$$\mathcal{W}^{\varepsilon \star} \quad \le \quad \beta^{-1} \max\_{\rho\_f \mid \text{diag}\,\rho\_f = \text{diag}\,\rho\_i} \left( \mathcal{S}^{\varepsilon}[\rho\_f] - \mathcal{S}^{\varepsilon}[\rho\_i] \right)$$

$$= \quad \beta^{-1} (\mathcal{S}^{\varepsilon}[\text{diag}\,\rho\_i] - \mathcal{S}^{\varepsilon}[\rho\_i]) \tag{45}$$

since, among all states with a given energy distribution, the Shannon entropy is maximized by the homogeneous state (This follows from the fact that Shannon entropy increases under coarse-graining, which in turn is a consequence of Jensen's inequality, ln *x*- ≤ ln*x*-).

Similarly to the quantum case (Equation (19)), the bound given by Equation (45) is saturated [27,28] by the protocol

$$H(\Gamma, t) = \begin{cases} H\_0(\Gamma) & t \le 0 \\ -\beta^{-1} \ln[(1 - \lambda)\rho\_i(\Gamma) + \lambda \operatorname{diag} \rho\_i(\Gamma)] & 0 < t < \tau \\ H\_0(\Gamma) & \tau \le t \end{cases} \tag{46}$$

with *λ* = *t*/*<sup>τ</sup>*, and *τ* sufficiently long that the process is effectively quasistatic. The protocol begins with a classical quench at *t* = 0. Immediately after this quench, the system's state *ρi*(Γ) is in equilibrium with its instantaneous Hamiltonian, *<sup>H</sup>*(<sup>Γ</sup>, 0+) = −*β*−<sup>1</sup> ln *ρi*(Γ). During the interval *t* ∈ (0, *<sup>τ</sup>*), the quasistatic switching of the Hamiltonian drags the system through a sequence of equilibrium states from *ρi* to *ρf* = diag *ρi*, and at *t* = *τ* the cyclic process is completed by suddenly returning the Hamiltonian to *H*0. The evolution of the system's state *ρ*(<sup>Γ</sup>, *t*) is entirely analogous to that given by Equation (20). Summing over the work extracted during the initial quench, the quasistatic driving, and the final quench, we find (see Appendix A) that the total extracted work is

$$\mathcal{W}^{\varepsilon \star} = \mathcal{S}^{-1}(\mathcal{S}^{\varepsilon}[\text{diag}\,\rho\_{i}] - \mathcal{S}^{\varepsilon}[\rho\_{i}]),\tag{47}$$

i.e., the bound in Equation (45) is saturated. By Equation (44), any protocol satisfying Equation (42) that would bring the system to a final state *ρf* = diag *ρi* would necessarily result in less work extracted.

#### **4. Quantum–Classical Comparison**

We have seen that the maximum work extracted in the quantum case, subject to the isoenergetic constraint, diag *ρ*ˆ*i* = diag *ρ*ˆ*f* , is achieved by quasistatically removing all energy coherences from the system's initial state: *ρ*ˆ*i* → diag *ρ*<sup>ˆ</sup>*i*. Similarly, the maximum work extracted in the classical case is achieved by quasistatically removing all energy-shell inhomogeneities. The optimized work values W*q* and W*c* are given by Equations (21) and (47). The close similarity between these results supports our view that classical energy-shell inhomogeneities are thermodynamic counterparts of quantum energy coherences. Both are resources that can be leveraged to extract work.

While the expressions for W*q* and W*c* are nearly identical, it still remains to compare them quantitatively. Ideally, we would like to compare the values of W*q* and W*c* for a given quantum reference Hamiltonian *H*ˆ 0 and initial state *ρ*<sup>ˆ</sup>*i*, and appropriately defined classical counterparts *<sup>H</sup>*0(Γ) and *ρi*(Γ). To this end, throughout this section and the next we assume that *H*ˆ 0 is a function of position and momentum operators (*x*ˆ1, ··· *<sup>x</sup>*<sup>ˆ</sup>*N*) and (*p*ˆ1, ··· *p*<sup>ˆ</sup>*N*), and we further assume that *H*ˆ 0 has a well-defined counterpart *<sup>H</sup>*0(Γ). This condition is satisfied, for instance, by Hamiltonians of the kinetic-plus-potential form *H*ˆ 0 = *<sup>K</sup>*(*p*ˆ1, ··· *p*<sup>ˆ</sup>*N*) + *<sup>V</sup>*(*x*ˆ1, ··· *<sup>x</sup>*<sup>ˆ</sup>*N*), for which *<sup>H</sup>*0(Γ) is obtained by replacing momentum and position operators with classical momentum and position variables.

Identifying a correspondence between quantum and classical states *ρ*ˆ and *ρ*(Γ) is trickier. Common approaches that map density operators into phase space distributions [45,46] suffer from undesirable properties. For instance, neither the Wigner [47] nor Husimi [48] function representation of the quantum thermal state corresponds to the classical thermal phase space distribution. Additionally, the Wigner function in general can become negative while the Husimi function depends on the choice of coherent states.

To circumvent such issues, we will compare quantum and classical *energy distributions* rather than individual states. Instead of focusing on the maximum work that can be extracted from a particular initial state, we will consider the maximum work that can be extracted given a particular initial energy distribution. We begin by defining *energy equivalence classes* in Section 4.1, then in Sections 4.2 and 4.3 we compare maximum work values for corresponding quantum and classical energy equivalence classes.

#### *4.1. Energy Equivalence Classes*

We define a quantum energy equivalence class to consist of all states *ρ*ˆ that share a particular energy distribution, that is, a particular set of diagonal density matrix elements, with respect to *H*ˆ 0. An example is the thermal energy equivalence class given by

$$
\Pi^{\emptyset} = \{ \not\phi \, | \, \text{diag} \, \not\rho = \hat{\pi}\_0 \} \tag{48}
$$

where *π*ˆ 0 = exp(−*βH*<sup>ˆ</sup> 0)/*Zq* is the thermal equilibrium state. In addition to the state *π*ˆ 0, the set Π*q* includes exotic non-equilibrium states with significant energy coherences such as the pure state |*<sup>π</sup>*0- *<sup>π</sup>*0|, where

$$|\pi\_0\rangle = \sum\_n \sqrt{\frac{\varepsilon^{-\beta\varepsilon\_n}}{Z^q}} |n\rangle. \tag{49}$$

Examples of this state arise in quantum optics [49,50].

More generally (that is, not restricting ourselves to the thermal energy equivalence class, Equation (48)), every quantum state *ρ*ˆ belongs to a unique energy equivalence class Σ*q*

defined by the diagonal elements of *ρ*ˆ in the *H* ˆ 0 basis. Within this class, the von Neumann entropy is maximized by the state *σ*ˆ ≡ diag *ρ*ˆ = ∑*n pn*|*nn*|:

$$\max\_{\boldsymbol{\beta}\in\Sigma^{\mathcal{G}}} \mathcal{S}^{\boldsymbol{q}}(\boldsymbol{\beta}) = \mathcal{S}^{\boldsymbol{q}}(\boldsymbol{\beta}) = -\sum\_{\boldsymbol{n}} p\_{\boldsymbol{n}} \ln p\_{\boldsymbol{n}} \geq 0 \tag{50a}$$

where *pn* = *n*|*ρ*<sup>ˆ</sup>|*n*-. The von Neumann entropy is minimized within Σ*q* by pure states such as |*ψψ*|, where |*ψ*- = ∑*n* √*pn*|*n*-, and for these states the entropy vanishes:

$$\min\_{\boldsymbol{\uprho}\in\Sigma^{\mathcal{G}}} \mathcal{S}^{\boldsymbol{\uprho}}(\boldsymbol{\uprho}) = \mathcal{S}^{\boldsymbol{\uprho}}(|\boldsymbol{\uppsi}\rangle\langle\boldsymbol{\uppsi}|) = 0. \tag{50b}$$

A classical energy equivalence class contains all phase space distributions *ρ*(Γ) with a given energy distribution *η*(*E*). An example is the thermal energy equivalence class

$$\Pi^{\mathfrak{c}} = \{ \rho(\Gamma) \, | \, \text{diag}\, \rho(\Gamma) = \pi\_0(\Gamma) \}\tag{51}$$

where *<sup>π</sup>*0(Γ) = exp[−*βH*0(Γ)]/*Zc*. While the state *<sup>π</sup>*0(Γ) is homogeneous, the class Π*c* contains states with substantial energy-shell inhomogeneities. For instance, if the system is a one-dimensional harmonic oscillator, the thermal equivalence class Π*c* includes the state

$$\rho(E,T) = \underbrace{\left(\beta e^{-\beta E}\right)}\_{\eta\_{\pi\_0}(E)} \underbrace{\left(\omega \frac{e^{\delta \cos(\omega T)}}{2\pi I\_0(\delta)}\right)}\_{\xi(T)}\tag{52}$$

where *E* and *T* are the canonical energy and tempus (angle-like) coordinates [51] defined by *x* = √2*E*/*mω*<sup>2</sup> cos(*ω<sup>T</sup>*) and *p* = √2*mE* sin(*ω<sup>T</sup>*), with *ωT* ∈ (−*π*, <sup>+</sup>*π*]; *δ* is a non-negative parameter; and *I*0 is the modified Bessel function of order zero. For this example, it is convenient to use (*<sup>E</sup>*, *T*) rather than (*<sup>x</sup>*, *p*) to identify a point in classical phase space. *ζ*(*T*) is the von Mises distribution [52], an analogue of a Gaussian distribution for an angular coordinate. In Equation (52), the mean of *ζ*(*T*) is zero and its variance is controlled by *δ*. For *δ* = 0, *ρ*(*<sup>E</sup>*, *T*) reduces to the canonical distribution, which is homogeneous over every energy shell. With increasing *δ*, the distribution becomes more and more concentrated on the positive *x*-axis of phase space (where *T* = 0) and as a result its Shannon entropy S*c*[*ρ*] decreases, with no lower bound. Specifically, for large *δ*, we have

$$\mathcal{S}^{\varepsilon}[\rho(E,T)] \approx -\ln\left(h\beta\omega\sqrt{\frac{\delta}{2\pi}}\right) + \frac{1}{2} \quad , \quad \delta \gg 1. \tag{53}$$

Every classical state *ρ*(Γ) belongs to a unique energy equivalence class Σ*<sup>c</sup>*, defined by its energy distribution *η*(*E*) (Equation (37)). Within this class, the Shannon entropy is maximized by the diagonal state *σ*(Γ) = diag *ρ*(Γ) = *η*(*<sup>H</sup>*0)/Ω(*<sup>H</sup>*0), but there is no lower bound on the minimum entropy, as the phase space distribution can be concentrated to an arbitrary degree without affecting the energy distribution:

$$\max\_{\rho \in \Sigma^{\mathbb{C}}} \mathcal{S}^{\varepsilon}[\rho(\Gamma)] = \mathcal{S}^{\varepsilon}[\sigma(\Gamma)] = -\int dE \,\eta(E) \ln \left[ h^{N} \frac{\eta(E)}{\Omega(E)} \right] \tag{54}$$

$$\min\_{\rho \in \Sigma^{\varepsilon}} \mathcal{S}^{\varepsilon}[\rho(\Gamma)] = -\infty. \tag{55}$$

These extrema are illustrated by the values *δ* = 0 and *δ* → ∞ in the example in the previous paragraph.

To take another illustrative example—which will prove useful in the next section— consider an ideal gas of *n* particles inside a three-dimensional cubic box of volume *V* = *L*3, oriented parallel to the *x*-, *y*-, and *z*-axes, with one corner at the origin—see Figure 3. A point in phase space is given by Γ = (**<sup>r</sup>**1 ···**r***n*; **p**1 ··· **<sup>p</sup>***n*). For 0 < *α* ≤ 1, let *ρα*(Γ) denote

the distribution for which the momenta **p***k* are sampled from the Maxwellian distribution at temperature *β*−1, and the positions **<sup>r</sup>***k* are sampled uniformly within the region defined by 0 < *x*, *y* < *L* and 0 < *z* < *αL*. This distribution belongs to the thermal energy equivalence class Π*<sup>c</sup>*, and *ρα*=<sup>1</sup>(Γ) is exactly the (homogeneous) thermal distribution, whereas *ρα*<<sup>1</sup>(Γ) is an inhomogeneous, non-equilibrium distribution, in which the gas is entirely located within a fraction *α* of the volume of the box. For arbitrary *α* ∈ (0, 1], we have

$$\mathcal{S}^{\varepsilon}[\rho\_a(\Gamma)] = n \ln \left(\frac{\kappa V}{\lambda\_{\text{th}}^3}\right) + \frac{3n}{2} \tag{56}$$

where *λ*th = )*βh*2/2*πm* is the thermal de Broglie wavelength. The value S*c*[*ρα*] is maximized at *α* = 1, that is, for the homogeneous state, and it has no lower bound as *α* → 0.

**Figure 3.** An ideal gas inside a box of volume *L*3. The value of *α* ∈ (0, 1] parametrizes a family of energy equivalence classes, with *α* = 1 corresponding to thermal class Π*<sup>c</sup>*. See text for details.

In both of the above examples, by "squeezing" *ρ*(Γ) into an arbitrarily small region of phase space (*δ* → <sup>∞</sup>, *α* → 0) we obtain a distribution with arbitrarily large, negative entropy.

#### *4.2. An Unfair Comparison*

We now determine the maximum amount of work that can be extracted in a cyclic isoenergetic process where all states in the quantum equivalence class Σ*q* are considered. Using Equations (21) and (50b), we have

$$\begin{split} \max\_{\boldsymbol{\hat{\rho}}\_{i} \in \Sigma^{q}} \mathcal{W}^{q\*} (\boldsymbol{\hat{\rho}}\_{i}) &= \boldsymbol{\beta}^{-1} \max\_{\boldsymbol{\hat{\rho}}\_{i} \in \Sigma^{q}} [\mathcal{S}^{q} (\text{diag} \, \boldsymbol{\hat{\rho}}\_{i}) - \mathcal{S}^{q} (\boldsymbol{\hat{\rho}}\_{i})] \\ &= \boldsymbol{\beta}^{-1} \Big[ \mathcal{S}^{q} (\boldsymbol{\mathcal{\sigma}}) - \min\_{\boldsymbol{\hat{\rho}}\_{i} \in \Sigma^{q}} \mathcal{S}^{q} (\boldsymbol{\hat{\rho}}\_{i}) \Big] = \boldsymbol{\beta}^{-1} \mathcal{S}^{q} (\boldsymbol{\mathcal{\sigma}}) \end{split} \tag{57}$$

where *σ*ˆ = diag *ρ*ˆ*i* is the unique diagonal state belonging to Σ*q*. The minimal value of S*q*(*ρ*<sup>ˆ</sup>*i*) on the second line is achieved for any pure state |*ψψ*| ∈ Σ*q*, an example of which can always be constructed using the same argumen<sup>t</sup> as in Equation (50b). Hence, the maximum work is obtained by starting in a pure state, then quasistatically removing the coherences (e.g., following the protocol given by Equation (19)) so as to end in the diagonal state *σ*ˆ. This result has a simple interpretation in terms of the bound W*q* ≤ *<sup>β</sup>*−<sup>1</sup>S*qf* − S*qi* (see Equation (16)): we maximize the extracted work by starting in a state with the lowest entropy and ending in the state of highest entropy, within Σ*q*. By Equation (50) these are, respectively, any pure state and the unique diagonal state in Σ*q*. Equivalently (since U*qf* = U*qi* by Equation (13)), the maximum extracted work is obtained when starting in the state of highest free energy and ending in the state of lowest free energy. We emphasize that, here, free energy and entropy are defined by Equations (7) and (9), which apply to generic (not necessarily equilibrium) quantum states *ρ*ˆ.

The analogous classical calculation, using Equations (47) and (55), gives

$$\begin{split} \max\_{\rho\_{i}\in\Sigma^{\varepsilon}} \mathcal{W}^{\varepsilon\star}[\rho\_{i}(\Gamma)] &= \beta^{-1} \max\_{\rho\_{i}\in\Sigma^{\varepsilon}} (\mathcal{S}^{\varepsilon}[\text{diag}\,\rho\_{i}(\Gamma)] - \mathcal{S}^{\varepsilon}[\rho\_{i}(\Gamma)]) \\ &= \beta^{-1} \Big( \mathcal{S}^{\varepsilon}[\sigma(\Gamma)] - \min\_{\rho\_{i}\in\Sigma^{\varepsilon}} \mathcal{S}^{\varepsilon}[\rho\_{i}(\Gamma)] \Big) = +\infty \end{split} \tag{58}$$

where *σ*(Γ) = diag *ρi*(Γ). In other words, for a given classical energy distribution, there is no upper bound on the amount of work that can be extracted, as there is no lower bound on the entropy of the initial state. By "squeezing" a given phase space distribution *within* each energy shell, without altering the distribution of probability *among* energy shells, we can construct a distribution *ρi*(Γ) that is compressed within an arbitrarily small volume of phase space, hence we can make the value of S*c*[*ρi*(Γ)] arbitrarily small. This idea is illustrated by Equation (52) for the harmonic oscillator example of the previous section: as *δ* → <sup>∞</sup>, the von Mises distribution *ζ*(*T*) becomes ever more concentrated around *T* = 0, and the entropy of the distribution becomes arbitrarily large and negative.

The example of the ideal gas discussed at the end of Section 4.1 provides further intuition for Equation (58). For that example, consider the thermal equivalence class Π*<sup>c</sup>*, and imagine an initial inhomogeneous distribution *ρi*(Γ) = *ρα*(Γ) at *t* = 0, with *α* < 1, that is, with all gas particles initially located in the region 0 < *z* < *αL*. To maximize the extracted work, we first suddenly insert a partition at the location *z* = *<sup>α</sup>L*, and then quasistatically move this partition to the location *z* = *L*, while the system remains in contact with a thermal bath at temperature *β*−1. The process ends with the system in the homogeneous, thermal state *<sup>ρ</sup>f*(Γ) = *ρα*=<sup>1</sup>(Γ). The total work extracted during this process of removing inhomogeneities is

$$\mathcal{W}^{\mathbb{C}} = n\beta^{-1}\ln\frac{1}{\alpha} > 0\tag{59}$$

which follows from a well-known expression for the reversible isothermal expansion of an ideal gas: W = *nβ*−<sup>1</sup> ln(*Vf* /*Vi*). It is easy to see why there is no upper bound on the extractable work: at *t* = 0+, just after the insertion of the partition, the gas is an equilibrium state, confined within a volume *<sup>α</sup>V*, with free energy F*<sup>c</sup>*(*t* = 0+) = −*<sup>n</sup>β*−<sup>1</sup> ln(*αV*/*λ*3th). The smaller the value of *α*, the larger the initial free energy and therefore the greater the amount of work that can be extracted through reversible, isothermal expansion. In this idealized example, we can begin with an arbitrarily dense initial state, i.e., arbitrarily small *α* > 0.

In both the quantum and classical cases, the extracted work is maximized by evolving quasistatically from the state of lowest entropy to the state of highest entropy, within the equivalence class Σ*q* or Σ*<sup>c</sup>*. Thus, there appears to be an inherent quantum thermodynamic *disadvantage*, since S*q* is bounded from below by 0, while S*c* is unbounded from below.

The comparison, however, is unfair. Quantum mechanics obeys the Heisenberg uncertainty principle, a loose semiclassical interpretation of which states that every quantum state occupies a cell of volume *h<sup>N</sup>* in phase space. If we view classical mechanics as an approximate model of an underlying quantum reality, then when considering initial distributions *ρi*(Γ) we should allow only such distributions as are consistent with the

uncertainty principle. To impose this constraint, let us imagine dividing phase space into cells of volume *hN*. A distribution *ρ*(Γ) that is consistent with the uncertainty principle is one that is uniform within any such cell, but whose value differs from cell to cell: any finer-grained structure is offensive to the uncertainty principle. For such a distribution, we have *pk* = *<sup>h</sup>Nρ*(<sup>Γ</sup>*k*), where Γ*k* is a representative point in cell *k* and *pk* = <sup>Γ</sup>∈cell *k d*Γ *ρ*(Γ) is the probability to find the system in that cell. The Shannon entropy of this distribution is given by

$$\mathcal{S}^{\mathbb{C}}[\rho] = -\int \rho(\Gamma) \ln[h^N \rho(\Gamma)] d\Gamma = -\sum\_{k} p\_k \ln p\_k \ge 0 \tag{60}$$

where S*c* = 0 if and only if *pk* = *δkl* for some cell *l*.

If we thus reject distributions with negative entropy as being incompatible with the uncertainty principle, then Equation (55) is replaced by min*ρ*∈Σ*<sup>c</sup>* S*c*[*ρ*(Γ)] = 0, and Equation (58) becomes

$$\max\_{\rho\_i \in \Sigma^{\varepsilon}} \mathcal{W}^{c \star} [\rho\_i(\Gamma)] = \beta^{-1} \mathcal{S}^{\varepsilon} [\sigma(\Gamma)] \,. \tag{61}$$

Thus, after imposing consistency with the uncertainty principle (in an admittedly heuristic fashion), we conclude that for both the quantum equivalence class Σ*q* and the classical equivalence class Σ*<sup>c</sup>*, the maximum extractable work is given by the entropy of the diagonal or homogeneous state, multiplied by *β*−<sup>1</sup> (Equations (57) and (61)).

Throughout the following section, and in Section 5, we impose the constraint S*c*[*ρ*] ≥ 0 on the initial classical phase space distribution, to exclude states that are incompatible with the uncertainty principle.

#### *4.3. A Fair Comparison*

The final step in making a fair comparison between quantum and classical work extraction is to establish a correspondence between equivalence classes Σ*q* and Σ*<sup>c</sup>*. That is, we want to establish a correspondence between quantum and classical energy distributions. There is no unique way to do this, as energy takes on discrete values in one case and continuous values in the other. As a reasonable way to proceed, let us choose a real function *<sup>κ</sup>*(·) ≥ 0 with the property that both K*q* = Tr *κ*(*H*<sup>ˆ</sup> 0) and K*c* = *d*Γ *<sup>κ</sup>*(*<sup>H</sup>*0(Γ)) are finite. We then define the diagonal quantum and homogeneous classical states

$$\mathcal{O}\_{\mathbb{K}} = \frac{\kappa(\varprojlim \varprojlim \varprojlim)}{\mathbb{K}^{\mathfrak{q}}} = \sum\_{\mathfrak{n}} \frac{\kappa(\mathfrak{e}\_{\mathfrak{n}})}{\mathbb{K}^{\mathfrak{q}}} |\mathfrak{n}\rangle\langle\mathfrak{n}| \quad , \quad \sigma\_{\mathfrak{k}}(\varGamma) = \frac{\kappa(H\_{0}(\varGamma))}{\mathbb{K}^{\mathfrak{c}}} \tag{62}$$

along with the associated energy equivalence classes

$$
\Sigma^{\eta}[\mathbf{x}] \quad = \quad \{\not\rho \,|\,\text{diag}\,\rho = \vartheta\_{\mathbf{x}}\} \tag{63a}
$$

$$\Sigma^{\mathfrak{c}}[\mathfrak{x}] \; := \; \{ \rho(\Gamma) \, | \, \text{diag} \, \rho = \sigma\_{\mathfrak{x}} \; ; \; \mathcal{S}^{\mathfrak{c}}[\rho] \ge 0 \}. \tag{63b}$$

The equivalence class <sup>Σ</sup>*<sup>q</sup>*[*κ*] contains all quantum states with diagonal density matrix elements *ρnn* = *<sup>κ</sup>*(*n*)/K*q*, whereas <sup>Σ</sup>*<sup>c</sup>*[*κ*] contains every classical state with energy distribution

$$
\eta(E) = \frac{\kappa(E)\Omega(E)}{\mathcal{K}^c}.\tag{64}
$$

Thus, a given choice of *<sup>κ</sup>*(·) specifies both a quantum and a classical energy distribution. As an example, for the choice *κ*(*x*) = *<sup>e</sup>*<sup>−</sup>*β<sup>x</sup>*, the reference states are *σ*ˆ*κ* = *π*ˆ 0 and *σκ*(Γ) = *<sup>π</sup>*0(Γ), and the energy equivalence classes are the thermal sets defined earlier: <sup>Σ</sup>*<sup>q</sup>*[*κ*] = Π*q* and <sup>Σ</sup>*<sup>c</sup>*[*κ*] = Π*<sup>c</sup>*.

In the semiclassical limit *h* → 0, as the level spacing between adjacent energy eigenvalues approaches zero, the normalized energy distribution associated with <sup>Σ</sup>*<sup>q</sup>*[*κ*] is conveniently written as *ξ*(*E*) = *<sup>κ</sup>*(*E*)*g*(*E*)/K*q*, where *g*(*E*) = ∑*n δ*(*E* − *n*) is the quantum

density of states. In turn, *g*(*E*) *dE* is approximated by the number of cells of volume *h<sup>N</sup>* that fit into the classical phase space volume between *E* and *E* + *dE*, for small *dE*. Equivalently,

$$\lim\_{h \to 0} h^N g(E) = \Omega(E) \tag{65}$$

where Ω(*E*) is the classical density of states, Equation (39). Hence, the quantum energy distribution is, semiclassically,

$$\xi(E) = \frac{\kappa(E)\Omega(E)}{h^N \mathbb{K}^q}.\tag{66}$$

Since both the classical and quantum energy distributions *η*(*E*) and *ξ*(*E*) (Equations (64) and (66)) are normalized to unity, we have

$$\lim\_{h \to 0} h^N \mathcal{K}^q = \mathcal{K}^c. \tag{67}$$

ˆ

From Equations (64), (66) and (67), we conclude that in the semiclassical limit *h* → 0, the discrete energy distribution associated with the equivalence class <sup>Σ</sup>*<sup>q</sup>*[*κ*] approaches the continuous distribution associated with <sup>Σ</sup>*<sup>c</sup>*[*κ*]. In this sense, we view <sup>Σ</sup>*<sup>q</sup>*[*κ*] and <sup>Σ</sup>*<sup>c</sup>*[*κ*] as having equivalent energy distributions.

Now, finally, for a given quantum reference Hamiltonian *H* 0 and its classical counterpart *<sup>H</sup>*0(Γ), and for a given choice of the function *<sup>κ</sup>*(·), let

$$\mathcal{W}^{\mathfrak{q}\star}\_{\max}[\mathbf{x}] = \max\_{\rho\_i \in \Sigma^{\mathfrak{q}}[\mathbf{x}]} \mathcal{W}^{\mathfrak{q}\star}(\rho\_i) \quad \text{and} \quad \mathcal{W}^{\mathfrak{c}\star}\_{\max}[\mathbf{x}] = \max\_{\rho\_i \in \Sigma^{\mathfrak{c}}[\mathbf{x}]} \mathcal{W}^{\mathfrak{c}\star}[\rho\_i(\Gamma)] \tag{68}$$

denote the maximum quantum and classical work that can be extracted during a cyclic, isoenergetic (in the sense of Equations (13) and (42)) process, for initial energy distributions determined by *<sup>κ</sup>*(·). We assert that by comparing the values of <sup>W</sup>*q*max[*κ*] and W*<sup>c</sup>*max[*κ*], in the semiclassical limit *h* → 0, we make a fair comparison between quantum work that can be extracted from coherences, and classical work that can be extracted from inhomogeneities.

From Equations (57), (61) and (63), we have

$$\mathcal{W}^{\mathfrak{f}^{\bullet}}\_{\max}[\kappa] = \mathcal{J}^{-1}\mathcal{S}^{\mathfrak{q}}(\hat{\sigma}\_{\mathfrak{k}}) \quad , \quad \mathcal{W}^{\mathfrak{c}\star}\_{\max}[\kappa] = \mathcal{J}^{-1}\mathcal{S}^{\mathfrak{c}}[\sigma\_{\mathfrak{k}}(\Gamma)] \, . \tag{69}$$

therefore, let us inspect the difference between these two values,

$$
\Delta \mathcal{W}^\star = \mathcal{W}^{q\star}\_{\text{max}}[\mathbf{x}] - \mathcal{W}^{c\star}\_{\text{max}}[\mathbf{x}]\_\prime \tag{70}
$$

in the limit *h* → 0. Following the semiclassical approach used above, we obtain

$$\lim\_{h\to 0} \Delta \mathcal{W}^\* \quad = \beta^{-1} \lim\_{h\to 0} \left\{ -\sum\_{\overline{\lambda}\in \mathcal{I}} \frac{\kappa(\epsilon\_0)}{\overline{\lambda}!} \ln \left[ \frac{\kappa(\epsilon\_0)}{\overline{\lambda}!} \right] + \int \frac{\kappa(\mathcal{U}\_0)}{\overline{\lambda}!} \ln \left[ h^N \frac{\kappa(\mathcal{U}\_0)}{\overline{\lambda}!} \right] d\Gamma \right\} $$

$$= \beta^{-1} \lim\_{h\to 0} \left\{ -\int g(E) \frac{\kappa(E)}{\overline{\lambda}!} \ln \left[ \frac{\kappa(E)}{\overline{\lambda}!} \right] dE + \int \Omega(E) \frac{\kappa(E)}{\overline{\lambda}!} \ln \left[ h^N \frac{\kappa(E)}{\overline{\lambda}!} \right] dE \right\} \tag{71} $$

$$= \beta^{-1} \lim\_{h\to 0} \left\{ -\int \Omega(E) \frac{\kappa(E)}{\overline{\lambda}!} \ln \left[ h^N \frac{\kappa(E)}{\overline{\lambda}!} \right] dE + \int \Omega(E) \frac{\kappa(E)}{\overline{\lambda}!} \ln \left[ h^N \frac{\kappa(E)}{\overline{\lambda}!} \right] dE \right\} $$

$$= 0.$$

Here, Equation (62) has been combined with the expressions for von Neumann and Shannon entropy (Equations (9) and (31)) on the first line; the sum over energy eigenstates and the integral over phase space have been replaced by energy integrals on the second line; and Equations (65) and (67) have been used to ge<sup>t</sup> to the third line.

For *κ*(*x*) = *<sup>e</sup>*<sup>−</sup>*β<sup>x</sup>*, Equation (71) can alternatively be established from the result (see Equations (4) and (26))

$$\Delta \mathcal{W}^\* = \beta^{-1} \left\{ \mathcal{S}^q(\text{\textit{\textit{\textit{\textit{\textit{X}}}}}) - \mathcal{S}^c[\text{\textit{\textit{\textit{\textit{\textit{X}}}}}] \text{\textit{\textit{\textit{\textit{\textit{\textit{X}}}}}}}] \right\} = \beta \frac{\partial}{\partial \beta} (\mathcal{F}^{\text{q,c}\text{q}} - \mathcal{F}^{c,\text{eq}}) = \left( -\frac{1}{\beta} + \frac{\partial}{\partial \beta} \right) \log \frac{Z^c}{h^N Z^{\text{\textit{\textit{\textit{\textit{\textit{\textit{\textit{\bullet}}}}}}}} \tag{72}$$

where *Zq* and *Zc* are equilibrium partition functions. Taking the limit *h* → 0 and using the known result [47,53–55] that (for kinetic-plus-potential Hamiltonians) *h<sup>N</sup> Zq* can be expanded in a power series of *h* whose first term is exactly the classical partition function *Zc*, the right side of Equation (72) vanishes.

From Equation (71), we conclude that in the semiclassical limit, the maximal work that can be extracted from the energy coherences of a quantum state *ρ*ˆ*i* ∈ <sup>Σ</sup>*<sup>q</sup>*[*κ*] is the same as the maximal work that can be extracted from the energy-shell inhomogeneities of a classical state *ρi*(Γ) ∈ <sup>Σ</sup>*<sup>c</sup>*[*κ*]. In both situations, the work is maximized by starting in the state of least entropy within Σ*q* or Σ*<sup>c</sup>*, then quasistatically removing the coherences or inhomogeneities. This result leads us to conclude that, within our framework for comparing quantum and classical systems, quantum coherences offer no particular thermodynamic advantage over classical inhomogeneities.

#### **5. Dropping the Isoenergetic Constraint**

In the previous sections, we have imposed the isoenergetic constraint, namely that the initial and final energy distributions are identical (Equations (13) and (42)). Let us now drop this constraint and pose the following question. For a quantum or classical system described by an initial Hamiltonian *H*ˆ 0 or *<sup>H</sup>*0(Γ), in the presence of a thermal bath at temperature *β*−1, what is the maximum work that can be extracted during a cyclic process if the energy distribution of the initial state is determined by a given function *<sup>κ</sup>*(·)?

In the quantum case, we first let <sup>W</sup>*q*†(*ρ*<sup>ˆ</sup>*i*) denote the maximum work extracted for a given initial state *ρ*<sup>ˆ</sup>*i*—this quantity is analogous to W(*ρ*<sup>ˆ</sup>*i*) (Section 2) but without the constraint diag *ρ*ˆ*f* = diag *ρ*<sup>ˆ</sup>*i*. From Equations (5) and (10) and the non-negativity of the Kullback–Leibler divergence, we have

$$\begin{split} \mathcal{W}^{q\dagger}(\boldsymbol{\hat{\rho}\_{i}}) &\leq \mathcal{F}^{q}(\boldsymbol{\hat{\rho}\_{i}}, \boldsymbol{\hat{H}\_{0}}) - \mathcal{F}^{q}(\boldsymbol{\hat{\rho}\_{f}}, \boldsymbol{\hat{H}\_{0}}) \\ &\leq \mathcal{F}^{q}(\boldsymbol{\hat{\rho}\_{i}}, \boldsymbol{\hat{H}\_{0}}) - \mathcal{F}^{q}(\boldsymbol{\hat{\tau}\_{0}}, \boldsymbol{\hat{H}\_{0}}) \\ &= \mathcal{U}\_{i}^{q} - \boldsymbol{\beta}^{-1} \mathcal{S}\_{i}^{q} - \mathcal{F}\_{0}^{q, \text{eq}} \end{split} \tag{73}$$

where the inequality on the first line is valid for any final state *ρ*ˆ*f* , and F*q*,eq 0 ≡ <sup>F</sup>*q*,eq(*H*<sup>ˆ</sup> 0). As shown in Appendix A, the bound obtained in Equation (73) is saturated by the protocol

$$\hat{H}(t) = \begin{cases} \hat{H}\_0 & t \le 0 \\ -\beta^{-1} \ln \left[ (1 - \lambda) \hat{\rho}\_i + \lambda \, e^{-\beta \hat{H}\_0} \right] & 0 < t \le \tau \\ \hat{H}\_0 & \tau \le t \end{cases} \tag{74}$$

where *λ* ≡ *t*/*τ* and the process is quasistatic: *τ* → ∞. (Note that there is no quench at *t* = *τ*.) Since the bound can be saturated, and <sup>W</sup>*q*†(*ρ*<sup>ˆ</sup>*i*) was defined as the maximum work that can be extracted, we simply write

$$\mathcal{W}^{q\dagger}(\not p\_i) = \mathcal{U}\_i^q - \beta^{-1} \mathcal{S}\_i^q - \mathcal{F}\_0^{q, \text{eq}}.\tag{75}$$

Now, maximizing this quantity over all *ρ*ˆ*i* ∈ <sup>Σ</sup>*<sup>q</sup>*[*κ*], we have

$$\begin{split} \mathcal{W}\_{\text{max}}^{q\dagger}[\kappa] &= \max\_{\hat{\rho}\_{\hat{i}} \in \Sigma^{q}[\kappa]} \left( \mathcal{U}\_{i}^{q} - \beta^{-1} \mathcal{S}\_{i}^{q} - \mathcal{F}\_{0}^{q, \text{eq}} \right) \\ &= \max\_{\hat{\rho}\_{\hat{i}} \in \Sigma^{q}[\kappa]} \left( \mathcal{U}^{q}[\kappa] - \beta^{-1} \mathcal{S}\_{i}^{q} - \mathcal{F}\_{0}^{q, \text{eq}} \right) \\ &= \mathcal{U}^{q}[\kappa] - \mathcal{F}\_{0}^{q, \text{eq}} \end{split} \tag{76}$$

 where U*<sup>q</sup>*[*κ*] ≡ (1/K*q*) ∑*n <sup>κ</sup>*(*n*)*n* is the average energy for every state *ρ*ˆ*i* ∈ <sup>Σ</sup>*<sup>q</sup>*[*κ*], and we have used Equation (50b) to arrive at the third line.

As a consistency check, we combine Equations (69) and (76) with Equations (7) and (10) to obtain

$$\begin{split} \mathcal{W}\_{\text{max}}^{q\dagger}[\kappa] - \mathcal{W}\_{\text{max}}^{q\star}[\kappa] &= \mathcal{U}^{q}[\kappa] - \mathcal{F}\_{0}^{q\text{eq}} - \mathcal{\beta}^{-1} \mathcal{S}^{q}(\hat{\sigma}\_{\text{k}}) \\ &= \mathcal{F}^{q}(\hat{\sigma}\_{\text{k}}, \hat{H}\_{0}) - \mathcal{F}\_{0}^{q\text{eq}} = D(\hat{\sigma}\_{\text{k}} | \hat{\pi}\_{0}) \ge 0 \end{split} \tag{77}$$

where *σ*ˆ*κ* is the unique diagonal state belonging to <sup>Σ</sup>*<sup>q</sup>*[*κ*]. Thus, <sup>W</sup>*<sup>q</sup>*†max[*κ*] ≥ <sup>W</sup>*q*max[*κ*], which makes sense: the maximum work that we can extract without imposing the constraint diag *ρ*ˆ*f* = diag *ρ*ˆ*i* must be no less than the maximum work we can extract with the constraint.

In the classical case, essentially identical calculations—which we do not reproduce here—lead to the result

$$\mathcal{W}^{c\dagger}\_{\text{max}}[\kappa] = \mathcal{U}^c[\kappa] - \mathcal{F}^{c,\text{eq}}\_0 \tag{78}$$

where <sup>W</sup>*<sup>c</sup>*†max[*κ*] is the maximum work that can be extracted over all initial states *ρi*(Γ) ∈ <sup>Σ</sup>*<sup>c</sup>*[*κ*], without imposing Equation (42), and U*<sup>c</sup>*[*κ*]=(1/K*c*) *d*Γ *<sup>κ</sup>*(*<sup>H</sup>*0)*<sup>H</sup>*0 is the average energy for every state in <sup>Σ</sup>*<sup>c</sup>*[*κ*]. Following steps similar to those of Section 4.3, we obtain

$$\begin{aligned} \lim\_{h \to 0} \mathcal{U}^{q}[\kappa] &= \lim\_{h \to 0} \frac{1}{\mathcal{K}^{q}} \sum\_{n} \kappa(\epsilon\_{n}) \epsilon\_{n} \\\\ &= \lim\_{h \to 0} \int dE \,\mathcal{g}(E) \frac{\kappa(E)}{\mathcal{K}^{q}} E \\\\ &= \int dE \,\Omega(E) \frac{\kappa(E)}{\mathcal{K}^{c}} E = \frac{1}{\mathcal{K}^{c}} \int d\Gamma \,\kappa(H\_{0}) \,H\_{0} = \mathcal{U}^{c}[\kappa] \end{aligned} \tag{79}$$

and

$$\lim\_{h \to 0} \left( \mathcal{F}\_0^{\varepsilon, \text{eq}} - \mathcal{F}\_0^{q, \text{eq}} \right) = -\beta^{-1} \lim\_{h \to 0} \ln \frac{Z^{\varepsilon}[H\_0]}{h^N Z^q(\hat{H}\_0)} = 0 \tag{80}$$

using Equation (67), with *κ*(*x*) = *<sup>e</sup>*<sup>−</sup>*β<sup>x</sup>*.

> Defining ΔW† ≡ <sup>W</sup>*<sup>q</sup>*†max[*κ*] − <sup>W</sup>*<sup>c</sup>*†max[*κ*], Equations (76) and (78)–(80) give us

$$\lim\_{\text{hr}\to 0} \Delta \mathcal{W}^{\dagger} = 0 \tag{81}$$

which is the counterpart of Equation (71), after abandoning the constraint of equal initial and final energy distributions. We again conclude that quantum coherences provide no inherent thermodynamic advantage over classical inhomogeneities, in the semiclassical limit.

## **6. Conclusions**

In Sections 2 and 3 of this paper, we argued that quantum energy coherences (as shown earlier [9]) and classical energy shell inhomogeneities represent thermodynamic resources, which can be leveraged to deliver work. In Sections 4 and 5, we argued that a fair comparison shows these resources to be equivalent: in the semiclassical limit, and for a

given initial energy distribution, the amount of work that can be extracted from quantum coherences is the same as the amount that can be extracted from classical inhomogeneities.

Our study has focused on processes during which the system of interest is in contact with a thermal reservoir, and here (as we have seen) the free energy F plays an important role. Sone and Deffner [18] have recently carried out a similar investigation for isolated quantum and classical systems, in which case *ergotropy* (defined in Ref. [1] for quantum systems and in Ref. [18] for classical systems) plays a role analogous to free energy in our paper. In Ref. [18], as in our paper, energy-shell inhomogeneities are classical counterparts of quantum energy coherences.

In making our comparison in Sections 4 and 5, we invoked a quantum–classical correspondence based on canonical quantization, in which the system of interest is described by coordinates *x*1, *x*2, ··· and conjugate momenta *p*1, *p*2, ··· , which are either quantum operators or classical observables. For such systems, the classical phase space is unbounded and the quantum Hilbert space is infinite-dimensional.

However, in the quantum thermodynamics literature one often encounters systems with finite-dimensional Hilbert spaces, such as the illustrative qubit example analyzed in Ref. [9]. It then seems natural to take, as the quantum system's counterpart, a discrete-state classical system of equal dimensionality. Thus, a qubit's counterpart may be taken to be a classical bit. For such discrete-state systems there is no opportunity to introduce a classical analogue of quantum coherences, as the statistical state of a classical *D*-state system is specified *entirely* by the probabilities *P*1, ··· *PD*, and these are in one-to-one correspondence with the diagonal elements of the corresponding quantum system's density matrix *ρ*ˆ. In this situation, it seems that quantum coherences really do provide a unique thermodynamic resource that is unavailable to classical counterparts.

This conclusion, however, is misleading, as an apparently discrete-state classical system is in reality a coarse-grained version of a more microscopically detailed system. For example, an effective classical bit can be obtained by coarse-graining a classical particle in a double-well potential, such that the location *x* of the particle in the left (right) well indicates a bit value of 0 (1). The apparent quantum thermodynamic advantage—due to coherences—arises in this case because potentially useful classical information (e.g., how the particle's potential energy depends on its location *x*) has been thrown out in the process of coarse-graining from the double well to the bit. Comparing a qubit—an intrinsically twostate quantum system—with an effective classical two-state system obtained by discarding microscopic information, is an apples-to-oranges comparison.

There is no generally applicable procedure for identifying a proper classical counterpart of a quantum system with a finite-dimensional Hilbert space. It is instructive, however, to consider the simplest case of a spin-1/2 particle (qubit) in a magnetic field, governed by a Hamiltonian *H*ˆ = *g* **B** · **s**ˆ, where **s**ˆ = (*h*¯ /2)( *<sup>σ</sup>*<sup>ˆ</sup>*x*, *<sup>σ</sup>*<sup>ˆ</sup>*y*, *<sup>σ</sup>*<sup>ˆ</sup>*z*). In the absence of a thermal bath, the unitary dynamics in the Heisenberg representation are given by the equations of motion

$$i\hbar \frac{d}{dt} \mathfrak{s}\_{\mathbb{H}} = [\mathfrak{s}\_{\mathbb{H}}, \hat{H}] \tag{82}$$

where the right side is evaluated using the commutation relations

$$\left[\mathcal{S}\_{\rangle}, \mathcal{S}\_{k}\right] = i\hbar \,\varepsilon\_{jkl}\,\mathcal{S}\_{l} \tag{83}$$

and *<sup>ε</sup>jkl* is the Levi-Civita symbol. Kammerlander and Anders [9] showed how work can be extracted from energy coherences in such a system, using a protocol involving quenches and the quasistatic variation of **B**, along with coupling to a thermal bath.

As a possible classical counterpart, instead of a two-state bit let us consider a system whose microscopic state is described by a vector **S** = (*Sx*, *Sy*, *Sz*) of fixed magnitude, governed by a Hamiltonian *H* = *g* **B** · **S**, evolving under the Poisson bracket formulation of Hamiltonian dynamics,

$$\frac{d}{dt}\mathbf{S} = \{\mathbf{S}, H\} \tag{84}$$

with

$$\{\mathcal{S}\_{\dot{l}\prime}\mathcal{S}\_{k}\} = \varepsilon\_{jkl}\mathcal{S}\_{l}.\tag{85}$$

The phase space for this classical system is bounded: it is the two-dimensional surface of a sphere of radius |**S**|. An energy shell is represented by a circle on that sphere, oriented along the **B**-direction. The dynamics given by Equation (84) describe an isolated system, and would have to be supplemented by appropriate terms in order to include the effects of contact with a thermal bath. It would then be interesting to investigate classical protocols designed to extract work from an initial distribution that is inhomogeneous on the energy shells, and to compare this classical situation with the quantum case of Ref. [9].

We note that the approach described in the previous paragraphs is readily extended to a system composed of *N* > 1 spins, interacting both with external fields and among themselves, e.g., through Hamiltonian terms of the form *cmn***s**<sup>ˆ</sup>*m* · **s**ˆ*n* or *cmn***S***m* · **S***<sup>n</sup>*. Thus, comparisons between quantum and classical work extraction can be extended to multi-spin systems, within this framework. For example, it has been demonstrated that quantum correlations within a many-body system can be utilized for extracting work [56–59], and it would be pertinent to study whether one can leverage classical correlations and inhomogeneities in a similar way. Such comparisons may further elucidate whether thermodynamic advantages can be identified that are unique to quantum systems.

**Author Contributions:** Conceptualization, A.S., K.S. and C.J.; Investigation, A.S., K.S. and C.J.; Writing—original draft, A.S. and C.J.; Writing—review & editing, K.S. All authors have read and agreed to the published version of the manuscript

**Funding:** This research was funded by the United States National Science Foundation under gran<sup>t</sup> DMR-1506969 (AS, CJ).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We gratefully acknowledge stimulating discussions with Janet Anders, Sebastian Deffner, David Goldwasser, Wade Hodson, Paul Riechers, and Nicole Yunger Halpern.

**Conflicts of Interest:** The authors declare no conflict of interest.
