Collective operations can extremely reduce work fluctuations

We consider work extraction from $N$ copies of a quantum system. When the same work-extraction process is implemented on each copy, the relative size of fluctuations is expected to decay as $1/\sqrt{N}$. Here, we consider protocols where the copies can be processed collectively, and show that in this case work fluctuations can disappear exponentially fast in $N$. As a consequence, a considerable proportion of the average extractable work $\mathcal{W}$ can be obtained almost deterministically by globally processing a few copies of the state. This is derived in the two canonical scenarios for work extraction: (i) in thermally isolated systems, where $\mathcal{W}$ corresponds to the energy difference between initial and passive states, known as the ergotropy, and (ii) in the presence of a thermal bath, where $\mathcal{W}$ is given by the free energy difference between initial and thermal states.

Fluctuations play a crucial role in microscopic systems, which has motivated enormous efforts to build a theory of thermodynamics for small fluctuating systems, both classical [1] and quantum [2][3][4]. Particularly relevant is the question of how to extend the celebrated second law of thermodynamics. A standard formulation of this law is that the average work consumed or extracted from a system in a process is bounded by its change of free energy. Since this law concerns average work, one can naturally ask what can be said about its fluctuations. This question has lead to very interesting insights, which notably include fluctuation theorems [4] and extensions of the free energy to the nanoscale regime [5,6]. This rich behaviour is smoothed out in the macroscopic limit, where fluctuations can usually be neglected with respect to average quantities. To formalise this limit, consider N identical copies of a quantum system, ρ ⊗N , and assume that the same work extraction process is implemented on each ρ. Then, the law of large numbers ensures that work fluctuations have a size proportional to √ N , whereas average work scales as N . The relative size of fluctuations will hence decay as 1/ √ N , regardless of the specific protocol being implemented. The universality of this behaviour builds a clear intuition that (work) fluctuations disappear in the limit 1/ √ N → 0. In this article, we also consider work extraction from ρ ⊗N , but away from the thermodynamic limit 1/ √ N → 0, so that N can be small. Crucially, we consider processes where the N copies of ρ can be processed collectively, hence assuming a high level of control. In this case, using concentration results from probability and information theory [7,8], we find that work fluctuations can disappear exponentially fast in N . This is first derived in thermally isolated systems, where the amount of extractable work from a state ρ is given by its ergotropy W erg (ρ) [9], i.e., the energy difference between ρ and its passive state [10,11]. While the concepts of ergotropy and passivity have been explored for average quantities [9-14], our results show that a large proportion of W erg (ρ ⊗N ) can be extracted without fluctuations by means of global operations. We note that the relevance of work fluctuations has been characterised in the creation of correlations from work [15], and in the charging of a quantum battery consist- * marti.perarnau@mpq.mpg.de ing of an harmonic oscillator [16]. We also study the decay of work fluctuations in the presence of a thermal bath. In this case, fluctuation-free protocols have been extensively studied in the last years within the field of single-shot thermodynamics [5, 6, 17-25], and our results provide new insights on the mesoscopic regime [23,24], where N is finite and possibly small. Our considerations concern states that are diagonal in the energy basis, as the definition of work fluctuations in coherent systems is subtle and an active area of research [26].
Work extraction from N two-level systems. We start by considering an illustrative scenario where work is extracted from a system S made up of N identical qubits (spins). This is modelled by putting S in contact with a a work repository W, which can extract or give energy to S. The total Hamiltonian , and H W = L j=−L w j |w j w j | with w j = jν and L ≥ N . Energy is transferred from S to W via unitary operations that satisfy which ensures that energy extracted/consumed from S comes solely from W . The initial state of SW is taken to be a product state, with ρ = (1 − p)|0 0| + p|1 1| being the single-qubit density matrix. We also assume p > 1/2 so that the qubits have initially population inversions and can transfer energy to W. Work is quantified as energy exchanges in the state of W, which starts with a well-defined energy [27]. Given (2), the probability to extract an amount of work w through U reads P (w) = w|Tr S U σU † |w , i.e., the probability that W has raised from 0 to w. Our goal is to maximise P (w) over all protocols, and hence we definẽ where the maximisation runs over U 's that satisfy (1). The computation of (3) is provided in Appendix A, and here we provide a sketch of it. We take w = kν and consider transitions in degenerate energy levels of SW, which is imposed by (1). Considering a global energy E SW = jν, with j = 0, ..., N , the populations of SW are distributed as: • Initial state: W is at |0 . Then there are C j N ≡ N j states of SW at energy E SW = jν, each with probability p(j) ≡ p j (1 − p) j .
• Target state: W is at |kν , and there are C j−k N states at energy E SW = jν. Note that initially these states are not populated.
The protocol that achievesP (kν) is the one that moves as much probability as possible from the initial to the target state in each degenerate energy E SW = jν. Since we only consider diagonal states, it is enough to consider permutations within each subspace. This leads tõ By noting that min we can use Hoeffding's inequality [8] to finally obtain (see Appendix A for details) To interpret this result, note that W qbit in (5) corresponds to the extractable work on average from ρ, i.e., its ergotropy [9] ( see also discussion in (9)). Hence, this result shows that we can extract deterministically a proportion 1 − γ of the total average extractable work from ρ ⊗N with a failure probability that decays exponentially with γ 2 N . The result (6) is a direct consequence of concentration results in probability theory, in particular of Hoeffding's inequality. Yet, it is important to distinguish it from the standard thermodynamic limit in which the relative size of fluctuations disappears as 1/ √ N . To illustrate this point, consider a simple protocol where work is extracted from each qubit individually. Energy is transferred from S to W by iteratively applying the energy-preserving unitary that swaps |0, j + 1 ↔ |1, j ∀j, where |a, b ≡ |a S ⊗ |b W , successively for each of the N copies. After N iterations, the state of W is described by σ W = N j=0 C j N p(j)|2j − N 2j − N |. The work distribution P (w) = w|σ W |w is hence a binomial distribution centred at w = N (2p − 1)ν with standard deviation µ = ν N p(1 − p). Hence, this protocol is essentially probabilistic and becomes only reliable by considering a confindence interval N W qbit ± O( . This is contrast with the fast convergence to a deterministic process of global protocols suggested by (6). This difference appears in protocols that extract less than N W qbit , see Fig. 1.
This simple example demonstrates the capabilities of global protocols for reducing work fluctuations. This is illustrated in Fig. 2, where we plotP (ω γ ), with ω γ = N W qbit (1 + γ) , as a function of γ and for different N 's. Note, for example, that one can extract 2/3 of the average extractable work from ρ ⊗N with success probabilityP ≈ {0.90, 0.92, 0.97, 0.99} by collectively processing N = {10, 25, 50, 100} copies of the state. On the other hand, the success probability of protocols that attempt to extract more than N W qbit rapidly decays to zero. This is made rigorous in Appendix A where we show thatP (ω +γ ) ≤ O exp −N γ 2 .
The goal of this article is to show that the behaviour explained for this example is in fact a generic property of workextraction protocols. In what follows, we show that with where γ ∈ (0, 1) and W is the maximal extractable work on average from a diagonal state ρ; whereas c and d in (7) are parameters that depend on the initial conditions (the initial state ρ and Hamiltonian H S ). Expressions (7) and (8) naturally connect average bounds with almost-deterministic protocols for work extraction. A direct consequence of this result is to provide a bound on the minimal number of copies that we need to globally process in order to extract ω −γ with success probability 1 − , N ≥ ln(d/ )/(γ 2 c 2 ). Disappearance of work fluctuations in thermally isolated systems. Consider that S is made up of N identical (dlevel) qudit systems, each of them with a Hamiltonian h = d i=1 ν i |i i| and an initial diagonal state ρ = d i=1 p i |i i|. The total state reads ρ ⊗N and the total Hamiltonian is non- The state is thermally isolated, so that work can only be extracted via controlled operations. Then, a crucial quantity is the ergotropy of ρ [9] where ρ pas is the passive state associated to ρ [10, 11], with p ↓ i being the eigenvalues of ρ ordered in decreasing order. W erg (ρ) quantifies the average extractable work from ρ by means of unitary operations as Tr(h(ρ − U ρU † ) ≤ W erg (ρ), ∀U [29]. That is, it quantifies how much accessible energy on average is stored in ρ.
Now we move to probabilistic work-extraction protocols. For that, as in the previous section, we consider an auxiliary work-repository W, which can accept energy from S. The possible work values are given by energy differences of H S , which can be written as with n ≡ (n 1 , n 2 , ..., n d ), n i ∈ [0, N ], i n i = N and similarly for m. The spectrum of W is taken to be non-degenerate and dense enough to accept all possible w n,m . We now proceed to find a lower bound onP (w) as defined in (3), where the maximisation runs over all unitaries satisfying (1). The possible values of w are given by (11) with all possible w n,m . We focus on a subset of those, defined by where γ ∈ (0, 1), and the ... ensures that the k i are natural numbers. This can be understood as a work-extraction protocol which moves the probability of the S' initial states with energy i j i ν i to target states of S with energy i (j i − k i )ν i , ∀j. As a consequence, W raises from 0 to w = . Given this protocol, the calculation ofP (w) is conceptually similar to the previous qubit example, but becomes more involved and is carried out in Appendix B. There, using arguments based on typicality [7, 8, 30], we show that (7) and (8) are satisfied with where c, at leading order in 1/N and γ, is given by where h i ≡ p ↓ i /p i and where S(ρ||σ) is the relative entropy. The expression (13) for c is correct up to corrections of order An exciting property of passive states is the phenomenon of activation [10-14, 31], which means that W erg (ρ ⊗N ) can be larger than N W erg (ρ). That is, by collectively processing ρ ⊗N more average work can be extracted. Interestingly, Alicki and Fannes showed in [31] that for N → ∞, where ω β (h) is a thermal state, ω β (h) ≡ e −βH /Tr(e −βH ), whose temperature is defined by S(ρ) = S(ω β (h)) with S(ρ) = −Tr(ρ ln ρ) the Von Neumann entropy [32]. Note that (14) can be also obtained through catalytic transformations [14], and that the state after the transformation is a completely passive state, i.e., a Gibbs state [10-12]. A natural question then arises: Can one extract a large proportion of the fundamental bound (14) almost deterministically? In Appendix C, we answer this question positively. More precisely, we consider a protocol defined by where γ ∈ (0, 1), and p th i are the populations of ω β (h), i.e., ω β (h) = i p th i |i i|. For such k th i 's, we then show that (7) and (8) are satisfied with W = W glob erg (ρ), d = d S , and where c can be inferred from (13) by replacing ρ pas ↔ ω β (h) and This value of c holds up to corrections O(1/N ), which are expected to appear as (14) can only be achieved for large N .
We illustrate these results in Fig. 3, where we compute ex-actlyP (w) for 50 copies of a qutrit system. From the figure it becomes clear that there are many possible protocols corresponding to different choices of w n,m , and that the choices (12) and (15) perform particularly well, especially the latter.
In fact, numerical results suggest that (15) rapidly becomes optimal with increasing N .
Disappearance of work fluctuations for systems in contact with a thermal bath. We now move to the case where S may be put in contact with a thermal bath B at temperature 1/β. In this case, the average extractable work from ρ with internal Hamiltonian h is given by the non-equilibrium free energy change [33][34][35][36][37], where and ω β (h) = e −βh /Tr(e −βh ), the temperature being now provided by B. The expression (16) bounds the extractable work on average of a generic ρ with B, and can be interpreted as a formulation of the second law of thermodynamics with a single bath. Note also that N W th (ρ) = W th (ρ ⊗N ) when the total Hamiltonian H S is a sum of non-interacting h's, so there is no possible activation here. Now we move to probabilistic work-extraction protocols on Our goal is to findP (w) in (3), and in particular to prove (7). In Appendix D, we boundP (w) in (3) by combining the techniques developed in Ref.
Our results hence show that by sacrificing a proportion (1− γ) of the extracted work allows for exponentially enhancing its purity. This is possible through collective operations. To give a feeling of these results, in Appendix E we compute exactlyP (w) for a qubit system, which in particular shows that one can extract 2/3 of the free energy of ρ ⊗N with a success probability {0.92, 0.96, 0.98, 1.0} by processing ρ ⊗N collectively with N = {10, 25, 50, 100} (see also Fig. 4 in Appendix E).
While so far we have assumed γ to be a finite fixed number, let us now consider the limit γ N →∞ −→ 0. In particular, by choosing γ 2 = ln(N )/(c 2 N ), we obtain simultaneously: Hence bothP (w) → 0 and w → N W th (ρ) in the thermodynamic limit 1/  (18), between the success probability and the extracted work. These different asymptotic scenarios are known as small, large, and moderate deviation regime in the context of information theory (see [41] and references therein).
Finally, let us briefly comment on the notion of work used here, based on a transition 0 → w in the state of W [5,18]. When transitions between more energy levels of W are considered, apparent violations of the second law can appear [20, 42,43], which can be avoided by either accounting for the entropy increase in the state of W [20,42] or by restricting to protocols that act on W in a translationally invariant manner [37]. In this sense, we note that (i) sincẽ P (w) → 1 the entropy gain can be neglected in most cases, and that (ii) it is possible to adapt the protocols derived here to satisfy translational invariance in W without modifyingP (w) [44]. Another relevant question is the presence of quantum coherence in ρ. While the coherent part of W th (ρ) can not be extracted for U 's that commute with the total Hamiltonian, it can by considering global protocols on ρ ⊗N [45][46][47] -note that, since such protocols do not rely on B, they can also be applied to W glob erg (ρ). This is yet another interesting effect of collective protocols.
Conclusions. We have shown that collective operations can extremely reduce work fluctuations, which can decrease exponential with N -the number of copies being processed. This is in contrast to the standard "thermodynamic limit", where the relative size of fluctuations decays slowly as 1/ √ N . The exponential decay (7) of work fluctuations has been proven for generic systems both in the presence and absence of a thermal bath, providing a bound on the exponential decay for each case. This result contributes to our understanding of global effects in quantum thermodynamics, and in this sense we note that previous results have shown that global operations can extract more average work [31, [48][49][50][51][52][53][54][55] and in a faster manner [56][57][58][59][60][61][62].
A high level of control has been assumed to obtain these results and arguably one of the most interesting questions is to which extent they can be observed experimentally. As recently proposed in Ref. [59], a promising candidate are architectures based upon the Dicke model [63], which allow for generating genuinely collective interactions. Another interesting proposal is the one of Ref. [ [27] (), at this point, the separation between work and heat may be a bit ambiguous, since energy acquired by W might have fluctuations, and hence its entropy can increase. However, in our results these fluctuations turn out to be exponentially small.
[28] (), other protocols that act locally on each copy of ρ ⊗N have a similar scaling w ∝ N and µ ∝ √ N due to the law of large numbers.
[29] (), it is interesting to note that the ergotropy also bounds the average extractable work in the presence of W, as given by the average energy change in W. For that, one needs an extra assumption on the unitaries allowed, namely they must commute with the translation operator on W [37]. In this case, the effective action of the operations on S corresponds to a mixture of unitary operations [65], from which no more work than the ergotropy can be extracted. We start by describing in detail how to reach (6) in the main text, where the goal is to extract w = kν. Recall that, at global energy E = jν (j = 0, ...,): • Initial state (ρ W = |0 0|). There are C j N , j = 0, ..., N states with population p j (1 − p) j .
• Target state (ρ W = |w w|). There are C j−k N , j = k, ..., N + k states with zero population. We wish to move as much probability from the initial state to the target state. The success probability can be written as We first focus on the case kν = W qbit (1 − γ), with W qbit ≡ (2p − 1)ν so that k = N (2p − 1)(1 − γ) with γ ∈ (0, 1). We will use Hoeffding's inequality [8], which in the particular case of a Binomial distribution with probability of success p and N runs, states that the probability p(k) of obtaining k ≤ N p successes is bounded by p(k) ≤ exp(−2(N p − k) 2 /N ). Here, this implies for kν = W qbit (1 − γ), as desired.
Let us now investigate the case where kν = W(1 + γ), so that k = N (2p − 1)(1 + γ). Note that here γ = (0, 1−p 2p−1 ). The second term in (A1) can be easily bounded as where we used again Hoeffding's inequality [8]. On the other hand, the first term in (A1) can be upper bounded as where we used well-known results on the tails of binomial distributions [66], and with By plotting the function f (p, γ) we find that f (p, γ) > 0 for p ∈ (0.5, 1) and γ = (0, 2 1−p 2p−1 ). Furthermore, we can expand f (p, α) for small γ obtaining Putting everything together, we obtain the desired resultP (W +γ ) ≤ O(e −γ 2 N ). We consider that S to be made up of N identical qudit systems, with a Hamiltonian The possible values of work are given by energy differences within S, where n i , m i ∈ [0, N ], are natural numbers with i n i = i m i = N . We assume that the spectra H W of the work repository W system is dense enough so that it can accept the values of work (B1). One may take it to be continuous, H W = dx x|x x|, but in principle this is not necessary. It is convenient to introduce the vector notation k ≡ {k 1 , k 2 , ..., k d }, and The average extractable work from ρ (under unitary operations on S, but note also [29]) is bounded by its ergotropy [9], where p ↓ i are the probabilities p i ordered in decreasing order, p ↓ i+1 ≤ p ↓ i , ∀i. We wish to extract deterministically an amount of work w = (1 − γ)N W by acting globally on ρ ⊗N and W. For that we take the n i , m i in (B1) to satisfy where the is introduced to ensure that the n i are natural numbers. In what follows we use that and assume (B7) for the sake of mathematical tractability (when we consider numerical simulations we use (B6)). It is important to recognize that there are many other choices of n that could lead to the desired w = (1 − γ)N W. For example, one can simply increase n 1 while taking n j = 0 with j > 1. The choice (B7) will however turn out to be crucial for the proof. We now consider the calculation ofP (w). Proceeding as in the main text, we have that at a global energy E = i n i ν i of SW: • Initial state (ρ W = |0 0|). There are C n N states with population p(k) for each k.
• Target state (ρ W = |w w|). There are C n−k N states, all of them with zero population.
Then, since we can only move energies in degenerate levels, we obtaiñ In what follows, we find min C n N , C n−k N within the typical subspace of S (see, e.g. [30]). To define the typical subspace, let |i 1 , ..., i N with i j ∈ {0, 1, ..., d} be an energy eigenstate of S, and n l = N j=1 δ ij ,l the number of l's in |i 1 , ..., i N . The (strong) δ-typical subspace T δ S is defined by eigenstates whose n typ i satisfy p i − n typ i /N ≤ δ ∀i (and assuming p i = 0). For convenience, let us also make the formal choice where c is still to be defined. Then, we can write a typical state as with i c i = 0 and c i ≤ c. Combining (B10) with (B7) we obtain In order to find min C n typ N , C n typ −k N , we consider where we applied Stirling's approximation N ! = √ 2πN N N /e N (1 + O(1/N )). Expanding for small γ we further obtain In order to find a bound that is independent of the α i 's, we can maximise the denominator of (B15). For that, we can define the vector h i = p ↓ i /p i , i = 1, ..., d and order it decreasingly Defining, we finally arrive at the choice for δ (B9) Here we multiplied the term O(1/N 2 ) by γ because for γ = 0 we have that C k typ N = C k typ −n N , and hence it must vanish. Note that the choice (B18) only depends on γ and on the initial state, defined by p ↓ .
Coming back to (B8), this result allows us to write, Hence, all that is left is to evaluate the population in the δ-typical subspace, i.e. the probability that |n i /N − p i | ≤ δ, ∀i. For that, we note thatP where in the first inequality we used union bound of probability theory, and in the second inequality we used the Hoeffding's inequality. Because δ is proportional to γ as in (B18), we obtain the desired decay in N γ 2 with corrections of order N O(γ 3 ) and 1 N O(γ 2 ). Absorbing the factor 2 into c, we obtained the desired result from the main text.
Appendix C: Passivity, complete passivity, and work fluctuations.
In this section, we extend our previous results when (B4) is substituted by the global bound where ω β (h) is a thermal state, ω β (h) ≡ e −βh /Tr(e −βh ), whose temperature is defined by S(ρ) = S(ω β (h)) with S(ρ) = −Tr(ρ ln ρ) the Von Neumann entropy. It is convenient to expand ω β (h), In order to extend our previous considerations, the key idea is to chose the n i in (B5) as and again we assume for the sake of analytical tractability. Proceeding as in the previous section, one then obtains, where we now keep only highest orders in N . Expanding over γ as before, where we used that S(ρ) = S(ω β (h)). At this point, the proof proceeds in complete analogy with the previous one, so one easily arrives at the final result Appendix D: Lower bounds onP in the presence of a bath Let us start by making explaining H B following [5]. Denote by e B an energy of B and by g(e B ) its degeneracy, and recall that B is assumed to be in a Gibbs state ω β (H B ). Then, we assume that there exists a subset E of the whole spectrum {e B }, in which g B (e B ) increases exponentially and satisfies g B (e B − e S ) = g B (e b )e −βe S , where e S is any energy of Ssee [5] for more details. For short-range interacting systems, the subset E contains 1 − of the population of ω β (H B ), with approaching 0 with the size of B [67]. Therefore, in what follows we assume that B is sufficiently large so that ≈ 0, and hence H B is well described by the subset E.
Having defined B, we now proceed as in the last section and consider work extraction from ρ ⊗N , with ρ = d i=1 p i |i i| and H S = d i=1 ν i |i i|. The initial state of SBW being ρ ⊗N ⊗ ω β (H B ) ⊗ |0 0|. FindingP (w) essentially boils down to moving states in the degenerate energy levels of SBW. Let us fix a total energy E for SBW and note that (using the notation (B2) and (B3)), • Initial state. Given W at |0 , e S = i n i ν i and e B = E − i n i ν i , there are g(E)e −β i niνi C k N states at total energy E, each with probability Z −1 e −βE e β i niνi p(n) and ∀n.
• Target state. Given W at |w , e S = i n i ν i and e B = E − i n i ν i − w, there are g(E)e −βω e −β i niνi C n N states at total energy E (each of them initially with 0 probability), and ∀n.
We now wish to move the probabilities of the initial to the target state, which has e −βw less states, in each degenerate energy level E. Since the scenario is identical for each global energy, we can focus on a particular E with a degeneracy g(E) = g, which we take to be a very large number. By normalising the probabilities within this subspace, we can reformulate the setting as: In a subspace of SBW with fixed global energy, the initial state of SBW has normalised populations of the form p ini. (n) = g −1 e β i niνi p(n) (D1) with degeneracies with C n N = N !/ i n i ! and n i = 0, .., N with i n i = N . On the other hand, in the target state, where W is at ω, has degeneracies g tar (n) = ge −β(ω+ i niνi) C n N .
The total number of states in the target state is then given by To findP (w), in principle one has to move the N tar biggest populations of the initial state to the target state. This leads to a set of inequalities that are directly related to the concept of thermomajorisation [5]. In order to obtain a simple and analytic result, here we will instead move the typical subspace of S, hence finding a lower bound onP (w). Recall that the (strong) δ-typical subspace T δ S is defined by eigenstates that satisfy n l ∈ N [p − δ, p + δ], ∀l with δ > 0 where n l = N j=1 δ ij ,l is the number of l's in the state |i 1 , ..., i N . The number of states in T δ S can be upper bounded by e N (S(ρ)+Cδ) with C = − i ln p i for p i ∈ (0, 1) The goal now is to find the largest δ that guarantees N tar ≥ N (δ) typ.ini. , so that we can move the whole population of the typical space to the target space. This easily leads to Putting everyting together, and proceeding exactly as in (B20), we obtaiñ with w = N W th (ρ)(1 − γ), and where we used that W th (ρ) = F (ρ) − F (ω β (h)). This shows the desired exponential decay with N γ 2 , which is multiplied by a function that depends on β and on S through H S and ρ S .