Random Quantum Batteries

Quantum nano-devices are fundamental systems in quantum thermodynamics that have been the subject of profound interest in recent years. Among these, quantum batteries play a very important role. In this paper we lay down a theory of random quantum batteries and provide a systematic way of computing the average work and work fluctuations in such devices by investigating their typical behavior. We show that the performance of random quantum batteries exhibits typicality and depends only on the spectral properties of the time evolving operator, the initial state and the measuring Hamiltonian. At given revival times a random quantum battery features a quantum advantage over classical random batteries. Our method is particularly apt to be used both for exactly solvable models like the Jaynes-Cummings model or in perturbation theory, e.g., systems subject to harmonic perturbations. We also study the setting of quantum adiabatic random batteries.

In a closed quantum system, a battery can be modeled by a time-dependent Hamiltonian H(t) evolving from an initial H 0 to a final H 1 . The system is initialized in a state ρ and, given that the entropy of the battery is constant under unitary evolution, the work extracted is given by the difference between the initial and final energies as measured in H 0 [2].
In this paper, we lay down the theory of Random Quantum Batteries (RQB). The randomness lies in the initial state ρ, the Hamitonian defining the units of the energy H 0 , and the time-evolution operator U t . We are concerned with the average work extractable by (or storable in) such a device and its fluctuations.
The main results of this paper are: (i) proving a typicality result for the extracted work in a large class of time dependent quantum systems. We show that -as the dimension n of the Hilbert space becomes large -the extracted work is almost always given by the difference in energy between the initial state and the completely mixed state, amplified by a quantum efficiency factor 1 + Q t /n 2 that depends solely on the distribution of the eigenvalues of the exponential of the time-dependent perturbation operator K. For Q t = 0, this result can be obtained by a classical system at infinite temperature. A random quantum battery can do it with limited energy resources. A non vanishing Q is a contribution that is purely quantum and depends on the constructive interference between different eigenvalues of K. The second main result is (ii) to provide a general method to study the average extractable work and its fluctuations in perturbation the-ory, which is essential to obtain results for physically relevant systems beside few exactly solvable models. We study as an example the Jaynes-Cummings model with a harmonic perturbation. Finally (iii), we study the case of adiabatic random quantum batteries, that is, batteries that operate slowly, so that there is no inversion of the populations of the energy levels. We show that also adiabatic random quantum batteries feature typicality in the large Hilbert space dimension n limit.
There is a large interest in typical properties in batteries due to the effect of disorder and the environment. In [34] a model of quantum battery based on a spin chain is studied where randomness is introduced as disorder in the couplings of the Hamiltonian H 0 . In [35] the disorder is introduced in the interaction Hamiltonian which is chosen to be in the Many-Body Localized (MBL) phase. In [36,37], the work statistics in the scenario of a random quantum quench are computed, and it is shown that the knowledge of the work statistics in this setting yields information on the Loschmidt echo dynamics. The importance of work fluctuations in quantum thermodynamics in a different setting than ours was also studied in [38].

