Charging assisted by thermalization

A system in thermal equilibrium with a bath will generally be in an athermal state, if the system-bath coupling is strong. In some cases, it will be possible to extract work from that athermal state, after disconnecting the system from the bath. We use this observation to devise a battery charging and storing unit, simply consisting of a system, acting as the battery, and a bath. The charging cycle—connect, let thermalize, disconnect, extract work—requires very little external control and the charged state of the battery, being a part of global thermal equilibrium, can be maintained indeﬁnitely and for free. The efﬁciency, deﬁned as the ratio of the extractable work stored in the battery and the total work spent on connecting and disconnecting, is always (cid:2) 1, which is a manifestation of the second law of thermodynamics. Moreover, coupling, being a resource for the device, is also a source of dissipation: the entropy production per charging cycle is always signiﬁcant, strongly limiting the efﬁciency in all coupling strength regimes. We show that our general results also hold for generic microcanonical baths. We illustrate our theory on the Caldeira-Leggett model with a harmonic oscillator (the battery) coupled to a harmonic bath, for which we derive general asymptotic formulas in both weak and ultrastrong coupling regimes, for arbitrary Ohmic spectral densities. We show that the efﬁciency can be increased by connecting several copies of the battery to the bath. Finally, as a side result, we derive a general formula for Gaussian ergotropy, that is, the maximal work extractable by Gaussian unitary operations from Gaussian states of multipartite continuous-variable systems.


I. INTRODUCTION
The second law of thermodynamics, as per the Kelvin-Planck formulation [1], states that no work can be extracted in a cyclic manner from a system in thermal equilibrium.On the other hand, in the presence of strong interactions, the reduced state of a subsystem of a thermal system will not typically be thermal [2,3].Thereby, not limited by the second law anymore, it may be possible to cyclically extract nonzero work from such a subsystem, when manipulated in separation from the rest of the system [4,5].
In precise terms, the cyclic processes referred to above are processes where the system's state evolves according to an externally driven time-dependent Hamiltonian, the value of which at the end is equal to that at the beginning.The maximal amount of work extractable from a system by such processes is called ergotropy [6] (see Appendix A for a detailed definition).A system with zero ergotropy is called passive [7,8], and that with nonzero ergotropy is called active.In these terms, the Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
previous paragraph reads: Thermal states are passive; however, the reduced state of a subsystem of a thermal system can be active with respect to the local Hamiltonian.
Inspired by these basic observations and the setup of Ref. [9], we introduce a battery charging cycle, the central idea of which is to connect a system (the "battery") in a passive ("depleted") state to a thermal bath, and wait until they thermalize.This will prepare the system in an active ("charged") state, from which we will be able to extract work after the system is disconnected from the bath.Of course, this cycle does not violate the Kelvin-Planck formulation of the second law since connecting and disconnecting the system will cost work, which will have to be provided by external agents [4,5,10].
Since large thermal baths generically thermalize finite-size systems that come in contact with them (see the discussion in Sec.II and Appendix B), our device offers two main advantages: (i) the creation of the battery's charged state requires no fine external control and, since thermalization is the process preparing that state, is robust against minor variations of the system-bath interaction; (ii) it costs nothing to maintain the charged state for as long as might be needed-it is the stationary state of system-bath interaction.These advantages are not simultaneously met in other battery designs which either require finely tuned external fields to perform unitary charging operations on the depleted state of the battery and assume that the battery is isolated after it is charged [11,12]; or require system-bath interaction engineering [9]; or, to prevent the battery from leaking the charge, either rely on fragile symmetries of the system-bath interaction [13,14] or actively manipulate the battery [15][16][17].
Although, at the beginning of the cycle, the state of the total system is not thermal-especially when the bath is not in a canonical (a.k.a., Gibbs) state, e.g., when it is in a microcanonical state-and therefore may be active, we show for a generic thermalizing bath that, due to cyclicity, the total work one has to spend on connecting and disconnecting the system from the bath, W c:d , is always larger than the maximal work one can extract from the system after it is detached from the bath, i.e., its ergotropy E. This means that the efficiency, defined as the ratio of the energy one is able to extract and the energy one has to invest for that: η := E/W c:d , is 1.We show that this is nevertheless a consequence of the passivity of Gibbs states, even in those cases when the bath is microcanonical, and hence active [18], due to the so-called equivalence of canonical and microcanonical states [19][20][21].
We support our general findings with a detailed calculation of all relevant quantities for the Caldeira-Leggett model [22], where the system is a harmonic oscillator and couples to a bath made of harmonic oscillators.We study the operation of the device in all relevant parameter regimes.In particular, we show that, other parameters fixed, the device is most efficient in the intermediate range of couplings.On the route of exploring different parameter regimes, we derive exact asymptotic expansions for the covariance matrix of the oscillator-for general Ohmic spectral densities-in the weak-coupling and ultrastrong-coupling.In the high-temperature limit, we prove an equipartition result for both the kinetic and potential energies of the oscillator, for arbitrary spectral densities and strengths of coupling.Found general expansions can be useful beyond the type of problems studied in this work.
Let us note in passing that, taking a different standpoint, our device can be viewed as a single-bath machine, and our above definition of its efficiency is standardly used in other types of single-bath machines [9,[23][24][25][26][27][28].Importantly, however, the working principle of our device is fundamentally different from that of other single-bath machines such as molecular motors or force-to-force converters [23][24][25][26][27][28] in that, as opposed to those machines, in our case, the thermalizing effect of the environment has a constructive role in generating the output work.More generally, our device utilizes an energy storage mechanism which is not characteristic of the conventional molecular motors and force-to-force converters.

II. THE CYCLE
Let us formalize the description of the device and its charging cycle in the introduction.First, we introduce the total system-bath Hamiltonian, where H s is the system Hamiltonian, H B is the bath Hamiltonian, H I is the interaction Hamiltonian (possibly containing a system renormalization term), and H 0 = H s + H B is the bare, noninteracting total Hamiltonian of the system and the bath.We assume the bath to be large, i.e., consisting of a large N 1 number of, generically, interacting constituents, and "complex" enough to thermalize nonmacroscopic systems that come in contact with it (see below for precise definitions).We denote the state of the system before interacting with the bath by ρ 0 and the initial state of the bath by R B , so that the initial state of the total system is (2) Also, we introduce ∞ = e −iH T t ∞ e iH T t ∞ , where t ∞ is a duration long enough to equilibrate the overall system.With these, the (Hamiltonian) cycle we consider consists of the following strokes: Here, the unitary operator, acts solely on the system, and is such that U E extracts the full ergotropy from ρ ∞ s = Tr B ∞ , with respect to H s , and leaves the system in the passive state, which is the state in which the system enters the next cycle.(See Appendix A for precise definitions of ergotropy and passivity.)As we mentioned, connecting and disconnecting the system to the bath has a work cost W c:d := W c + W d , where W c is the cost of stroke 1 and W d of stroke 3 .
Remarkably, the stroke 2 can last as long as one wishes.Being the stage at which the active state of the system is prepared, it gives one considerable flexibility in practical situations, as there is no need of fine control-one just leaves the system in contact with a bath, and, whenever work is needed, one abruptly detaches the system (stroke 3 ) and performs the extraction process (stroke 4 ), which, when the bath is thermalizing (see below), is conveniently independent of the initial state of the system.Moreover, the time interval between 3 and 4 is also arbitrary, since the ergotropy of a system evolving under the influence of its own Hamiltonian is time-independent (see Appendix A).However, the unitary operation extracting the ergotropy does depend on time, so, although one does not need to perform 4 immediately after 3 , fine tuning is necessary in any case.Importantly, in contrast to the lack of upper bound on the duration of stroke 2 , there is a lower bound: the equilibration time.Although it can sometimes be quite short [29][30][31], relaxation generically takes nonnegligible amount of time which depends on the details of the system-reservoir interaction and the initial state.This time-the charging time in our case-is generally expected to be a decreasing function of the coupling strength [22,31,32], especially in those situations when the initial state commutes with the bare Hamiltonian [31], which is what we will typically have in our protocol.Indeed, starting from the second cycle, the initial state is ρ p ⊗ R B , and passive states always commute with the Hamiltonians with respect to which they are passive ([ρ p , H s ] = 0; see Refs.[6][7][8] or Appendix A).Moreover, R B will generically-although not necessarily-be either a canonical or a microcanonical state, and both commute with H B [cf. Eqs. ( 9) and (10) below].This is a fortunate situation for our device: the need for strong (but not too strong-see Sec.III A) coupling in order to operate comes with the benefit that it also ensures fast charging.
In our further analysis, we will assume-and this is a crucial assumption-that the joint evolution of the system and bath is thermalizing in the following sense: Say, S is a finite subsystem of the system-plus-bath composite belonging to s ∪ supp(H I ), then the long-time limit of the reduced state of S should be where More precisely, Eq. ( 6) should be understood as a statement about the long-time average of the state: lim where t = e −iH T t e iH T t and • is the trace norm [33], should be small and tend to zero as the size of the bath (N) increases (see Refs. [29,31,34,35] for examples of such bounds).By the Markov inequality, this means that the system's state will be close to Tr T\S τ T for most of the time.This is a rather weak assumption as thermalization is ubiquitous in macroscopic systems [19][20][21]29,30,[34][35][36][37][38][39].When the bath starts in a canonical (Gibbs) state, thermalization in the above sense is sometimes referred to as "return to equilibrium."It has been rigorously proven in several generic scenarios.First is when, to an infinitely large, strictly continuous bosonic bath, is linearly coupled a continuous-variable system (e.g., the Caldeira-Leggett model [22]) [40,41] or a system with a finite-dimensional Hilbert space (e.g., the spin-boson model [42]) [37,43].The other scenario is when the bath is a many-body system with shortrange interactions, away from criticality (i.e., with a finite correlation length), satisfying some transport conditions, e.g., nonzero Lieb-Robinson velocity ("speed of sound") [44] or absence of many-body localization [39].Note that exponential decay of correlations is very common in short-range interacting systems and is guaranteed for arbitrary Fermi systems at nonzero temperature [45], vacuum states of gapped lattice Hamiltonians [46], general lattice systems above a critical temperature [47], and is often related to the presence of a finite Lieb-Robinson velocity in the system [48][49][50].Early results about thermalization in systems with short-range interactions deal with the thermalization of locally perturbed translationally invariant infinite chains satisfying certain ergodicity-like properties (see, e.g., Ref. [36] and references therein).More recent results for lattice systems are based on the ideas of equilibration [29,31,34] and the equivalence of canonical and microcanonical ensembles and Berry-Esseen-type concentration bounds [19][20][21]35,51], and therefore require only few generic assumptions in order to hold.The most general rigorous proof of thermalization in the sense of Eq. ( 6) holds for any finite-range H T (i.e., containing at most k-body interaction terms, where k is some finite number) such that both τ T and τ B have exponentially decaying correlations [35], with the only additional (weak) requirement being that H T needs to have not too many repeating gaps in its spectrum (see Appendix B for more details).Using the fact that, with these conditions, microcanonical and canonical states are locally equivalent [21,51], in Appendix B, and combining results from Refs.[21,35], we show that subsystems thermalize (i.e., Eq. ( 6) holds) also when the bath starts in the microcanonical state [20,21,52]: where e n and |e n are, respectively, the eigenvalues and eigenvectors of H B , O ln 2d N O √ N is the microcanonical energy window, d (E , ) is the amount of energy levels in the interval [E − , E + ], and E (chosen to be macroscopic, i.e., ∝N) is the energy of the microcanonical state.Temperature is prescribed to the microcanonical state through E = Tr(H B τ B ), which, in view of Tr(H B τ B ) also being ∝N and monotonic with respect to β, defines a unique function β = β(E /N ).
In view of the thermalization assumption, let us point out that, whereas, similar to stroke 2 , stroke 1 can last as long as one wishes (W c will of course depend on the details of the switching protocol, but the end result of stroke 2 will always be the same), stroke 3 has to be a nonstationary process.Indeed, if performed in a slow, "quasi-equilibrium" manner, the system will at all times remain thermalized with the bath.Therefore, by the end of the disconnection process, when the interaction vanishes, it will simply be in a Gibbs state with respect to H s , hence, no work will be possible to extract during stroke 4 , rendering the cycle useless.
We also note that, since the Ohmic Caldeira-Leggett model with a Lorentz-Drude cut-off function can be mapped into a gapless harmonic lattice (where the system maps to a node) with polynomially decaying interactions [53], the thermalization result for the CL model [40,41] can be thought of as an extension of the above-described paradigm of [short-range]+[noncritical]→[thermalization] to critical systems with long-range interactions; in this context, the result about the microcanonical bath also being thermalizing is unlikely to hold.

A. The energetics of the cycles
Here we will study the energetics of the device and find its efficiency.Directly reading from Eqs. (3a)-(3d), where ρ ∞ s is given by Eq. ( 6).Mind the sign convention: W c and W d are works performed on the (composite) system, whereas E is the work extracted from the system.
Relying on the thermalization results discussed in the previous subsection, we will henceforth assume that the system-bath evolution is thermalizing for all finite subsystems of s ∪ supp(H I ), and therefore, in Eqs.(11a) and (11b), we will substitute R B by τ B [if R B is not τ B , e.g., when it is the microcanonical state (10)] and ∞ by τ T : Here we need to maintain caution, since, as is detailed in Appendix B, Eq. ( 6) generally holds only up to a correction O(N − ) (with some > 0) if supp(H I ) is finite.However, when supp(H I ) scales with N, the corrections may potentially accumulate into something nonnegligible.We discuss this further in Appendix B 1, where we show the conditions on H I that guarantee the correctness of Eq. ( 12); for the Caldeira-Leggett model, these reduce to an explicit condition on the spectral density.Now, defining "dissipated work" as which, according to our sign convention, corresponds to the total work performed during the cycle, and putting Eqs.(11a)-( 12) together, we find that where and the reason for this notation will become clear in what follows.
Starting from the second cycle, the information about the initial state of the system, ρ 0 , is lost: the system starts and ends in the state ρ p [see Eqs.(3e) and ( 5)], meaning that where i counts the cycles.In fact, in most physically relevant situations, w (1) will also be =0.Indeed, generically, , where S κ are some system operators and B κ are bath operators.The latter will typically be some sort of "field operators," and, above a critical temperature (and at any temperature in one-and two-dimensional systems with continuous symmetries), will generically have zero thermal averages: Tr[B κ τ B ] = 0 (see the discussion of the Mermin-Wagner theorem in, e.g., Ref. [54]), meaning that w (1) = 0 irrespective of ρ 0 .This is obviously the case for all quadratic bosonic [22,42] and fermionic baths [55].That said, note that we need not and will not make such an assumption in what follows.In fact, we could get rid of w (1) Note that this formula holds both when R B = τ B and R B = μ B (E , ) (and whenever the bath is thermalizing).Furthermore, introducing where β is determined from where S(ρ) = − Tr(ρ ln ρ) is the von Neumann entropy and, in the second equality, we noted that ρ p and ρ ∞ s have the same entropy because one is a unitary transformation of the other.Next, by a series of algebraic manipulations, presented in Appendix C, we prove the following identity: where S(ρ||τ ) = Tr[ρ(ln ρ − ln τ )] is the relative entropy and ) 0 is the mutual information [33] between the system and the bath in the state τ T .
All the terms in the right-hand side of Eq. ( 21) are nonnegative, therefore, 0, (22) which, given that T is the dissipated work, allows us to loosely interpret 0 as the entropy production of the cycle [56].In the following subsection, we will corroborate this interpretation by putting the energetics of the cycle in the context of the second law.

B. Discussion of the energetics
Ergotropy, E, is by definition nonnegative, and since T 0, the connection-disconnection work will also be nonnegative: In those cases when w (1) happens to be nonzero, it can be both positive and negative, therefore, (only) for the first cycle, W c:d 0 may not hold.Therefore, in the general spirit of thermodynamics, and particularly the example of other single-bath machines [9,23,24], we will define the efficiency of our device as an output/expenditure ratio: System-bath coupling is what powers our device, which is expressed in the fact that as follows from Eqs. (11a)-(11c).Suppose for a moment that there exists a single parameter, g, that controls the strength of the interaction (say, H I = gV ).Then, assume that is differentiable in g (this is not guaranteed as U E is a permutation matrix that depends on eigenvalue ordering).Now, since both the relative entropy and mutual information are nonnegative, their Taylor expansion for small g's cannot start with a ∝g term, because otherwise the g → −g transformation would change the sign of the quantity, for sufficiently small g's; therefore, T = O g 2k , where k 1 is a natural number.
For the same reason, since W c:d 0 and → 0 as g → 0, we also conclude that W c:d = O g 2m , where m 1 is another natural number (see Ref. [57] for an example of an explicit calculation of connection-disconnection work in a harmonic chain, with continuous switching, where m = 1).This means that, although both the nominator and denominator in Eq. ( 24) go to zero, their ratio, and therefore η, will either go to zero (if k > m), or remain finite but <1 (if k = m), or go to one (if k < m).The latter case is obviously the most interesting, however, we could not find such an example.Now we will link Eq. ( 22) to the second law.With this purpose, we note that, by virtue of Eqs. ( 12)- (14), the total work performed on the total system in n cycles is nT (when w (1) = 0, we need to add it as well, however, it will play no role for sufficiently large n's, therefore, we will omit it).By definition, −nT is the net extracted work in n cycles.On the other hand, when the bath is canonical, the bird-eye view on the cycle reveals a simple picture: the total system starts with (ρ 0 ⊗ τ , where (the subsequent argumentation does not rely on the baths being identical; we made that choice to merely simplify notation).Then, the overall Hamiltonian, B , undergoes a cyclic variation in time, thereby driving the overall system unitarily, and, at the end of the cycle, we are left with (U , where U is the unitary evolution operator generated by the cyclic variation of the Hamiltonian.Now, a well-known formulation of the second law (which was first noted in phenomenological thermodynamics in Refs.[1,52] and rigorously proven in the quantum regime in Ref. [58]), states that the maximal work that can be extracted by a cyclic variation of the Hamiltonian (that is, ergotropy) from (ρ 0 ⊗ τ, H s + H ) (where H is an arbitrary Hamiltonian and τ ∝ e −βH ) is upper-bounded by the difference of nonequilibrium free energies of the system: where τ s ∝ e −βH s and is the nonequilibrium free energy of the system with respect to a T -temperature bath.This fact is a consequence of the simple identity where U is the unitary evolution operator generated by the cyclic variation of the Hamiltonian, τ sB ∝ e −β(H s +H ) , and is the average extracted work.For sufficiently large baths, the upper bound in Eq. ( 26) is always reachable through a sequence of small quenches and thermalizations (namely, S(U ρ 0 ⊗ τU † τ sB ) can be made arbitrarily small; see, e.g., Ref. [59]).Note that, in view of Eq. ( 27), extracting T S(ρ 0 τ s ) amount of work necessarily leaves the system in the state τ s .It follows from Eq. ( 26) that, even in the presence of n copies of the bath, the maximal work extractable from all the systems altogether is upper-bounded by a fixed quantity, T S(ρ 0 τ s ), meaning that, if the net extracted work in a single cycle would be positive, then, after sufficiently many cycles, one would surpass W max , which is impossible.For our cycle, this indeed means that 0, which we now derived as a consequence of the second law.Interestingly, Eq. ( 26) also implies that, since at the beginning of each cycle the Hamiltonian is H 0 and the state is ρ p ⊗ τ B (or ρ 0 ⊗ τ B for the first cycle), the ergotropy of the system-plus-bath is T S(ρ p τ s ) > 0 (if the bath is large).However, as we noted above, spending that resource in one cycle would leave the system in the state τ s , thereby trivializing all subsequent cycles.
Note that the argument about no positive net extracted work applies to an arbitrary cycle, not just the one defined by Eqs.(3a)-(3d).More formally, any [attach]-[operate]-[detach]-type cycle (after the first one) of a single-bath machine can be thought of as a unitary transformation ρ p ⊗ τ B → = U ρ p ⊗ τ B U † such that Tr B = ρ p (U is the unitary evolution operator generated by the cyclic variation of the Hamiltonian, which depends on the particularities of the cycle).Now, it immediately follows from I (s : B) 0 that S(Tr s ) S(τ B ). Therefore, for the invested ("dissipated") work we have: , which, using the H B = −T ln τ B − T ln Z B identity, we rewrite as This relation of course applies to our cycle as well.The merit of Eq. ( 21) for canonical baths is in the nuance that, in general, Eq. ( 6) applies only to small subsystems of the whole, which means that the final state of the bath, Tr s , does not coincide with Tr s τ T , and accessing the state of the bath is usually an intractable problem.Whereas Eq. ( 21) does not require knowledge of Tr s .The content of Eq. ( 21) is much more nontrivial when the baths are microcanonical.In that case, the baths themselves are active, in the sense that the ergotropy of the bath (or the system-plus-bath) can be ∝ √ N [18].Moreover, it is possible to extract ∝ √ N work from them in a (global) process that is cyclic in terms of both the Hamiltonian and the state of the system, in other words, there exist [attach]-[operate]-[detach] cycles which extract O √ N amount of work on each run.In this context, Eq. ( 22) [and Eq. ( 29) below, in the worst-case scenario for the first cycle] pose very strong restrictions on the operation of the device.These restrictions are the price we pay for designing our cycle in such a way that the maintenance of the active state comes for free.
We observe, in passing, that Eq. ( 26) also allows us to lower-bound w (1) .Indeed, since one can extract at most T S(ρ 0 τ s ) net amount of work in the first cycle, Importantly, this bound holds also when the bath is microcanonical, since the expression for w (1) , Eq. ( 15), is the same for canonical and microcanonical baths.
Remarkably, in the context of Eq. ( 25), Eq. ( 21) provides yet another nontrivial insight: on the one hand, the device needs nonzero coupling in order to function; on the other hand, nonzero coupling inevitably leads to dissipation.This conflict is an analog of the power-efficiency trade-off in ordinary thermal machines (see, e.g., Refs.[60][61][62]).
Finally, let us comment on the necessity of attaching the system to a new bath at the beginning of each cycle.If one would have to literally keep many copies of the bath in store to operate the device, then the whole setup would have virtually no practical significance.However, we expect a generic many-body, short-range-interacting bath which couples to the system locally will, upon the start of each cycle, appear to the system as if it were new.The idea is that, as we argued above, for the bath to be thermalizing, it needs to have exponentially decaying correlations, which generically means that there is nonzero Lieb-Robinson velocity in the bath [48][49][50].Therefore, once the system is detached from the interaction site of the bath, the area around the interaction site, perturbed by the interaction with the system, will diffuse away from the site, carrying the system-bath correlations along with it.Thus, by the time the next interaction session starts, the interaction site will be only slightly perturbed, as if the bath were fresh.An example explicitly illustrating such a behavior, on the specific model of the bath being a linear chain of harmonic oscillators with nearest-neighbor interactions, was reported in Ref. [57].Of course, when dealing with finite-size baths, after a certain amount of cycles, the perturbations will travel back to the interaction site.The rigorous formalization and justification of the described picture and the study of finite-size-bath effects delineate a motivating class of problems for future studies, but are beyond the scope of the present work.

III. CALDEIRA-LEGGETT MODEL FOR THE DEVICE
We will now illustrate the theory developed above on the example of a harmonic oscillator (the system) linearly coupled to a bath consisting of independent harmonic oscillators and starting in a Gibbs state.This system-bath model is widely used to study quantum Brownian motion and is known as the Caldeira-Leggett (CL) model [22].Its total Hamiltonian, H (CL)  T , is thus defined through where the renormalization frequency of the system, ω R is given by and ensures that the Hamiltonian is bounded from below.
In the continuous limit, namely, when the bath frequencies, ω k , are very close to each other (ω k+1 − ω k ω 0 ) and range from near-zero values to those significantly larger than ω 0 (see, e.g., Ref. [53] for a careful treatment of the discretecontinuous transition), there exists a unique steady state to which the oscillator evolves [41].Moreover, this state is Gaussian, since it is given by Eq. ( 6) [41] and τ (CL)   T is a Gaussian state (as H (CL)   T is quadratic).This means that the oscillator's steady state can be fully described by the first moments, q ∞ , p ∞ , and second moments, where x = (q, p) and is the infinite-time average and σ ∞ i j comprise the so-called covariance matrix of the oscillator in the steady state.The covariance matrix can be defined for any state [see Eq. (E2) in Appendix E], and it fully determines the state if the state is Gaussian.
Writing the Heisenberg equations, and introducing the so-called spectral density, one can show that (see, e.g., Ref. [22] or Appendix D), as long as where i is either 1 or 2 and ω −ω (the symbol P in front of an integral means the Cauchy principal value; see Appendix D for details).
In this paper (except for Appendix F), we will work with so-called Ohmic spectral densities, which are characterized by a linear scaling of J (ω) with respect to ω, for small ω's [22].Also, since a given oscillator with frequency ω 0 cannot couple to bath modes with much higher frequencies, J (ω) must be a decaying function for ω ω 0 .With these in mind, we write the spectral density as where ω c is the cutoff frequency and f (z) is a well-behaved dimensionless "cutoff" function that decays for z > 1 (i.e., when ω is above the cutoff frequency ω c ) and f (0) > 0. In this notation, taking into account Eqs. ( 34) and (31), therefore, for ω R to be finite (so that the model has a physical meaning), f (z) should decay fast enough for the second integral in Eq. ( 37) to be convergent.Here, γ is a dimensionless constant defining the coupling strength.The most common choices for the cutoff function in the literature are the Lorentz-Drude, and exponential, f (exp) (z) = π e −z , cutoff functions [22].Equation (35) clearly shows that ρ ∞ (CL) , the thermal state of a free (i.e., not interacting with any other system) harmonic oscillator with frequency ω 0 is characterized by and only in the zero-coupling limit do σ ∞ 11 and σ ∞ 22 coincide with σ (free) 11 and σ (free) We emphasize that Eq. ( 6), i.e., ρ ∞ (CL) T , holds in all coupling regimes.

A. Energetics of the single-oscillator device
For the Caldeira-Leggett model described above, we can derive explicit expressions for all the relevant quantities: the ergotropy, E (CL) [Eq.( 11c Starting with the ergotropy, we show in Appendix E 3 that the ergotropy of the state ρ ∞ (CL) s [the covariance matrix of which is σ ∞ in Eq. ( 35)], with respect to H (CL) s , is given by and the extraction of the ergotropy leaves the system in the Gaussian state ρ (CL) p characterized by the covariance matrix which is in fact a thermal state at inverse temperature ).From Eq. ( 40) we immediately see that the ergotropy vanishes if and only if the steady state of the system, ρ ∞ (CL) s , is characterized by energy equipartition, namely, the average kinetic and potential energies are equal: As we will discuss below, this occurs in both weak-coupling and high-temperature limits.In the following, we will assume that the device has already performed the first cycle, so that the state of the total system at the beginning of the cycle (stroke 1 ) is where we emphasize again that ρ (CL) p is the "exhausted" state, characterized by the covariance matrix in Eq. (41).With this premise, invoking the fact that Tr [q k τ (CL) B ] = 0, ∀k, and taking into account the definitions ( 30) and ( 32 where, in the second equality, we used Eq. ( 41).[cf.Eq. (11b)], we have where we have used Eq. ( 33) to get rid of the q k g k q k term in H (CL)

22
+ ω 0 T 2π γ , where T and T are dimensionless functions of ω 0 , ω c , and T [see Eqs.(G19) and (G22)].Using these in Eqs. ( 40) and ( 44), we find that W (CL)   c:d ∝ ω 0 γ and see Appendix G 4 for details (as well as a discussion on the low-temperature limit).The message of Eq. ( 45) is that the device is essentially useless in the weak coupling limit: not only the output is small-which is indeed expected in this regime-but also the efficiency is vanishing.Interestingly, the ultrastrong-coupling limit (γ → ∞) is not much better: the efficiency also goes to zero, this time, , which, plugged into Eqs.( 40) and ( 44), yield W (CL)   c:d ∝ ω c γ and With these observations, we expect that the efficiency, as a function of γ , is maximized at some intermediate value of γ , which is in fact what we observe in Fig. 1(a), where, for the Lorentz-Drude cutoff function [Eq.(38)], the numerically calculated η (CL) is plotted against γ , for three different values of the temperature, T .Figure 1(b) shows the dependence of E (CL)  on γ for the same choice of parameters, and it can be seen how E (CL) changes its convex (γ 2 ) behavior, at small γ 's, to concave behavior ( √ γ ), at large γ 's.Plausibly assuming that the thermalization (hence, charging) time, t c , monotonically decreases with γ (see the corresponding discussion in Sec.II) for all values of γ , and viewing the device as a thermal machine, we notice that the power per cycle, E (CL) /t c , increases with γ [see Eq. ( 46)].Therefore, the assumption holding, there is a certain power-efficiency tradeoff: for sufficiently large γ 's, more power means less efficiency.
An interesting observation from Fig. 1 is that, after its peak, η (CL) decays rather slowly with the increase of γ , which means that the device can maintain a close-to-maximum efficiency while producing a significantly larger amount of output work (ergotropy) than that at maximum efficiency.Indeed, e.g., for T = 0.1, the maximum of η (CL) (≈6.5%) is reached at γ max ≈ 3.8, with E (CL) (γ max ) ≈ 0.7, whereas at γ = 14.4, the device still operates at ≈80% of the maximal efficiency, while outputting 3.2 times as much ergotropy compared to that at γ max .
We also observe from Fig. 1 that both η (CL) and E (CL) are decreasing as the temperature increases.It can also be  40) and ( 44), these yield and therefore, η We see that the device's figures of merit, E (CL) and η (CL) , decay with T very quickly, which means that low temperatures are essential for the efficient operation of the device.
We can see in Fig. 1 that the device operates with a rather small efficiency.Looking to find high-efficiency regimes, we turn to ω c -the only parameter we have not explored yet.Generically, it is assumed that ω c is large enough to be greater than ω 0 , and it is the choice made in Fig. 1: ω 0 = 2 and ω c = 4.However, it turns out that the maximal-over all the parameters-efficiency the Caldeira-Legget model can allow, with the condition that ω c ω 0 , is ≈10.19%, which is achieved at T = 0, γ ≈ 5.94, and ω c = ω 0 , and does not depend on the value of ω 0 , as long as it is >0.Indeed, we already knew that η (CL) is a monotonically decreasing function of T , so T = 0 follows trivially.With T = 0, the only dimensional parameters are ω 0 and ω c , therefore, the dimensionless η (CL) FIG. 2. Figures of merit of the device as functions of the cutoff frequency.On the main plot, efficiency, η (CL) is shown as a function of ω c , for different values of the coupling constant, γ .Observe that the higher the desired efficiency, the sharper the peak and the smaller the ω c have to be.The inset, where E (CL) and W (CL)  c:d are plotted against ω c , for γ = 15, is to illustrate that, in contrast to the efficiency, E (CL)  and W (CL)   c:d are monotonically increasing, respectively, concave and convex functions of ω c .The cutoff function is of Lorentz-Drude form and ω 0 = 2, T = 0.1.can depend only on ω 0 = ω 0 /ω c .We find numerically that, for sufficiently large values of γ , η (CL) is an increasing function of ω 0 , for ω 0 1; this fact is illustrated in Fig. 2. The only free parameter left now is γ , and we find the above-mentioned optimal values of γ and η (CL) also numerically.To sum up, if we want efficiencies higher than 10%, we need to consider ω c < ω 0 .Such situations may occur when the bath is a harmonic (Rubin) chain with nearest-neighbour interactions [53], and the individual frequencies of the oscillators in the chain are smaller than ω 0 .In Fig. 2, it is shown that, for sufficiently high γ 's, the efficiency does go above the 10% value.The nonmonotonic behavior of η (CL) is due to the fact that E (CL) is concave, whereas W (CL)   c:d is convex, although both E (CL) are W (CL)   c:d monotonic in ω c (see the inset of Fig. 2).Importantly, high efficiencies are obtained only at small values of ω c , which means that the output work also has to be small.On the other hand, with the increase of γ , the peak around the maximum becomes increasingly sharper, which means that the cutoff frequency of the bath needs to be fine-tuned to achieve higher efficiencies.The efficiency can be further increased by taking larger ω 0 's and γ 's; we numerically found that the highest efficiency possible within this model is 50%, which is achieved in the ω 0 → ∞ and γ → ∞ limit.

B. Several identical oscillators attached to a common bath
As we saw in the preceding subsection, for the device to produce a significant output ergotropy with a reasonably high efficiency, large values of the coupling constant are necessary, which may be challenging to achieve in practice.Here, we discuss a collective enhancement effect appearing when the system is comprised by n copies of the oscillator, which are simultaneously coupled to a common bath.Below, we will show that this system is essentially equivalent to a single oscillator coupled to the bath, but with a rescaled coupling constant: nγ .Indeed, the n-oscillator Caldeira-Leggett Hamiltonian with a shared bath consists of the system Hamiltonian, n a=1 , and the bath Hamiltonian, H (CL)

B
, comprising together the bare Hamiltonian, and the interaction term: By introducing we can rewrite the total Hamiltonian as where with I n being the n × n identity matrix and the symbol ⊕ denoting the direct sum.Note that the indices a and b will, in this subsection, distinguish the system degrees of freedom from those of the bath, labeled by k.Now, we introduce a change in the system's variables, according to Q a = b ℵ ab q b and P a = b ℵ ab p b , with ℵ being an orthogonal matrix (ℵℵ T = I n ) such that Q 1 = a ℵ 1a q a coincides with that defined in Eq. ( 47).(Note that this condition fixes only the first row of ℵ.) Given the special structure of M [Eq.(48)], this transformation leaves M unchanged.Therefore, in terms of the new system variables, the total system-bath Hamiltonian will be which means that the original total system is equivalent to a collection of n − 1 free oscillators and a single oscillator coupled to a bath.As per the initial state of the system, given the symmetry of the problem, we choose it to be a product of identical single-oscillator states, and the initial system-bath state is, as usual, given by Eq. ( 2).Moreover, we will assume the single-oscillator states to be purely quadratic Gaussian, so that where A, B, and C are real, subject only to the condition that the operator Aq 2 a + Bp 2 a + 1 2 C{q a , p a } > 0. Noting that, due to the orthogonality of ℵ, a q 2 a = a Q 2 a , a p 2 a = a P 2 a , and a {q a , p a } = a {Q a , P a }, the initial state looks the same from the perspective of the new variables: Let us note at this point that, although all states of the form of Eq. ( 50) have identical spectra, they all live in different subspaces of the system's Hilbert space.Summing up Eqs. ( 49) and ( 51), the model consists of n identical oscillators with frequency ω 0 , n − 1 (numbers 2 to n) of which are uncorrelated and uncoupled from the first one and from the bath, whereas the first one evolves as a standard quantum Brownian particle with rescaled coupling: The states of the uncoupled oscillators evolve under the influence of their internal Hamiltonian [and therefore (i) they do not relax and (ii) their ergotropy is constant in time] and remain uncorrelated to the state of the first particle (which thermalizes with the bath).
Moreover, as we show in Appendix E, the maximal work cyclic (in Hamiltonian) Gaussian operations can extract from a (multimode) Gaussian state-the "Gaussian ergotropy"-is given by where {s ↑ a } are the symplectic eigenvalues of the covariance matrix, σ [defined in the multimode case identically to the single-mode case given in Eq. ( 32)], taken in increasing order, and {m ↓ a } are those of M, taken in decreasing order [see Appendix E for detailed definitions and a proof of Eq. ( 52)].Now, since all n symplectic eigenvalues of M are equal to each other (all are ω 0 ), the ergotropy of the n-oscillator system is equal to the sum of individual ergotropies.Therefore, our protocol [Eqs.(3a)-(3e)] will independently process the first (Brownian) particle and the rest of the n − 1 particles.Thus, the latter will have no role in the energetics at all, because we will extract all the erogtropy from them during the first cycle-which will be nonzero only if ‫(ג‬A, B, C, Q a , P a ) is active with respect to 1  2 P 2 a + 1 2 ω 2 0 Q 2 a -and for the rest of the time, they will just remain in their passive states.The energetics for the Brownian particle (Q 1 , P 1 ), on the other hand, will be identically as described in Sec.III A, with the only difference that, since G k = g k √ n and J (ω) is quadratic in g k [Eq.(34)], the coupling constant is now nγ .

IV. DISCUSSION
In this paper, we investigated the arguably simplest model of charging a battery: a system in strong contact with a bath, jointly evolving towards global thermal equilibrium (or a nonthermal state appearing to be thermal when viewed locally).The idea is based on the basic observation that the reduced state of a subsystem of a globally thermal system is not generally thermal and may thus harbor extractable work when disconnected from the global system.This setup presents a number of benefits, not simultaneously met in any other battery-charger setup, such as no need of fine control over the preparation of the charged state and no necessity to exert any effort to maintain the battery in that state.Indeed, the existing setups either require isolating the system or engineering a special system-bath evolution [9,11,12] or maintaining fragile internal symmetries [13,14] or active external stabilization [16] to preserve the charged state.
We studied this setup in full generality, revealing fundamental limitations on the process, encapsulated in a nonnegative entropy production [Eq.( 21)], quantifying the amount of work that needs to be dissipated in order for the device to function.A key peculiarity of the charging cycle is that, nonweak coupling, being responsible for preparing the charged state of the battery, is also an inevitable source of dissipation.Indeed, at least the disconnecting stroke must be fast, because, otherwise, by the end of the disconnecting stroke, the system will be in thermal equilibrium with the bath while being only weakly coupled to it, i.e., it will be in a Gibbs state, which is passive.Therefore, the device is expected to function optimally in the moderate-to-large coupling regime, which is what we saw on the example of Caldeira-Leggett model in Sec.III A.
Rephrasing the previous paragraph, the autonomy and robustness of our device come at the cost of limited efficiency.However, it is important to note that η is merely an upper bound to the "de facto efficiency," which we define to be one that compares E to the actual energy spent on the charging cycle, including the energy spent on control and stabilization.To illustrate this point, observe that one can simply take a system in the ground state and call it a depleted battery; then, one can use a unitary transformation to rotate that ground state to the highest energy eigenstate of that system, thereby charging the battery.The "formal efficiency" of that charging cycle is 1-no energy is dissipated and all energy transferred to the battery is unitarily extractable [63].Immersed into reality, this simple picture breaks down by the necessities of isolating the battery and performing tailored unitary operations on it.All this requires an intervention of macroscopic laboratory equipment that consumes macroscopic amounts of energy.Moreover, keeping an excited system from spontaneously emitting a photon (i.e., leaking the charge) also requires complicated equipment, thus also requiring macroscopic energy expenditures.With all these taken into account, the de facto efficiency of that simple battery will essentially be zero.Although our design is not completely autonomous in that it is sensitive to the bath's temperature and certain, albeit little, external control is required, it is free from two major sources of macroscopic energy expenditure described above.Therefore, the de facto efficiency of our design is arguably higher than that of other existing designs, which are at best free from only one of such sources of energy expenditure and therefore suffer close-to-zero de facto efficiency.Moreover, to the best of our knowledge, except for this work, the thermodynamic efficiency of charging and discharging process has been addressed only in Ref. [9], albeit in a slightly different setting.
In our cycle, the connection and disconnection steps are instantaneous, therefore, the minimum time needed to run the cycle is determined by the relaxation time, which in turn depends on the strength of the coupling.As we briefly mentioned in Sec.II B, other time-dependent protocols for connecting and disconnecting the system could be considered, exploring whether (and ultimately to what extent) the efficiency can increase and whether that occurs at the expense of the time it takes to run the cycle.
Another aspect of the device is that the "battery" system has to be small.Indeed, the athermality of a macroscopic system coupled to a bath will generically be a boundary effect tending to thermalize away once the system is detached from the bath.Moreover, generically, the effect of systembath correlations due to interaction are more pronounced at low temperatures.Therefore, although the working principle of our device is not inherently quantum, it is better suited for the quantum regime.This aspect is well illustrated by the Caldeira-Leggett model, for which the machine is most efficient in the low-temperature regime, where its quantum features are most prominent.In the classical (high temperature) regime, the system exhibits energy equipartition, independently from the details of the bath and the strength of the coupling, so that the ergotropy vanishes.We note, however, that in this work we forewent studying quantum effects such as system-bath entanglement and coherence in the system's state, leaving the study of these important questions to the future.
In a simplistic scenario of n noninteracting oscillators coupled to a common bath, in Sec.III B, we showed that enhancement of both the output and efficiency takes place in the Caldeira-Leggett model.Situations where collective effects are more pronounced, e.g., when the "working medium" is critical [64,65], constitute a promising direction of future research.
Lastly, we envision that our battery-charger setup can be experimentally realized on certain quantum optical platforms where strong and ultrastrong coupling regimes are reachable [66,67], such as, e.g., cavity QED [68,69].Another possible direction to look for practical realizations of our protocol, is using chemical bonds, which are a good example of naturally occurring strongly coupled microscopic systems.For example, the cycle can be realized through creating chemical bonds, which can store the energy for a long time, and then breaking them and using the athermal states of the components to extract work.the initial moment of time t = t in , one drives the system Hamiltonian according to some time-dependent protocol H (t ) in such a way that, at the end of the process (the final moment of time t = t fin ) the system Hamiltonian is back to its original value: H (t in ) = H (t fin ) = H.Such processes are called cyclic Hamiltonian processes, and their fundamental importance in thermodynamics is that, at the end of the process, we end up with the same system as we had at the beginning of the process.The Kelvin-Planck formulation of the second law refers to these very processes.Indeed, if one considers processes at the end of which the Hamiltonian is allowed to differ from the initial Hamiltonian, then extracting work from a thermal system becomes trivial: imagine an initially thermal gas in a chamber expanding adiabatically and pushing against a piston-this process extracts work from a system in thermal equilibrium, as long as one does not require the piston to be back at its original position at the end of the process.
The time evolution of the state under the influence of the time-dependent Hamiltonian is unitary: for ∀t ∈ [t in , t fin ], where the unitary evolution operator is standardly written as the symbol T signifying chronological ordering.The inverse is also true: any unitary evolution can be generated by a time-dependent Hamiltonian process.Indeed, given a unitary U (t fin , t in ), we can generate it by the cyclic Hamiltonian process where we, at the moment of time t in abruptly change the Hamiltonian from H to i t fin −t in ln U (t fin , t in ), let it run until the moment of time t fin , and then abruptly change the Hamiltonian back to H.In view of this, the ergotropy of a system in state ρ with respect to Hamiltonian H can be defined as therefore, E is sometimes also called unitarily extractable work.
The unitary evolution operator delivering the minimum in Eq. (A3), (A4) thus takes the system to a state from which no more work can be extracted (by means of cyclic Hamiltonian processes).The latter, ρ p = U E ρU † E , is called passive [7,8], and it can be shown that, in the eigenbasis of H with the basis elements chosen such that the eigenvalues of H are in the increasing order: where r ↓ k are the eigenvalues of ρ in the decreasing order (r Obviously, since ρ p is diagonal in the eigenbasis of H, [ρ p , H] = 0. Also note that, since all positive-temperature Gibbs states over H are already in the form of Eq. (A5), they are all passive; however, not all passive states are Gibbs states (see the discussion about complete passivity in Refs.[7,8]).Lastly, let us note that the ergotropy of a system is invariant under the internal dynamics.Indeed, Tr [He −iHt ρe iHt ] = Tr[Hρ] and, since e −iHt

APPENDIX B: THERMALIZATION IN NONCRITICAL MANY-BODY SYSTEMS WITH SHORT-RANGE INTERACTIONS
Generally, the results on thermalization are obtained in two steps.One first proves that the system equilibrates, i.e., ν S ( ) := lim To establish that, one uses the general bounds on equilibration obtained in Refs.[29,31,34]; most appropriately in our situation [34], where d S is the Hilbert-space dimension of the subsystem S, G D is the highest "gap degeneracy" of H T (gap degeneracy is the number of times a given gap is repeated in the spectrum, and G D is the largest such number), and d eff , the "effective dimension," is where P a are the eigenprojectors of H T (H T = a E a P a , where E a are H T 's eigenvalues).When, as in our setting, d S is a fixed finite number, the thermalization is guaranteed as long as G D d eff is 1 and tends to zero with N → ∞.Generically, for canonical and microcanonical R B 's, one expects d eff to be exponential in N in view of the exponential (in N) density of eigenvalues of H T and H B [38].However, that is not always the case; instead, in Ref. [35], it was shown that, whenever H T is k-local and has exponentially decaying correlations, which, since in our case = ρ 0 ⊗ R B , is guaranteed as long as R B has exponentially decaying correlations, where d is the spatial dimension of the lattice, meaning that equilibration is guaranteed as long as G D does not scale with N (or scales slower than √ N), which is expected to not hold only in extremely exotic cases as it holds in all known models and, in general, Hamiltonians with degenerate gaps are of zero measure in the space of all Hamiltonians and any given degeneracy can be lifted by an infinitesimal perturbation.See also Ref. [31] for a discussion on equilibration timescales.
Another important bound proven in Ref. [21], and adapted in Ref. [35], is the following.Say, τ is some state with exponentially decaying correlations on the lattice and 0 < α < 1.Then, for any state ρ satisfying where E m denotes the arithmetic mean over all subsystems S m of diameter m lattice units.For this bound to be valid, m must be • N d+1+α (d+1)(d+2) , which is the case in our analysis as we will always work with subsystems that do not scale with N, i.e., m will always be a finite number.Note that the bound in Eq. (B5) holds on average, whereas in this paper we need it to hold only for all subsystems in s ∪ supp(H I ).Importantly, we do not require Eq. (B5) to be valid for all subsystems of T, so we do not need to require translation invariance of T. However, H T must satisfy certain transport properties such as nonzero Lieb-Robinson velocity [44].Indeed, in systems with inhibited energy/information transport, e.g., those that are many-body localized, there will be subsystems which "remember" their initial states (see Ref. [39] for a discussion); here, we require that no localization phenomena occur on (or near) s ∪ supp(H I ).
Combining Eqs.(B1), (B3), and (B5), and using the triangle inequality for the trace norm [33], one obtains [35] . Since is obtained from by erasing some of its nondiagonal elements in H T 's eigenbasis, S( ) S( ).Also, by 's definition, it holds that Tr(H T ) = Tr(H T ).Therefore, S( τ T ) S( τ T ), meaning that the bound in Eq. (B6) holds as long as Now, when the initial state of the bath is τ B , for any ρ 0 , we can write where the first inequality is due to the Peierls-Bogoliubov inequality [71], the second inequality is by the very definition of the trace norm [33], and the last equality is due to the fact that H I is at most k-local (because H T is).Hence, Eq. (B5) is satisfied with α = 0, which, in view of Eq. (B4), means that ν S (τ T ) goes to zero polynomially with N. In other words, when the bath starts in a canonical state (R B = τ B ), thermalization as in Eq. ( 6) takes place.Using the results in Refs.[21,51] about the equivalence of canonical and microcanonical ensembles, and combining them with results in Ref. [35] mentioned above, we will now show that, with the same requirements on H T as above, thermalization as in Eq. ( 6) takes place also when the bath starts in a microcanonical state, μ(E , ), as defined in Eq. (10).To ensure that τ B and μ For such a choice, as is proven in Refs.[21,51], where S is an arbitrary subsystem of B with a diameter O N 1 d (d+1) .In order to prove that the system equilibrates, we need to ensure that d eff is large.To see that it is, we note that, in view of local indistinguishability of μ B (E , ) and τ B [Eq. (B13)], the correlators on μ B (E , ) and τ B will coincide up to an correction, meaning that, up to distances O(ln N ), correlations in μ B (E , ) are guaranteed to decay exponentially, which is sufficient to prove the bound in Eq. (B3) [35].In systems where the density of eigenstates of the bath increases exponentially with N, d eff is obviously guaranteed to also be large (keeping in mind that H I does not scale with N, see also the discussion in Ref. [31]).
Again, using the Peierls-Bogoliubov inequality in the last term and the definition of trace norm as we did in Eqs.(B10)-(B11), we get that Finally, invoking the following inequality from Ref. [21] (Lemma 7): we obtain that S( mc τ T ) O ln 2d N , and therefore it satisfies Eq. (B4).In turn, this means that Eq. (B6) again holds, implying that thermalization takes place also when the bath starts in a microcanonical state.

The case when supp(H I ) is not finite
We have just shown that, for any H I , as long as supp(H I ) is finite, Eq. ( 6) generally holds up to a correction O(N − ), with some > 0 which is to be read from Eq. (B6) [although it is possible to be more specific, it is sufficient to note that, when G D scales at most as N y , with some 0 y < 1/2, O(N − ), with an arbitrary 0 < < min { 1 4 − y 2 , 1−α 2d+4 }, will be an upper bound for the corrections to Eq. ( 6)].
Let us now turn to the case when supp(H I ) is not finite, e.g., when it scales with N. In this case, the O(N − ) corrections may potentially sum up into something nonnegligible.Therefore, we consider only such H I 's that can be decomposed into a sum of local terms: H I = κ γ κ V κ , where γ κ 0 are the "coupling constants" and, for all κ's, the diameter of supp(V κ ) is < C, with C being some finite number not depending on κ.Moreover, γ κ are adjusted so that, for any κ, max{| Tr(V κ )|, | Tr(V κ τ T )|} K, where K > 0 is a finite constant independent on κ.Then, for each term, the substitution of R B by τ B and ∞ by τ T in Eqs.(11a) and (11b) introduces an f (γ κ )O(N − ) correction (here f is some function; see next).Therefore, to guarantee that Eq. ( 12) is correct, we require the series κ f (γ κ ) to scale slower than N .Note that this case formalizes the situation when the system is attached to a bosonic bath, as in the Caldeira-Leggett and spin-boson models, where the accumulation of corrections is forestalled by the introduction of a cutoff frequency (see Sec. III for definitions).There, V κ 's correspond to q ⊗ q k 's and γ κ 's to g k 's.Due to the fact that Tr(q k τ B ) = 0, any nonzero quantity in Tr(H I τ T ) will be due to system-bath correlations in τ T stemming from g k = 0 [72], and, given that H I is itself ∝g k , any such quantity will be ∝g 2 k , meaning that f , which, for the exponential cutoff, is finite and, for the Lorentz-Drude cutoff, is ∝ ln N, which is slow enough to guarantee the correctness of Eq. ( 12).

APPENDIX C: THE DERIVATION OF EQ. (21)
Here, we show how to get from Eq. ( 18) to Eq. ( 21).To do so, we start with Eq. ( 18), and adding and subtracting Tr[H T τ T ] to it, we rewrite it as where S(ρ||τ ) = Tr[ρ(ln ρ − ln τ )] is the relative entropy [33].Using the identities H T = −T ln τ T − T ln Z T and H B = −T ln τ B − T ln Z B , we can further transform Eq. (C3) into where, adding and subtracting T Tr[τ T ln τ T ], we obtain which, by noting Eq. ( 20) and adding and subtracting we finally bring to which is exactly the desired Eq. ( 21).

APPENDIX D: STEADY STATE OF THE HARMONIC CALDEIRA-LEGGETT MODEL
The harmonic Caldeira-Leggett is described by the Hamiltonian given in Eq. ( 30), and the Heisenberg equations of motion are given by Eq. (33).The solution of the set of equations for q k 's in Eq. ( 33) can be written as which, substituted into the first equation in Eq. ( 33), yields where Here, noting that t 0 is the point in time when the description starts, and so it is the point in time at which the Heisenbergpicture operators are initialized by their Shchrödinger-picture originals, we substituted q k (t 0 ) = q k and p k (t 0 ) = p k .Moreover, we also note that, since the choice of t 0 is arbitrary and we are primarily interested in the long-time behavior of the system, for convenience, while always keeping t 0 −|t| finite, we will take the limit t 0 → −∞ at relevant places.
Using the definition of the spectral density, Eq. ( 34), we can write ) Switching to the Fourier-transformed functions [q(ω) = ∞ −∞ dte iωt q(t )] in Eq. (D1), and keeping in mind that t 0 → −∞, we obtain and hence where and, through the Sokhotski-Plemelj theorem [73] (a.k.a., Kramers-Kronig relations), where, and henceforth, it is understood that J (ω) is extended to negative arguments as an odd function, in accordance with Eq. (D9).The sign P signifies that one takes the Cauchy principal value of the integral, namely, P Note that, as long as t 0 → −∞, the only solution of the free equation (D1), namely, one with ξ (t ) put to zero, is q(t ) ≡ 0, which is a consequence of α(ω) = 0 for any value of ω, provided that J (ω) = 0 only for ω = 0 and ω 0 > 0 (which is what we assume throughout this paper).Therefore, the unique solution of Eq. (D1) for any initial conditions is Note that, since t 0 → −∞, these solutions are the steadystate solutions (the time-dependence is simply because [H (CL) s , q] = 0).Now, since the initial state of the bath is a Gibbs state, and therefore Tr[ and, analogously, Furthermore, taking into account that the initial state of the bath is a thermal state, and hence, it is easy to show that the steady-state covariance matrix σ ∞ of the oscillator [see Eq. (E2)] is and σ ∞ 12 = σ ∞ 21 := 1 2 {q, p} t→∞ = 0. Lastly, let us find qq ∞ since we need it in Sec.III A. Calculating q using Eq.(D11), multiplying it by q(t ) from the same equation, and then using Eq.(D14), we find that

APPENDIX E: GAUSSIAN ERGOTROPY
Here we show how to calculate the maximal work cyclic Gaussian Hamiltonian processes can extract from a Gaussian state.We give this maximal amount of work the natural name of Gaussian ergotropy as all the involved states, Hamiltonians, and (therefore) unitaries are Gaussian.

Setting notation
Suppose we are given a d-mode bosonic system described by d pairs of coordinates and momenta (q i , p i ).All Hamiltonians are going to be quadratic, so, composing the column vector x = (q 1 , p 1 , • • • , q d , p d ) T , we introduce the symmetric matrix M via Note that M has to be positive-semidefinite because otherwise H will have a spectrum unbounded from below, which is unphysical.
Let us furthermore define the covariance matrix, σ , as where ρ is the state of the system [cf.Eq. ( 32)].Note that σ is a symmetric, positive-definite matrix.Thus, for the average energy, we obtain In order to proceed, we need to keep in mind the following two facts: (i) Any Gaussian unitary transformation of the state taking ρ to ρ = U ρU † amounts to the covariance matrix where U is a real symplectic matrix; symplectic is any matrix satisfying U J T U = J, where with is the symplectic identity; the real symplectic matrices constitute a group usually denoted by Sp(2d, R).The statement in Eq. (E4) works in both directions: any symplectic transformation of σ can be generated by a Gaussian unitary evolution of ρ.
(ii) The Williamson's theorem [74] holds.Namely, any 2d × 2d symmetric matrix σ 0 has d nonnegative symplectic eigenvalues s 1 , • • • , s d and there exits a symplectic matrix where In Appendix E 4 a below we show how to obtain the Williamson decomposition from canonical Schur decomposition.The latter is built into most major packages (such as Mathematica, SciPy, ALGLIB, etc.), thus the method provides a ready protocol for numerically computing the Williamson decomposition.Note that the presented method is a widely used one and is by no means an invention of ours.

The ergotropy
Gaussian ergotropy is the maximal amount of work cyclic Gaussian Hamiltonian processes can extract from a Gaussian state.Keeping in mind that any Gaussian Hamiltonian process generates a Gaussian unitary evolution operator and, vice versa, any Gaussian unitary evolution operator can be generated by a Gaussian Hamiltonian process, we can write the Gaussian ergotropy, G, as where the minimization is carried over the set of all Gaussian unitary operators.The state U ρU † delivering the minimum is called Gaussian-passive and such states were fully characterized in Ref. [75].Now, recalling Eqs.(E3), (E4), and (E7), we can write 2 Tr s ↑ T m ↓ , where m ↓ is the symplectically diagonalized M (with correspondingly ordered elements) (namely, denoting the symplectic eigenvalues of M by m i , σ is a real symplectic matrix.Thus, we can rewrite Eq. (E9) as whence it immediately follows that min Tr[ where the minimization is carried over the whole group of real symplectic matrices, i.e., ∈ Sp(2d, R).Note that there is no a priori guarantee that a single will minimize all S k ( )'s; however, invoking Lemma 1 of Ref. [76], stating that we see that = I 2d simultaneously minimizes all S k ( )'s, thereby showing that We have thus just proven the main result of this section, namely, that and the maximum is delivered by the unitary U that generates the real symplectic transformation, To the best of our knowledge, Eq. (E14) in this general form has never been reported in the literature.An analog of Eq. (E14) for free-fermionic systems is reported in Ref. [77].
It is a trivial consequence of Eq. (E14) that, when the system consists of d noninteracting modes, the Gaussian operation extracting maximal work is the one that evolves the initial state into s ↑ , if M = m ↓ , or s ↓ , if M = m ↑ .Note that, although the maximal work is always delivered by the symplectic transformation in Eq. (E15), in some cases, it will not be unique.Indeed, as can be observed by inspecting Eq. (E11), when some of the normal frequencies coincide, the symplectic transformation of the state delivering the minimal final energy might not be unique.In fact, for d = 2 and m 1 = m 2 (i.e., m ↓ = m ↑ = m), it was shown in Ref. [75] that the minimum of Tr[ σ T m] is delivered by both σ T = σ ↑ and σ T = s canon , where s canon is in the "canonical" form [78,79]: Importantly, in the canonical form, the nondiagonal block matrices are generally of the form diag(c 1 , c 2 ), whereas for s canon above to deliver the minimum of energy it is necessary that c 1 = c 2 [75].We elaborate on this in Sec.E 4 b below.

Ergotropy of a single oscillator
Let us illustrate the above theory on the simple example of a single oscillator with Hamiltonian where x = q p and that starts in a Gaussian state ρ described by some σ = σ 11 σ 12 σ 12 σ 22 (note that here we explicitly take into account that σ is a symmetric matrix).Now, in order to use Eq.(E14), let us find the symplectic eigenvalues of σ and M 1 .Using the observation below Eq. (E28), we immediately find the (only one) symplectic eigenvalue of σ to be s 1 = √ det σ and that of M 1 : m 1 = ω 0 .Eq. (E14) then implies Note that the passive state in which the system ends up, as a result of the ergotropy extraction, is the state which, when H = , has a covariance matrix s ↑ = s 1 0 0 s 1 , which means that the final passive state, ρ p , is a Gibbs state at some temperature: ρ p ∝ e −β p H , where the temperature can be determined from S(ρ p ) = S(ρ).It obviously follows from this that the Gaussian ergotropy coincides with the full ergotropy (i.e., one found by optimizing over all-not only Gaussianunitary operations), as the Gibbs state has the lowest energy for a given value of entropy.We note that the Eq.(E19) can be deduced from the analysis in Ref. [75]; see also Ref. [80].
Note also that, upon bringing H back to the p 2 2 + ω 2 0 q 2 2 form, i.e., switching to the covariance matrix of the final (Gaussian-)passive state will take the form The temperature T p = β −1 p can be determined by noting that the covariance matrix elements of a free harmonic oscillator in a Gibbs state are given by Eq. ( 39), therefore, comparing with Eq. (E21), we find (E22)

Technical nuances a. Williamson decomposition from Schur decomposition
Let us describe how, given a σ , to obtain the and the s in Eq. (E7) numerically.To that end, we are going to make use of the standard Schur decomposition for skew-symmetric matrices: Say, A is a 2d × 2d real skew-symmetric matrix (i.e., A T = −A), then there exists a real orthogonal matrix O such that where a i are real, nonnegative numbers and This decomposition is a built-in function in most numerical software packages.Notice that, as Eqs.(E23) and (E24) readily suggest, ±ia i are the eigenvalues of A. Now, noticing that the matrix σ 1/2 Jσ 1/2 is skewsymmetric, let us write its Schur decomposition as From there, we immediately see that ¯ = σ −1/2 Os 1/2 is a symplectic matrix and that ¯ T σ ¯ = s.Introducing the symplectic matrix we observe that Moreover, ±is i are the eigenvalues of σ 1/2 Jσ 1/2 , or, equivalently, of σ J.

b. Canonical form versus Williamson form
Say, we have two noninteracting modes with both frequencies equal to some ω 0 (i.e., m 1 = m 2 = ω 0 ).Let us now see how where we keep c 1 and c 2 general, compares to First of all, we note that the covariance matrix is a positivesemidefinite matrix, which is equivalent to Next, it is easy to calculate the symplectic eigenvalues of s canon , thereby relating s 1 and s 2 with a, b, c 1 , and c 2 : where Now, it is straightforward to see that where, using the inequality 2 √ xy x + y on the square-root term, we obtain with an equality if and only if c 1 = c 2 , in which case one obtains s 1 + s 2 = a + b, which, in turn [cf.Eqs.(E29) and (E30)] implies E canon = E min , meaning that the canonical form with (and only with) c 1 = c 2 does indeed have minimal energy, corroborating Theorem 1 of Ref. [75].

APPENDIX G: ASYMPTOTIC EXPANSION OF THE COVARIANCE MATRIX WITH RESPECT TO THE COUPLING CONSTANT
In the subsequent subsections, we will derive the weak coupling expansion of the covariance matrix of an oscillator coupled to a Caldeira-Leggett bath and use that expansion to calculate the figures of merit of the device.We will also explore the low-temperature limit.
In the weak coupling limit, one expects the steady state of the system to be close to the thermal equilibrium state (τ (CL) s ), characterized by the covariance matrix comprised by σ (free) 11 and σ (free) 22 in Eq. (39).Below, we will see that this is indeed the case and, moreover, will find the O(γ ) correction to σ (free) 11 and σ (free) 22 for γ 1.

Asymptotic expansion of σ 11 with respect to γ 1
Performing the following change of variable in the integral in Eq. ( 35) for σ ∞ 11 , invoking Eq. ( 36), and noticing that ω 2 R − χ (ω) is linearly proportional to γ , and therefore is independent on γ , we find where, for convenience, we have introduced In Eq. (G3), the dependence of σ 11 on γ is localized in one place.Inspecting the part of the integrand without the coth, we immediately recognize a similarity to the lim representation of Dirac's δ function, and thus expect the leading term in Eq. (G3) to be ∝ 1 2ω 0 coth ω 0 2T , which is what one should indeed obtain [cf.Eq. ( 39)].However, we are going to need the higher-order terms in the expansion of σ 11 around γ = 0, and to make further progress, we note that the behavior of the integrand in Eq. (G3) behaves differently when x ∝ γ (it is finite) and when x ∝ 1 (it tends to zero with γ → 0).To isolate different behaviors, we divide the integration region as Let us first deal with I 1 .Introducing, for convenience, and further denoting we obtain Keeping in mind that F (x) and K (x) are quickly decaying functions for x > (ω c /ω 0 ) 2 , we do not concern ourselves with the large-x behavior of g(x) (which, for some generic choices of f , is a low-degree polynomial of x) and consider it small ("finite") as compared to x/γ in the asymptotic limit of γ → 0. Expanding the integrand in Eq. (G9) around γ = 0, we find x 3 , (G10) plus higher-order terms.Noting that, for k > 1, we see that the first term in Eq. (G10) scales as γ 5/3 and the second term scales as γ 7/3 , meaning that, to σ 11 , these contribute as γ 2/3 and γ 4/3 , and, since we are interested in the next to the leading order term in the expansion of σ 11 with respect to γ , we will discard the second term and focus on the first one: Now, since K (x) is a regular, analytic function, we can write Hence, We further notice that So, setting A = 1, and noticing that K (0) = K (0) and K (x) = −K (−x), we obtain By taking identical steps for I 3 , we can immediately write to finally arrive at Turning to I 2 , let us switch the integration variable to x/γ , so that and notice that γ x 1 in the whole integration interval, which allows us to Taylor-expand the integrand around γ x (while not touching x): where, for simplicity, we denoted Now, performing the integration in Eq. (G13) and denoting and keeping in mind that L 1, we obtain where we have introduced and Next, Taylor-expanding these expressions around 1 L = 0, we obtain which means that and hence, keeping in mind Eqs.(G4), (G7), and (G8), we find that where Recall that F (x), K (x), Y (x), and g(x) are defined in Eqs.(G4), (G7), (G16), and (G2).From these equations, we find

Asymptotic expansion of σ 22 with respect to γ 1
Turning to σ 22 [see Eq. ( 35)], we immediately see that its analysis is identical to that for σ 11 , only multiplied by ω 2 0 and with K (x) substituted by K Thus, from Eqs. (G18) and (G19), we read where here and in Eq. (G16), we obtain (G22)

The low temperature limit: σ for γ 1 and T ω 0
To find the simultaneous low-γ and low-T expansion of σ , we will combine our results above with the low-T results obtained in Ref. [53]: Last, since T /ω 0 1, we will write the • terms as γ 1+ζ m with some ζ m > 0; note that ζ 2 and ζ 4 will generally be different.
To illustrate the above formulas on concrete examples, it is useful to write explicitly Recalling Eqs.(G2) and (D5), we also have where, for simplicity, we have introduced For a generic Lorentz-Drude spectral density, J (L) (ω) = 2γ ω 0 ω ω 2 c ω 2 c +ω 2 , the cutoff function will be [cf.Eq. (38)] for which, it is a simple exercise to show that Taking into account that it is also easy to calculate the integrals in Eq. (G26), which brings us to Let us now use Eqs.(G18) and (G21) in Eqs. ( 40) and ( 44) to obtain the corresponding expansions of E (CL) and W (CL)   c:d with respect to γ .In order to do so, we first note that, since the high-temperature limit is already covered by Eqs.(F13) and (F14), the pertinent regime to explore is temperatures bounded from above by a constant not ω 0 .For such temperatures, coth ω 0 2T is a number close to 1, therefore, γ 1 is equivalent to γ coth ω 0 2T .Therefore, for E (CL) , we can use the is a number of the order of 1 (it is exactly 1 for Lorentz-Drude and exponential cutoff functions).Typically, ω c is of the order of ω 0 , therefore, = O(γ ).However, in certain situations, it may happen that ω c ω 0 , so that γ ω c O(ω 0 ), in which case ω R will not be small anymore.Therefore, although γ ω 2 R /ω 2 o will typically be •(γ ), we will keep the term proportional to in the expression for W (CL)  c:d .So, In the low-temperature limit, which we described in Appendix G 3 above, for the formulas below, we will absorb the O(γ ω c /ω 0 ) and O(T 4 /ω 4 0 ) terms into the •(1): where 0 and 0 are given in Eq. (G26) and their values in the specific of Lorentz-Drude cutoff are given in Eq. (G36).

APPENDIX H: THE SCALING OF THE COVARIANCE MATRIX IN THE ULTRASTRONG-COUPLING LIMIT
Let us now study σ ∞ 11 and σ ∞ 22 in the γ → ∞ limit.Switching the integration variable in Eq. (35) to we transform the expressions for σ ∞ 11 and σ ∞ 22 into where we have introduced where the function f is defined in Eq. ( 36), ω 0 is given in Eq. (G31), and g, an analog of g, is defined as and we will, in this section, work under this assumption.Finally, since Eq.(H5) tells us that f γ (ω) → 0 as γ → ∞, we recall the Dirac's δ function's representation in Eq. (G5) and use it in Eqs.(H2) to conclude that, as γ → ∞, Here, since we are interested in only zero-order terms, we can neglect the 1 γ inside the δ functions, substitute the coth's with 1 and g γ (ω) with g ∞ , according to Eq. (H6).Hence, Coming back to g ∞ , it is easy to see from Eq. (G33) that, for the Lorentz-Drude cutoff function, We immediately notice that, in the ω/ω c → ∞ limit, the second integral in Eq. ( H10) is damped at least as fast as e − ω ωc .The first integral is clearly dominated by the values at t = 1.
Therefore, switching the integration variable to z = ω ω c (1 − t ), we find that Lastly, let us mention that, since the asymptotic analysis here relies on ω 0 ω √ γ 1 [cf.Eqs.(H3) and (H5)], the γ → ∞ limit is to be understood as We would also like to emphasize that the scalings in Eq. (H7) are not limited to Lorentz-Drude and exponential cutoff functions: they hold whenever the conditions (H5) and (H6) are satisfied [although note that the scalings in Eq. (H7) will hold also when the condition (H6) is weakened to 0 < lim ω→∞ g(ω) lim ω→∞ g(ω) < ∞], which, we believe, is what happens generically.
[Discard old bath, attach new bath, start the next cycle from ρ p ⊗ R B ].

For
the disconnecting work, W d = − Tr H (CL) I (CL) ∞

FIG. 1 .
FIG. 1. Figures of merit of the device as functions of the coupling strength.(a) Efficiency, η (CL) , and (b) ergotropy, E (CL) , as functions of γ , for different values of the temperature T .The cutoff function is of Lorent-Drude form and ω 0 = 2, ω c = 4.

) with I 2 denoting the 2 ×
2 identity matrix.If we want s 1 • • • s d , then we add a ↑ superscript to σ and s: σ