II. SETUP
In this section, we are studying the typical behavior of random batteries when the interaction Hamiltonian is a random operator. The importance of this approach lies in the fact that typicality is a powerful argument to establish general features in quantum thermodynamics. As an example, typicality of entanglement in Hilbert space can be used to explain thermalization in a closed quantum system [18]. On the other hand, this approach is useful to argue about the robustness of a model of quantum battery.
We model the quantum battery in the following way. Start with a finite dimensional Hilbert space H = C n , and timedependent Hamiltonians H(t) ∈ B(H), that is, a bounded Hermitian operator on H. The initial state of the system will arXiv:1908.08064v6 [quant-ph] 18 Feb 2022 be denoted by ρ and its time evolution by ρ t = U t ρ ≡ U t ρU † t , where the unitary evolution operator is given by the timeordered product U t = T exp(−i t 0 H(s)ds). We model the Hamiltonian in two ways. In the first scenario we consider the time-dependence as a perturbation of a time-independent Hamiltonian H 0 , that is, H G (t) = H 0 + V G (t). The subscript G indicates the randomness of the perturbation which we take to be V G (t) = G † V (t)G, where G is a unitary representation of the unitary group on C n . In the second scenario we consider the time evolution generated by adiabatic evolution induced by a Hamiltonian H G (t), where the G t is a family of unitary operators that rotates the projectors onto the subspaces of a given energy. The discussion of the adiabatic scenario is deferred to section IV C.
In both settings, we can similarly model randomness in the initial state ρ or Hamiltonian H 0 also by random rotations ρ G = GρG † and H G = G † H 0 G. Loosely speaking, we will refer to the spectra of the initial state, of the measuring Hamiltonian H 0 , and of the evolution operator collectively as the battery spectrum. Notice that all these randomizations preserve the battery spectrum. This is a crucial point in this paper, as we are interested in ensembles of quantum batteries with a given spectrum. Randomizing also over the spectrum will yield, as we shall see, trivial results. In our setting the system is closed and evolves unitarily and the entropy of the battery does not change. Thus the work extracted from the quantum battery is given by (or ergotropy [1,39]). As mentioned before, this approach is different from the type of disorder in the couplings considered in the literature [34][35][36][37], as for us disorder is a random rotation G that mantains the spectrum of the eigenvalues of the perturbation V G (t) (the interaction).
A simple example which clarifies how our disorder affects the extracted work W is the following single spin case inspired by nuclear magnetic resonance (NMR). We consider a Hamiltonian of the form is the external magnetic field and σ the Pauli matrices. The spectrum of the interaction is effectively dependent only on the norm of the external field B, which can however be directed in all directions. We focus on an average which keeps the spectrum of the interaction constant, but rotates its basis. A two-level system example is the Jaynes-Cummings model of optics, on which we focus our attention in a random electromagnetic background. Precise experiments in these systems exist and thus provide a good background for testing the typical behavior of (random) quantum batteries [20,21].
The Hamiltonian H 0 defines the energy measurement, that is, the amount of energy stored in the battery. If we had access to any possible random Hamiltonian H(t), we would expect that the average state ρ t after the evolution should be the completely mixed state, in which case the average work extracted would be W = E 0 − Tr H 0 /n. This work is positive (that is, the battery has discharged) if the initial energy was larger than the energy in the completely mixed state, or it has charged if the initial state was populating the lower levels of H 0 . Notice that this setting we have arbitrary hamiltonians H(t) that can access arbitrary high energies as measured by H 0 . Instead, we ask how much work can be extracted if we have limited energetic resources, that is, when the spectra of H 0 and V (t) are fixed. This motivates our setting in terms of rotations of the time dependent part of the Hamiltonians as In the following, we are interested in the average work obtained by averaging over initial states ρ, the measurement of energy Hamiltonian H 0 , and the time dependent Hamiltonian H G (t). The averages are performed according to the Haar measure on C n . The fluctuations of work are defined through the same Haar averaging as ∆W 2 = (W − W ) 2 . In the following, the symbol X will represent the Haar average where G U is the suitable representation of the unitary group. We use standard techniques for the Haar averaging (see e.g. [40][41][42]) to compute the average and variances according to the Haar measure.

A. Work and quantumness
A quick calculation shows that These expressions imply that the extractable work depends on the lack of commutativity between the initial state ρ, the evolution operator U t , and the Hamiltonian H 0 . Moreover, they show that the coherence of the initial state in the eigenbasis of the evolution operator is necessary to have non vanishing extractable work from a quantum battery [43][44][45][46][47][48][49][50][51]. In particular, if the initial state is a steady state for the unitary evolution, the work is identically zero and so are work fluctuations. It is interesting that coherence in two different bases plays a role, which calls for a multi-basis definition of coherence from the resource theoretic point of view. In the following, we will see that this lack of commutativity takes the form of out of time order correlators, which is a hint to the connection between performance of quantum batteries and quantum chaos [36]. Notice that these expressions are also valid in the interaction an expression that will become useful later. Bounds on the stored and extracted energy have been obtained recently in [52], also in terms of the quantum Fisher information for the power P t = d dt W . As we remarked above, with no limit on energetic resources one can bring the system on average in the completely mixed state. A quantum channel that just dephases the system and mixes up the populations can achieve the same final result. The same result can be obtained by a classical system working at infinite temperature. Consequently, we are also interested in whether quantum coherence plays a specific role in outperforming the mixed state case. As we shall see, partial revivals due to the build-up of quantum coherence provide a quantum advantage.

III. AVERAGE WORK AND FLUCTUATIONS IN RQBS
In this section, we show how the average work and its fluctuations behave in quantum random batteries when we randomize over the initial states ρ, the measuring Hamiltonian H 0 , or the interaction V (t). In all cases, this average is obtained by rotating these operators by a random unitary operator and by taking the Haar average. Let us start by computing the average work obtained by a generic quantum evolution and averaging over all the initial states. It should not be surprising that the average extracted work amounts to zero. Indeed, we have where we defined the traceless operator δH 0 ≡ H 0 − U † H 0 U and have used that the Haar-average state in the Hilbert space is ρ = 1/n 1l. However, the fluctuations are not trivial [53]. Details of the calculation are in VI A. We obtain It is remarkable that the maximum of the fluctuations are reached for a pure state whereas they decrease with the purity of the initial state, and are identically zero if the system is initialized in the completely mixed state. Similarly, fluctuations in the work are smaller the larger the fluctuations in the eigenvalues of H 0 . Notice that the time dependent part has the form of a (two-point) out-of-time-ordered correlator (OTOC) [30,54]. What happens instead if we choose randomly the measuring Hamiltonian H 0 ? As we said above, we model this family of Hamiltonians as H G = G † H 0 G. This is a sensible definition as it gives us results that still depend on the spectrum of the Hamiltonian. Again, it should not surprise that the average work is zero, since as the average of every operator in the trivial representation is proportional to the identity, and ρ − ρ t is traceless. Some tedious calculations in Appendix VI B show that the work fluctuations are given by where ∆H 2 0 = 1 n Tr H 2 0 − 1 n 2 (Tr H 0 ) 2 are the fluctuations of the eigenvalues of H 0 , namely the fluctuations of H 0 in the completely mixed state. Again, the time-dependent part Tr (ρρ t ) has the form of an OTOC. The connection between OTO correlators and Loschmidt echo has recently been investigated in [32]. In terms of the 2−norm fidelity F 2 (ρ, σ) = Tr (ρσ)/ max[Tr ρ 2 , Tr σ 2 ] and the Loschmidt echo L t = F 2 (ρρ t ), we have Notice that as L t is typically scaling as n −2 [55], the average fluctuations are determined only by the fluctuations in H 0 and the purity of the initial state. However, at specific, revival times, there is a spike in fluctuations. Moreover, if we consider the average work over a large time T , the average Loschmidt echo becomes the purity of the the completely dephased state in the basis of the Hamiltonian,ρ and the above expression reads where the time average over a time T is defined as f T ≡ We see that large fluctuations can be achieved if there are not only large fluctuations in the energy gaps of the Hamiltonian H 0 , but also if the initial state is pure enough, or if the time evolution is nontrivial. If the initial state is very mixed or the time evolution does not feature an exponentially decaying Loschmidt echo, then work fluctuations will be negligible regardless of H 0 .
At this point, we are ready to tackle our main goal, that is, to compute the work and its fluctuations in a quantum battery modeled by H G (t) = H 0 + V G (t). In this setup, one has perfect control on the measuring Hamiltonian, but the controlled quantum evolution is very noisy, as V G (t) = G † V (t)G. However, one has retained control on the spectrum of the driving Hamiltonian, which is an experimentally realistic situation. In the interaction picture, and by defining C ≡ Tr [U I ρU † I H 0 ], we see that work is given by We can write the above expression as Now recall that the interaction picture operator U I depends on the random rotations G as GKG † . The average work W (t) V over the noise G can then be computed (see Appendix VI A for details) to give with where λ k = exp(iθ k ) are the eigenvalues of the evolution operator K = T exp(−i 4 of the work is thus contained in the function Q t . For large dimension n, the average work reads At this point, averaging over the initial state ρ would give zero, while averaging over the Hamiltonian H 0 gives an exponentially small work ∼ n −1 . Let us comment on the meaning of the result Eq. (12). We are starting with an initial state ρ and evolving with a random evolution generated by V (t). So far we have averaged over rotations of the time dependent perturbation V G (t). Such rotations keep the eigenvalues of V G unchanged so that all the results are a function of spectral quantities like Q t . One could expect that, if the evolution were completely random, one would end up with the completely mixed state, and then the work extracted would have to be W = (E 0 − Tr H 0 /n). However, we have fixed the spectrum of V (t) in the randomization, so it is remarkable that one can achieve the infinite temperature result.
Moreover, in the average work W (t) V there is an amplifying quantum correction (1 + Q t /n 2 ). These corrections are quantum in nature because they correspond to the constructive interference that builds up in Q t = −2 j =k cos(θ j − θ k ). One expects that without a specific structure in the θ's, the factor Q t /n 2 would rapidly decay to zero. This means that on average (and typically) one can achieve in this setting the same result that would be attained with random arbitrary resources. However, we can do better than that. First, if fluctuations are not a concern, it is possible for nano-systems with small n to have large Q t . We are going to give an example in the following, using an optical cavity. Moreover, it is possible to design devices with a spectrum such that, for specific values of t, the term Q t is of order one, which can be exploited as quantum advantage in the construction of a battery. In the next section, we show how, in a specific example, revivals in Q t allow the battery to outperform the infinite temperature (and classical) behavior.
The question of what happens in the large n case is very interesting. In the optical cavity application shown in the next section IV.A, the quantum amplifying factor is washed out as n −2 . We think that this would happen for most models. In this sense, this is a sign of the loss of quantumness as the dimension of the Hilbert space grows. One wonders, though, whether for some specific model the amplifying factor Q t /n 2 might not disappear in the large n limit. Finding such a realistic model would be of enormous practical interest. Conversely, proving that no model can feature this advantage as n goes to infinity would be a very interesting result in quantum thermodynamics.
As mentioned, one expects that for a random matrix its spectrum should yield a vanishing Q t . A natural question to ask then is what the typical behavior of this quantity is when these eigenvalues are taken randomly, according to a CUE distribution, (see e.g. [36]). Let us define r k = λ k+1 /λ k . We prove in Appendix VII E that The behavior of Q, evaluated numerically, is depicted in Fig.  1. We see that for large n the peak of the distribution moves towards zero. That is, averaging over the spectra does not give any amplification Q t . How typical is the behavior of a random quantum battery in the large n limit? If there is typicality, an optimal strategy for random quantum batteries would consist in fixing the optimal spectrum of K and then knowing that the other details of the evolution will not matter in the large n limit. To this end, we need to compute the fluctuations which is far more challenging because they involve the fourth tensor power of the unitary representation. We find that and a lengthy calculation yields where the inequality follows from the fact that (T (2) ) ⊗2 ≤ 1. Putting together all the terms, we find in Appendix VII B that the fluctuations are upper bounded by where M (n) is an upper bound to the terms of the form | mnop e i(θm+θp−θn−θo) |. If one chooses spectral properties for K such that M (n) = O(1), then the fluctuations scale like n −2 and thus a many-body quantum battery would show exponentially small fluctuations. Moreover, this is the typical case. Indeed, by averaging over CUE to compute M (n), we see in Fig. 1 that this quantity is concentrated near zero for large n. More in depth numerical evidence is provided in Appendix VII B, where we analyze numerically every single term which contributes to the fluctuations, showing that indeed every single term converges to zero for large n's. This represents the first main result of this paper: random quantum batteries show typicality in allowing a work extraction given by the difference in energy between initial state and completely mixed state, amplified (or attenued) by the form factor 1 + Q t /n 2 . By thus choosing a suitable V 0 , one can obtain with probability almost one the desired behavior for work extraction in the sense of the Haar measure on GV 0 G † .
The specific behavior of Q determines whether the quantum advantage in a random battery is washed out or not in the large n limit. We now apply these findings in the case of an exactly solvable model and study the behaviour of Q. We consider a two-level system in an optical trap described by the Jaynes-Cummings model [56]. In the rotating wave approximation only two adjacent modes at time (n, n + 1) of the electromagnetic field couple with the two level system (details provided Appendix VII C). For this calculation, we assume that the atom couples with a finite set of modes of electromagnetic field, which we truncate at a number n = 2R, where R is a truncation of the number of modes of the electric field. At the end of the calculation we will send R → ∞.
The Hamiltonian reads where we define ∆(t) = Ω(t) − ω(t), and we assume g(t) = g 0 e M t . For this model, we find the eigenvalues exp(iθ k ) exactly and use them to evaluate Eq. (13). Following the calcula- We then obtain the average work Eq. (12) where, as seen above, the function Q(α t ) is a sum of trigonometric functions whose complete expression is given in Appendix VII C, Eq.(106). We note that when Q > 0, effectively the system extracts more work than the classical counterpark. In this sense, Fig. 2 (a) shows that there can be a quantum advantage in a specific model. In Fig. 2 we plot the time evolution of the extracted work from the random Jaynes-Cummings battery averaged over V . As we can see, for most times the quantum efficiency gets washed out. For small n, at specific revival times given by inverting Eq. (106), the value of Q becomes of order one, and thus providing a non-vanishing quantum efficiency. This is at the price of performing much worse at different times. One can design a quantum battery by an array of many random nano-batteries of small n and evolve to the revival time where the work extracted goes above that corresponding to the maximally mixed state [3,57,58]. The fact that non-vanishing Q is obtained as revivals in Eq.(13) is a sign that this amplification comes from the constructive interference coming from the complex eigenvalues of K and therefore of its quantum nature. On the other hand, for large n, the system almost always behaves like in the limit of the battery that completely mixes the state, though one has obtained this performance with limited, realistic resources that do not require to bring the system at infinite temperature.

B. Time dependent perturbation theory
In the case of the Jaynes-Cummings model we could solve for the time evolution exactly, finding expressions for the average work and its fluctuations via perturbation theory. We make use of the Dyson series for the evolution operator in the interaction picture, namely U We consider perturbations up to the second order in the Dyson series, and at this point we can average over G. Define the operator A = t t0 V 0 (t )dt . Again we need the fluctuations of A in the completely mixed state, namely n 2 ∆A 2 = n Tr A 2 − (Tr A) 2 .
Averaging over G requires a lengthy calculation (see Appendix VII D) yielding The second term is the difference between the initial energy and the energy in the completely mixed state.
As an example consider the case of an exactly solvable Hamiltonian H 0 subject to the Harmonic perturbation where we have defined f (t, ω) = 2sin( t−t0 2 ω)/ω and λ's are the eigenvalues ofV . As one can see, the average work decreases with n. We plot W (t) V in Fig.3. In this model it is easy to find the revival times at which the quantum efficiency is maintained also for larger values of n. One can indeed show (see Appendix VII D) that the work performed by a random harmonic perturbation of the form 2V cos(ωt) has always a single maximum at t k = (2k + 1) π ω on average.

C. Adiabatic Quantum Batteries
Now let us consider the case of a quantum battery performing an adiabatic evolution connecting the two Hamiltonians H 0 and H 1 and the two respective equilibrium states ρ 0 , ρ 1 , e.g., two eigenstates or Gibbs states for H 0 , H 1 (but also thermal or more general mixed equilibrium states). Adiabatic evolution as a method to perform quantum computation [59] or quantum control has been long an important tool in quantum information processing, see, e.g., [60]. Adiabatic evolution to perform work extraction was studied in [61]. A model for an adiabatic quantum battery based on a three-level system was studied in [63]. In this section, we deal with general adiabatic quantum batteries in which the adiabatic drive is rotated in a random direction as a function of time.
In general, two Hamiltonians are adiabatically connectible if and only if they belong to the same connected component of the set of iso-degenerate Hamiltonians [62]. By denoting H α = R i=1 i α Π i α (α = 0, 1) the spectral resolution of H 0 and H 1 , and ordering their eigenvalues in ascending order i.e., 1 α < ... < R α . We define the vectors where the continuous unitary family {U t } 1 t=0 is such that U 0 = 1 1 and U 1 = U . The work extracted after the adiabatic evolution thus reads because the populations in the i−th subspace are conserved by the adiabatic evolution. We now have Π i α Π j β = δ ij if α = β, but otherwise they are not necessarily orthogonal. We see that the work depends on the choice of U as We can now perform the average over the unitary transforma-tion U . We easily obtain To understand the role of the degeneracies, let us consider the case of a non degenerate Hamiltonian, so that d i = 1 for all i. We obtain W ad = E 0 − Tr H 0 /n, which again is the difference between the initial energy and the energy of the completely mixed state and thus the quantum efficiency is washed out (see [64]). More generally, as we show in Appendix VII F, we find an upper bound on the adiabatic work given by so that potentially random adiabatic quantum batteries could give an advantage over classical devices as well (even at infinite temperature), as c ≥ 0.
Let us now look at the fluctuations ∆W 2 ad . The calculation involves averaging the square of the work and thus the order two tensored representation of the Unitary group. This is also a lengthy calculation, whose details are given in VII F. We obtain For n 1, the terms of order 1/n 3 go to zero faster than 1/n 2 , and we obtain which shows that random adiabatic quantum batteries feature typicality. Fluctuations during adiabatic driving were studied in a different context also in [65].

V. CONCLUSIONS AND OUTLOOK
In this paper we provided a notion of quantum random batteries by means of Haar averaging initial states, Energy measurement Hamiltonian, and the time-dependent driving of the quantum battery. This method allows to study large classes of systems, including not-exactly solvable systems or adiabatic quantum batteries. The average work and fluctuations are systematically studied; we find that quantum batteries exhibit typical behavior in the large n limit given the spectral properties of the driving system. On average, the work extracted is found to be typically equal to the difference between the energy of the initial state and that of the completely mixed state, amplified by a quantum efficiency factor 1 + Q t /n 2 that only depends on the spectrum of the driving Hamiltonian. Quantum efficiency is not washed out at specific revival times for small systems. Our method allows for the computation of Q t in perturbation theory, therefore allowing for the treatment of realistic systems. We have also treated the case of random adiabatic quantum batteries, finding that amplification is lost for a non-degenerate Hamiltonian.
In perspective, our results put forward several questions that we would like to investigate in the immediate future. We have shown that for small systems there are revival times in which quantum coherence builds up and gives a quantum advantage. Typically, this is not the case for large n. However, it is an open problem whether there are random quantum batteries whose spectral properties allow for the build-up of coherence that outperforms the classical case. Conversely, showing the impossibility of such quantum amplification for large n would be an important result in quantum thermodynamics. This is a problem which we plan to explore in the near future in a realistic model. A second question relating to the effect of quantum coherence also arises. As we have seen, the extracted work can be related to the coherence of the initial state in two different bases, or of the operator U t in two different bases. This suggests that there is a non trivial interplay between coherence and work that involves more than one basis [66]. Also, the lack of commutativity between the initial state and the evolution operator or the measuring Hamiltonian and the evolution operator take the form of out of time order correlators. It would then be interesting to explore the connection between fast decays of these quantities, chaos, scrambling, and work statistics. One very intriguing insight comes from the fact that the narrowing of fluctuations does shrink the quantum efficiency but at specific revival times. These revival times correspond to spectral properties of the time evolution operator and one would be interested in understanding the connection between quantum efficiency of random quantum batteries and the integrability or chaotic behavior of the Hamiltonian. Using tools from local Haar averaging [41], we can explore whether the efficiency in a battery with a microscopic local drive is influenced by quantum chaos or integrability. The optimization of the path in a adiabatic quantum algorithm is related to the brachistochrone or geodesics in the space of the ground state manifold [67]. It would be very interesting to see if optimal paths correspond to bounds given by quantum thermodynamics. Finally, it would be important to generalize these results to the case of open quantum systems.
Acknowledgments Also, see the notes on spectral methods available at https://www.overleaf.com/read/wdxknyfmdnww.

B. Work fluctuations averaging on H0
Let us define R = ρ − ρ t . We consider the fluctuations on the work via the averaging on the operator H 0 . We have where now the coefficients of the projectors are λ ± = Tr (Π ± H ⊗2 0 )/Tr Π ± . Direct calculation gives, defining a ≡ (Tr H 0 ) 2 and b ≡ Tr The work fluctuations can thus be written as Now, consider the fluctuations ∆H 2 0 of the eigenvalues of the Hamiltonian H 0 , namely the fluctuations of H 0 in the completely mixed state 1l/n. We have we then obtain ∆W 2 H0 = 2 n(n 2 − 1) n 2 ∆H 2 0 Tr R 2 and thus finally which is the result we report in the paper.

C. Traces of K
A direct calculation of the coefficients yields Moreover, we use that We now see that, defining and thus Using the relationships 1 2n we get D. Calculation of W (t) V and ∆W 2

V
The work extracted W (t) reads We can write the above expression as The average work over the noise G can then be computed as The unitary operator K = T exp(−i t 0 V (s)ds) will be diagonalized in the form K = k exp(iθ k )|k k|. Using the usual technique, we find G ⊗2 (K ⊗ K † )G †⊗2 = ± λ ± Π ± , where now λ ± = Tr (Π ± K ⊗ K † )/Tr Π ± . Notice that in this setup, already the average work involves the average over the tensored representation of the unitary group. We obtain We finally obtain where In the above equation, exp(iθ k ) are the eigenvalues of the evolution operator K = T exp(−i t 0 V (s)ds). All the time dependence of the is thus contained in the function Q(θ j − θ k ).
The fluctuations are more challenging because they involve the fourth tensor power of the unitary representation. Let us set out to find them. We see that where C ≡ Tr [U I ρU † I H 0 ]. The relevant object to compute is then This time, the average reads with λ i = (Tr Π) −1 Tr (Π i K ⊗2 ⊗ K †⊗2 ). Now, the Π i are the projectors onto the irreps of S 4 . There are five irreducible irreps of S 4 . In the next subsection we show an explicit expression of these projectors. A lengthy calculation yields

E. Irreps of S4
Let us first recall the character table of S 4 in Table I. The last row of Table I gives the size of each conjugacy class in S 4 . Given a permutation σ ∈ S 4 , we denote by S(σ) the representation of S : S 4 → GL(H ⊗4 ) given by By the Schur-Weyl duality the projectors onto its irreps are where χ (r) is the character of the (r) irrep of S 4 and χ (e) is the dimension of the irrep in S 4 . The five projectors are given by: In the above, the symbol + . . . denotes a sum over all the members of the conjugacy class. As well known, the five conjugacy classes of S 4 are given by their cycle structure of Table II

A. Main definitions and projectors
First, let us note that for any operator A, we have |Tr(AT (2),⊗2 )| ≤ |Tr(A)| since T (2),⊗2 ≤ 1, and since we care about the scaling of the fluctuations we can focus on a simpler calculation which does not involve the swap. Before we begin the calculation, we start with a few definitions which will be useful in the following: We then start with the construction of the projectors in a basis, which we take as the computational basis: Since we are interested only in the scaling with n of the fluctuations, we focus on the structure of the traces and not on the proportionality constants. Using the definitions above, the projectors can then be written explicitly in the computational basis.
At this point, we can start the evaluation of the traces. First we note that (ρ ⊗ H 0 ) 2 = abcdef ρ ab c ρ de f |acdf bcef |. We then have: ρ aa c ρ dd f + ρ ad a ρ de f + ρ ac d ρ da f + ρ aa f ρ dc d + ρ aa d ρ df c + ρ af c ρ da d +ρ ad c ρ df a + ρ af a ρ dd c + ρ ac f ρ dd a + ρ ac a ρ df d + ρ ad f ρ da c + ρ af d ρ dc a + ρ ad a ρ df c + ρ af d ρ da c + ρ ac f ρ da d + ρ ad f ρ da a + ρ ac d ρ df a + ρ ac a ρ dd f +ρ af a ρ dc d + ρ ad e ρ da f + ρ af c ρ dd a + ρ aa d ρ dc f + ρ aa f ρ dd c + ρ aa c ρ df d +ρ ad c ρ df a + ρ af a ρ dd c + ρ ac f ρ dd a + ρ ac a ρ df d + ρ ad f ρ da c + ρ af d ρ dc a − ρ af a ρ dc d + ρ ad a ρ df c + ρ af d ρ da c + ρ ac f ρ da d + ρ ad f ρ da a + ρ ac d ρ df a +ρ ac a ρ dd f + ρ ad e ρ da f + ρ af c ρ dd a + ρ aa d ρ dc f + ρ aa f ρ dd c + ρ aa c ρ df d Tr[Π (st⊗sgn) (ρ ⊗ H 0 ) ⊗2 ] ∝ − ρ ac a ρ dd f + ρ ad c ρ da f + ρ af c ρ dd a +ρ aa d ρ dc f + ρ aa f ρ dd c + ρ aa c ρ df d + ρ af a ρ dc d + ρ ad a ρ dc d + ρ af d ρ da c + ρ ac f ρ da d + ρ ad f ρ da a + ρ ac d ρ df a − ρ ac a ρ df d + ρ ad f ρ da c + ρ af d ρ dc a (64) We can now evaluate the trace over the operator K ⊗ K † with the projectors, T r[ΠK ⊗2 ⊗ K †⊗2 ]'s. We have the following results: At this point we can calculate the average fluctuations, which can be written as · ρ aa c ρ dd f + ρ ad a ρ de f + ρ ac d ρ da f + ρ aa f ρ dc d + ρ aa d ρ df c + ρ af c ρ da d + ρ ad c ρ df a +ρ af a ρ dd c + ρ ac f ρ dd a + ρ ac a ρ df d + ρ ad f ρ da c + ρ af d ρ dc a + ρ af a ρ dc d + ρ ad a ρ df c +ρ af d ρ da c + ρ ac f ρ da d + ρ ad f ρ da a + ρ ac d ρ df a + ρ ac a ρ dd f + ρ ad e ρ da f + ρ af c ρ dd a +ρ aa d ρ dc f + ρ aa f ρ dd c + ρ aa c ρ df d Let us now assume that A is a projector with k non-zero eigenvalues. Then the inequality implies that where b n · · · b n−k are the highest k's eigenvalues values of B. We thus need to focus on the singular values of (ρ ⊗ H) ⊗2 . The eigenvalues of ρ ⊗ H, are e ij = p i j , and the eigenvalues (ρ ⊗ H) ⊗2 are e ijkl = p i j p k l . Since p i ≤ 1 in the most general case, e ijkl is upper-bounded by 2 max . We thus have that a conservative upper bound is given by where k is the dimension of the non-zero subspace of the projector operator. Since each term of the trace is divided by the dimension of the projector operator and we have four non-zero terms, we have in the most general case. However, the bound p i ≤ 1 is very loose. If the p i ≤ γ n , we have and thus there is concentration. For instance, we have concentration if we have that ρ is a mixed sate. A stronger bound can be done by using the expressions we derived. We see from the bound above that this is not enough to prove concentration. However, the concentration can be proven if the take advantage of the structure of the fluctuations in terms of the density matrix. We now provide an alternative proof of the same statement. From the previous subsection, we see that the work fluctuations F can be upper bounded as where C is a O(1) constant counting the number of all the terms in F , M (n) is an upper bound over al the terms of the type mnop e i(θo+θp−θm−θn) and kn 2 is the upper bound to the terms containing the ρ, H 0 : which is true because the projector operators are positive. Putting things together, we obtain

C. Jaynes-Cummings model
As seen in eqn. (12), the average work depends only on the value of the eigenvalues of the Unitary evolution operator K. Let us consider the case of an optical cavity interacting with a 2-state system. The optical cavity with the two state system (an atom) span(|g , |e ) can be described within the rotating-wave approximation using the Jaynes-Cumming Hamiltonian: It is immediate to see that [H, a + a + σ z ] = 0. Specifically, we focus on the interaction picture, in which H I = RHR † , where (in the rotating frame) we have R = e −iωt(a † a+ σz 2 ) , and one has a Hamiltonian described by H I = RHR † , with We define ∆ = Ω − ω. The operators a and a † act on the electromagnetic field, while σ's act on the two-level system. We have We now consider a wave function of the form where we will send R → ∞ at the end of the calculation. The time evolution of this system is given by the Schroedinger equation (in the interaction picture), which is of the form: Which is not hard to see that it can be written as whose solution is given by We note that V (t )V (t) = V (t)V (t ) in the case of a time dependent interaction Hamiltonian. In fact, we see that on the n−th subspace of the wave function, given the definition W (∆, g) = ∆g − g∆ of the wronskian of the functions ∆ and g, we have from which we observe that we can have a time dependent and commuting (at all times) Hamiltonian if we have the condition which can be satisfied if for a constant M . In this case, the time ordering can be removed and we can write The Stone operator in this case can also be written explicitly on each subspace. It can be shown that in each r-th subspaces and where Thus, the Stone operator which describes the time evolution on the rth subspace is given by We now focus on the eigenvalues of the matrix above, which must be of the form e iθ k . For a matrix of the type the eigenvalues are known exactly and are of the form λ ± = a ± i √ c 2 + d 2 . It is immediate to see that the eigenvalues are complex, and have norm 1. The phases are given by ±θ k ≡ ±β k (t). We thus find that which is what we need for the evaluation for the work in the main text. We can now plug this result into eqn. (12), which reads . We thus need to calculate j =k cos(α(j − k)). Thankfully, this sum is known, and is given bỹ from which we obtain: with We thus see that the time dependence of the work enters only inQ (α(t)).
It is not hard to see that revivals occur for z k = k with k ∈ N, thus for α k a multiple of 2π. We now have that for k ∈ N, from which we get the revival times as a function of M and g 0 .

D. Time dependent perturbation theory
In the case of the Jaynes-Cummings model we could solve for the time evolution exactly. This is rarely the case and we must resort to perturbation theory in most cases. Consider to start the definition of thw work: where we consider a Dyson expansion. In this case, the solution is given by the Dyson time ordering We are interested in the case in which we need to resort to perturbation theory to evaluate the unitary operator above. Up to the second order, we have In what follows, we can assume that V † 0 = V 0 . Given the expressions above, we have now to evaluate the average of using the average of the unitary matrix G: with Note that Tr(U 2 U † 2 ) = Tr(U † 2 U 2 ) = n + O(t 3 ). We can use at this point the eqns. (41) again. After a rapid calculation we see that (up to corrections of order t 3 ), we have and thus where A = t t0 V 0 (t )dt , where we used the fact that inside the traces one has Tr ( t t0 t t0 (: V 0 (t )V 0 (t ) :) † dt dt ) = Tr ( t t0 t t0 : V 0 (t )V 0 (t ) : dt dt ). We can now write As it could be seen from the beginning, we see again explicitly that the average work is the product of two terms, the first is adimensional and due to the perturbation, and the second term has the dimensions of energy, and due to the density matrix only: This shows that no work can extracted if the density matrix is the one of a completely mixed state.

Example: Harmonic perturbations
Let us now consider the example of a n-level system. At time t = 0, the system is described by the eigenvalue equation and thus the wavefunction as a function of time can be written as We consider now a harmonic perturbation of the form: whereV is a generic operator andV † its hermitean conjugate. Then, according to the formulae we have derived, the average work if we consider random rotations with respect to G of We now use: and thus, if we define f (t, ω) = 2 And thus the A dependent part of the average work is given by which is the expression for the performed work due to a harmonic perturbation. What we see is that the overall work is proportional to product of two functions, one is the square of function f (t, ω) = 2 sin( t−t 0 2 ω) ω and a factor which depends on the eigenvalues of the operatorV . The function f is periodic with period 2π ω and has a maximum for t k = (4k + 1) π ω + t 0 . IfV is self-adjoint, σ k = λ 2 k , and we have in the parenthesis the function k,k λ k λ k cos 2 ( t + t 0 2 ω) + 2n which can be rewritten as If we introduce the constants a 0 , b 0 , c 0 , the work is thus a function of the form : which is periodic. For t t 0 , the function above has two minima if c 0 < d 0 and only one for c 0 > d 0 . However, it is not hard to see that c 0 > d 0 is always true if is always true ∀A. However the identity above follows immediately from the fact that Tr (aA + bI) 2 ≥ 0 (134) is true for arbitrary a, b ∈ R, and it follows from the choice a = n, b = c ± Tr(A) with c ± = −1 ± 1 − n 2n + 1 .
Thus, the work performed by a (random) harmonic perturbation of the form 2V cos(ωt) has always a single maximum at t k = (2k + 1) π ω on average. This can be interpreted as the fact that there are specific moments at which we stop our process to have performed the maximum amount of work on the battery.
the average of Q, evaluated numerically, is provided in Fig. 1. We see that for large values of n the peak of the distribution moves towards zero.

F. Adiabatic Quantum Batteries
Here we give the details for the calculation of work fluctuations ∆W 2 ad for the adiabatic batteries. We first recall the calculation of the average. Let us start from the following protocol. The Hamiltonian, for α = 0, 1, is written for an adiabatic transformation as Consider i (t) : [0, 1] → R, with i (0) = i 0 , i (1) = i 1 . It can be shown that the evolution of the projector operators can be written as Thus, the time evolution of the Hamiltonian for an adiabatic system can be written as where the while the density matrix as ρ(t) = i p i U t Π i 0 U † t . It is important that the vector d α i ≡ (Tr(Π j α ) does not change with time, and thus can simply call d i these quantities, meanwhile n is the dimension of the Hilbert space.
Because these relationships are in a way independent from the intermediate states, we simply write these expressions for t = 0 and t = 1 without loss of generality. The work as We now have Π i α Π j β = δ ij if α = β, but otherwise they are not necessarily orthogonal. Let us write the work as We can now perform the average over the unitary transformation U . We obtain Since we will need it for the calculation of the fluctuations, we note that Let us now calculate the fluctuations. The square of the work reads