Trade-offs between precision and fluctuations in charging finite-dimensional quantum batteries

Within quantum thermodynamics, many tasks are modelled by processes that require work sources represented by out-of-equilibrium quantum systems, often dubbed quantum batteries, in which work can be deposited or from which work can be extracted. Here we consider quantum batteries modelled as finite-dimensional quantum systems initially in thermal equilibrium that are charged via cyclic Hamiltonian processes. We present optimal or near-optimal protocols for $N$ identical two-level systems and individual $d$-level systems with equally spaced energy gaps in terms of the charging precision and work fluctuations during the charging process. We analyze the trade-off between these figures of merit as well as the performance of local and global operations.


INTRODUCTION
The second quantum revolution [1] has brought about unprecedented access to technologies at the nanoscale, which are currently operating in what has been dubbed the noisy intermediate-scale quantum (NISQ) regime (cf.[2]).These advances go hand in hand with the desire to further improve the control over quantum systems and to better understand their potential and limitations for storing and processing information.At the same time, residual heat and noise are ever present adversaries in this endeavour, and moving systems away from thermal equilibrium with their surroundings requires sufficient control as well as the investment of time and energy.A framework that aims to address fundamental questions regarding the dynamics, interactions, energetics, and control of quantum systems in the presence of heat baths presents itself in the form of quantum thermodynamics [3,4].Indeed, from a fundamental thermodynamic perspective, pure states can only be prepared approximately since Nernst's unattainability principle [5][6][7] -the third law of thermodynamics-requires infinite resources to cool any system to its ground state.To accurately assess the resources required for a specific task, it must therefore in principle be assumed that the respective system is initially in a thermal state and that work must be invested to change this.
Quantum thermodynamics offers a broad spectrum of different scenarios to model such state transformations and corresponding work inputs, but two distinct paradigms can be identified as the conceptual polar opposites of this spectrum (cf.[7][8][9]): Work can be supplied to a target system (i) via a heat flow generated by a temperature gradient between two thermal baths, or (ii) via a direct supply from a coherent work source.The former scenario can be understood as the operation of a heat engine [10][11][12][13][14][15], where the heat flow supplies work incoherently to a working substance and the dynamics are globally energy-conserving.On the one hand, this paradigm is appealing from a thermodynamic point of view, since the system is overall closed and external control can be minimal in the sense that an external agent operating the machine is only needed to switch on (and off) interactions between the target and the heat baths.On the other hand, only a restricted class of state transformations is achievable within this paradigm (cf.[7]) and practical laboratory situations in which quantum technologies are employed are not typically operated using heat engines.
An all-encompassing understanding of possible state transformations and their resource costs must therefore include coherent work sources as in (ii).Although the specific realizations of these work sources are often not included explicitly in modelling state transformations, doing exactly this will ultimately be necessary to truly obtain fine-grained descriptions that will lead to a better understanding of quantum systems beyond thermal equilibrium.Such descriptions can be envisioned to provide insights, e.g., regarding the effects of finitetime transformations, finite-size reservoirs, and fluctuations of relevant quantifiers.A starting point for such a more general approach lies in modelling the work sources -commonly dubbed quantum batteries [16,17] -on their own, i.e., independently of the systems that they eventually supply work to.In other words, quantum batteries are considered as quantum systems in which work can be temporarily deposited and from which it can subsequently be extracted.
Here, we follow the approach of Ref. [46], and consider battery charging realized via cyclic Hamiltonian processes.In this case, the system Hamiltonian returns to its original form at the end of each cycle and the battery-system state, FIG. 1. Charging processes.Illustration of a unitary charging process for a battery in initial state τ and final state ϱ = U τ U † , showing the distribution of probability weights in the energy eigenbasis.The charging precision V (ϱ) only depends on the final state ϱ, whereas the fluctuations ∆W depend on the particular evolution from τ to ϱ.
initially assumed to be thermal as we have reasoned above, can be modelled to lie within the unitary orbit of the initial state [16].This has the advantage that it allows us to consider the charging process independently of the specifics of other potentially involved auxiliary systems (e.g., the charger systems as in [24,27], or external classical power sources).We can thus focus on the properties of the charging process and of the charged battery, and study fundamental bounds on the chosen figures of merit.
We further centre our attention on two particular quantities: the charging precision, quantified by the variance of the final battery charge, and the work fluctuations arising during the charging process.While the former concerns a property of the final state of the battery, independently of how this state was reached, the latter characterizes the particular charging process, as illustrated in Fig. 1.Nevertheless, at fixed final battery charge, both quantities cannot be simultaneously optimized by the same charging procedure [46] except for certain special cases 1 .In other words, for nonzero-temperature initial states, optimal precision generally implies non-optimal fluctuations, and vice versa.Consequently, it is of interest to derive optimal protocols for both precision and fluctuations, and to determine trade-offs between these figures of merit.
For quantum batteries realized by quantum harmonic oscillators, such optimal protocols have been derived in Ref. [46], but translating them to finite-dimensional quantum systems frequently considered in pertinent literature (cf.[19, 20, 23-30, 35-37, 42, 47-50]) has proven to be a formidable task [51].Here, we present advances towards closing this gap: we construct a general protocol that optimizes the charging precision as well as a protocol aiming to minimize the work fluctuations for quantum batteries consisting of N identical two-level systems or of individual d-level systems with equally spaced energy levels.We compare the performance of these protocols in terms of both figures of merit, charging precision and work fluctuations, to investigate potential tradeoffs between them.Our results represent a first step towards more all-encompassing future analyzes of the performance of work-storage strategies for quantum-thermodynamic systems that take into account more complicated energy-level structures as well as other relevant properties of the charging process (power and fluctuations) and of the battery itself (charging precision and stability of the charge).
This manuscript is structured as follows: in Sec. 2, we provide technical definitions for the systems under study and the relevant figures of merit.In Sec. 3, we then present the protocol achieving the fundamental precision limit for the considered systems, before turning to the problem of determining a similar protocol for fluctuations in Sec. 4. There, we provide a protocol that is motivated by insights previously gained for harmonic oscillators, evaluate it numerically and argue that it is a good approximation to the optimum.We then explore the trade-off between the two quantities in Sec. 5. Finally, in Sec.6, we study the role of local versus non-local operations in charging N -qubit quantum batteries by comparing the optimal global protocols to the worst-case local protocols, before we present our conclusions in Sec. 7.

FRAMEWORK
In this section, we first discuss the basic setup for the type of charging processes we consider and provide definitions for the relevant figures of merit, i.e., charging precision and work fluctuations, in Sec.2.A, before establishing some preliminaries and notation for the particular battery systems we will study in Sec.2.B, i.e., batteries consisting of N noninteracting identical two-level systems, here simply referred to as N -qubit batteries, and what we call qudit batteries, where the system to be charged consists of a single d-dimensional quantum system with equally spaced energy levels.

2.A. Precision and work fluctuations of charging processes
We consider a d-dimensional quantum system with associated Hamiltonian H = ∑ d−1 n=0 E n |n⟩⟨n| as a quantum battery, where E n and | n ⟩ represent the nth energy eigenvalue and its corresponding energy eigenstate, respectively.Without loss of generality, we assume these eigenvalues to be labelled such that energies are non-decreasing with increasing n, i.e., the E n are ordered such that E n ′ ≥ E n for n ′ > n ∀n, n ′ , and we set E 0 = 0. We further assume that the battery is initially uncharged, i.e., contains no work that is extractable via cyclic Hamiltonian processes.In such a process, the evolution of the system is given by a time-dependent Hamiltonian H(t) = H + V (t), where the cyclic potential V (t) includes all time-dependency of the Hamiltonian and satisfies the condition V (t f ) = V (0) = 0, in which t f represents the duration of the protocol.Such processes can be represented on the system Hilbert space as unitary operations [16].In general, systems from which no work can be extracted by such operations are in so-called passive states [52], that is, states whose average energy cannot be lowered by unitary operations.Any passive state must hence be diagonal with respect to the energy eigenbasis and its eigenvalues must be ordered non-increasingly with increasing energy.The notion of passivity can even be relaxed to restrict the operations used for work extraction, e.g., to Gaussian operations in continuousvariable systems [53].Here, however, we wish to further restrict the class of initial battery states to those that are completely passive: passive states for which any number of identical copies also remains passive.For a given Hamiltonian, the only completely passive states are thermal states, and we hence consider initial battery states in Gibbs form, where β = 1/T is the inverse ambient temperature, Z = Tr[e −βH ] is the partition function of the canonical ensemble, and we use units where ̵ h = k B = 1 throughout the manuscript.Besides being motivated by Nernst's principle as we have argued in Sec. 1, this choice of initial state thus ensures that no work is extractable before the charging process.
We are then interested in charging procedures realized by cyclic Hamiltonian processes and thus need to consider the unitary orbit of the initial thermal state, i.e., ( Since the initial state is passive, all unitary operations increase (or leave invariant) the average energy.We therefore consider increasing the average energy of the battery by a fixed amount For any given charge ∆E, there exists a set of unitaries transforming a fixed initial state to final states with the same average energy.This non-uniqueness of the final state and of the unitaries leading to it provides an opportunity for optimization of different quantities of interest, such as the charging precision, energy fluctuations, and charging speed, all of which may play important roles during charging processes.
The first quantity that we will focus on here is the charging precision.We describe this quantity by the energy variance of the final state, given by such that a smaller variance corresponds to higher precision and vice versa.We note that for any given unitary U , the final-state variance V (ϱ) is non-negative but should be viewed relative to the variance of the initial state τ (β), and the former may be smaller or larger than the latter.The unitarily achievable values of V (ϱ) for any fixed energy input depend on the temperature of the initial state [46].Moreover, for infinite-dimensional systems, the increase in variance is not bounded from above for any nonzero energy increase ∆E.
Here, however, we wish to analyze the fundamental precision limit for fixed values of ∆E for finite-dimensional systems and all quantities of interest are hence bounded.In addition, it is worth mentioning that the precision only depends on the initial and the final states, and the type of dynamics relating these states do not play any role in the characterization of the precision.This means that the optimization of the precision reduces to finding the state with the minimum variance for a given energy input.
The second quantity that we analyze here can be viewed as a figure of merit for the quality of the charging process: specifically, we are interested in minimizing the work fluctuations for a given energy increase.In general, there are a number of inequivalent ways to define work in the quantum domain 2 Here, we focus on work quantified by two-point measurement (TPM) scheme [61].There, two ideal3 projective measurements, one prior to and one after the action of the transformation (here represented by U ), are used to estimate the work performed on the system in terms of differences between the respective measurement outcomes E m and E n .The work value assigned to such a pair of outcomes is W m→n = E n − E m and the work estimate is then obtained as the average Here, p m→n is the transition probability of the energy eigenstate | m ⟩ to | n ⟩ starting from the initial state τ , and p m = ⟨ m | τ | m ⟩ is the probability to find the initial state in the eigenstate | m ⟩.For unitary processes such as those we consider here one finds that the average work matches the change in average energy, ⟨ W ⟩ = ∆E.The work fluctuations (∆W ) 2 of the charging process can then be quantified by the average squared distance from the average energy increase, As such, (∆W ) 2 represents the second statistical moment of the work distribution obtained in the TPM scenario.Thus, if one raises the TPM scheme to a definition of work, the work fluctuations provide a basis for confidence statements about the estimated (average) work input required to charge the battery, whereas the energy variance allows one to quantify one's confidence regarding the estimated charge itself.However, the relation between the work fluctuations and energy variances is generally complicated: in our case, the quantities satisfy where Based on Eq. ( 7), one can easily see that if the battery is initially in the ground state 4 , i.e., T = 0, the work fluctuations coincide with the final-state variance, i.e., (∆W ) 2 = (∆σ) 2 = V (ϱ), but in general, these quantities do not coincide and cannot be simultaneously minimized.Detailed comparisons of optimal protocols for precision and fluctuations have so far only been available for harmonic-oscillator batteries.Here, we therefore want to analyze and compare achievable performances for finite-dimensional batteries, whose specific compositions we will describe next.

2.B. N -qudit model: single-qudit and N -qubit batteries
As a specific realization of a quantum battery, let us now describe a multipartite system consisting of N identical noninteracting d-dimensional subsystems with equidistant energy levels (here referred to as 'qudits').We will then consider two special cases of such a system more closely: N -qubit quantum batteries (arbitrary N but d = 2) and single-qudit batteries (N = 1 but arbitrary d).
The N -qudit system is described by a set of local Hamiltonians Since the individual battery systems are not interacting with each other, the joint initial thermal state is an uncorrelated product state of the form τ tot (β) = τ d (β) ⊗N , where Taking into account the degeneracy of the energy levels of the joint system, the total Hamiltonian can be written as where m and i m indicate the mth distinct energy level and ith level with energy m, respectively, while g d (m) indicates the degeneracy of the mth distinct energy eigenvalue, such that E m,im = E m,jm for all i m , j m = 1, . . ., g d (m) and E m,im ≠ E n,jn for all i m , j n as long as m ≠ n.For instance, for an Nqubit system, the degeneracy of energies is given by g ).In terms of this notation for the eigenstates of the joint Hamiltonian, one can rewrite the initial thermal state as where the initial probability distribution on the diagonal with respect to the energy eigenbasis is independent of the index i m , i.e., p m,im ∶= e −βmω /Z d (β) N where p m,im indicates the initial probability weight of the eigenstate | m, i m ⟩.Using Eqs. ( 9) and ( 10), the initial average energy of the total system takes the form where we have made use of the fact that Due to the symmetry g d (n) = g d (N (d−1)−n) with regards to the sizes of the degenerate subspaces, identifying the maximal average energy within the unitary orbit of the initial state simply corresponds to rearranging the probability weights in such a way that the populations of the eigenstates | n ⟩ and /ω and ϵ ∶= E(ϱ)/ω to describe the initial and final energies, respectively, we then have ϵ 0 ⩽ ϵ ⩽ N (d − 1) − ϵ 0 .Since ϵ = ϵ 0 + ∆ϵ, this implies that the charge ∆ϵ of the battery satisfies In the following, we investigate the precision and work fluctuations achievable with unitary charging processes for two different special cases of the system described by this model: single-qudit batteries (N = 1) and N -qubit batteries (d = 2).

FUNDAMENTAL CHARGING PRECISION LIMITS FOR ARBITRARY TEMPERATURES
We are now in the position to determine the optimal protocol minimizing the variance for a specific charge ∆ϵ.To briefly reiterate, the (single-qudit or N -qubit) system is initially prepared in the state τ (β) with energy ϵ 0 = ϵ(τ (β)).Then the energy of the system is unitarily increased to reach the target energy ϵ = ϵ 0 +∆ϵ.The goal is to choose the unitary operation such that the energy variance V (ϱ) of the final state is minimized for fixed ∆ϵ and fixed inverse initial temperature β.However, direct minimization of V (ϱ) is generally difficult even for fixed initial temperature and input energy, owing to the number of parameters required to describe the involved unitaries [64].This problems is exacerbated by the fact that the unitaries achieving the minimum variance are not unique and our desire to specify the result as a function of both ∆ϵ and β.We therefore employ an optimization protocol for V (ϱ) that makes use of an auxiliary quantity Ṽ , the average Step I of the optimal-precision protocol.The four cases (i)-(iv) of reordering the initial probability weights (a) are illustrated in panels (b)-(e), respectively, by varying the target energy ϵ (and hence k) for fixed N = 3.The vertical axes show the sizes of the respective probability weights, labelled by s = 0, 1, . . ., 7, while the horizontal axes show the energy levels (degeneracies are indicated by groups of vertical lines) as well as the initial energy ϵ0 in (a) and the different target energies ϵ in (b)-(e) as dashed lines.
squared distance (ASD) from the target energy, which we define as for a given probability distribution {p n,in } n .In general, the quantities Ṽ (ϵ) and V (ϱ) do not coincide, but when the distribution {p n,in } n,in matches the probability distribution of the final state ϱ with respect to the energy eigenbasis, i.e., when In this way, the optimization can be carried out in terms of a protocol that aims to minimize the ASD with respect to a fixed target energy in every step, while the average energy changes throughout the protocol and only reaches the target value at the end.The ASD thus allows us to obtain the optimal protocol for the variance in a simple way.This protocol can be divided into two distinct steps: I. In the first step, illustrated in Fig. 2, the initial probability distribution is rearranged such that the larger probability weights are assigned to energy levels closer to the target energy ϵ.That is, the resulting distribution {p (I) n,in } n,in satisfies p (I) m,jm ≥ p (I) n,in for all m, n with |m − ϵ| < |n − ϵ|, and hence corresponds to the minimum ASD with respect to ϵ in the unitary orbit of the initial state.II.In the second step, unitary two-level rotations that minimally increase the ASD per unit of energy shift are used to adjust the average energy to match the target energy ϵ.
In the following, we first describe these steps in more detail for the N -qubit system.It is then straightforward to adapt the N -qubit protocol to a single-qudit system of dimension d by considering the former protocol for N = d − 1 with the additional replacement g 2 (m) ↦ 1 for all m.For ease of notation we will drop the subscript on the degeneracy factor for qubits from now on and use g(m) instead of g 2 (m).
For the sake of notational simplicity, we further define the new variable to label the eigenstates of the joint system using only a single index s = 1, 2, . . ., 2 N .Since each value of s uniquely identifies a pair of values {m(s), i m (s)}, we use the notation p(s) ∶= p m(s),im(s) such that p(s) ⩾ p(s ′ ) for all s ⩽ s ′ with s, s ′ ∈ {1, 2, ⋯, 2 N }.This allows us to order the probability weights with respect to the energy eigenstates in non-increasing order using only the parameter s.
Step I of the protocol: In the first step, we rearrange the initial-state probability weights p(s) to form a new probability distribution {p(s)} s , such that the largest value, p(1), is associated with the energy level closest to ϵ, the second-largest weight, p(2), is associated with the second-closest level and so on, to reach the minimal value of Ṽ (ϵ) in the unitary orbit of the initial state.In order to do so, we first need to find the closest energy level (labelled k) to the desired target energy ϵ, which is given by where ⌊ϵ⌋ and ⌈ϵ⌉ denote the floor and ceiling functions, i.e., the closest integers to ϵ that are smaller or larger than ϵ, respectively.
We can then identify four different cases, labelled (i)-(iv) in the following, depending on the signs of the quantities ϵ − k and ⌊ N 2 ⌋ − k, where the latter represents a constraint arising from the finite system dimension.That is, the details of how the probability weights are reordered depend upon whether k is closer to the lowest or to the highest energy level.For all four cases, the resulting density matrix after step I is diagonal with respect to the energy eigenbasis, and the corresponding diagonal probability weights are given by: Case (ii): Case (iii): Case (iv): From the resulting probability distribution we obtain the average energy which is generally not equal to the desired target energy, εI ≠ ϵ, but may be smaller or larger than ϵ in any of the four cases (i)-(iv).Consequently, the ASD from the desired energy is generally different from the energy variance for the distribution arising from step I. Therefore, the energy of the system must be changed to reach the target ϵ, which is the purpose of step II.
Step II of the protocol: Now, we want to adjust the average energy by using a sequence of unitary two-level rotations.Each of the transformations slightly alters the average energy to achieve ϵ, but since the ASD was globally minimized (within the unitary orbit of the initial state) by the first step of the protocol, step 2 will increase the ASD.We hence select the transformations such that each of them increases the ASD only minimally per unit of energy change.Here we need to find the optimal sequence of these two-level rotations.
To do so, we first consider a two-level rotation between two arbitrary levels m and n with weights pm and pn and energies E m and E n , respectively, and parameterize the rotation by an angle θ.Starting from a diagonal density matrix (with respect to the energy eigenbasis), this transformation can be represented as the map The associated energy change is given by Similarly, the change in the ASD is From this expression we see that we have to apply a twolevel rotation between levels n and m, chosen such that ( Em+En ω −2ϵ) is minimized while bringing the average energy closer to ϵ, in order to obtain the minimum ASD increase per unit energy.To identify these pairs of levels, it is convenient to choose a relabelling relative to the value k from Eq. (15).That is, instead of m and n, we introduce the variables l ∈ N 0 and j ∈ Z such that m = k − l and n = k + l + j.The average energy change associated to the two-level rotation in Eq. ( 21) can then be written as Using Eq. ( 23), we can also obtain the change of Ṽ per unit energy change, i.e., We thus see that the variable j determines a hierarchy of possible two-level rotations that increase Ṽ the least, while ensuring that ∆ε II and ∆ Ṽ ∆εII have the desired sign.That is, when εI < ϵ, the average energy needs to be increased, suggesting that we have to select index pairs such that pm,im > pn,in and 2l + j > 0 while minimizing j.In contrast, when εI > ϵ, we should select levels with pm,im < pn,in and 2l + j > 0, making the maximal values of j desirable.According to these rules, we can identify the optimal value j opt of j and the corresponding set of admissible values l opt,i of l for any given probability distribution p(s).Since the energy levels can be degenerate, for any fixed pair (j opt , l opt,i ), corresponding to some variables (m, n), one can then find two sets of labels, im=1 with the desired properties.In contrast to the non-degenerate case discussed in Appendix A.1.II of Ref. [46] this means that there are now x max,i ∶= min{g(k − l opt,i ), g(k + l opt,i + j opt )} possible pairs of levels between which one may rotate for any fixed choice of (j opt , l opt,i ), and we label these pairs by a subscript x, i.e., (j opt , l opt,i ) x with x = 1, . . ., x max,i .For a given probability distribution p(s) we can thus generate a set P opt given by For each pair of levels in P opt , we can then perform a twolevel rotation.In principle, one has the freedom to adjust the angles of all possible rotations specified by the pairs in P opt individually to approach the desired target energy.For instance, one may perform individual operations one after the other with maximal angles θ = π 2 until one is close enough to the target energy so that the adjustment of a single rotation angle to a suitable value 0 < θ < π 2 reaches the value ϵ exactly.Irrespective of the order or particular distributions of these angles, the resulting energy variance is always the same, as long as the target average energy is reached.However, we note that different choices of these angles may lead to different results as far as other figures of merit for the charging process are concerned.In particular, adjusting the angles individually can result in discontinuities of the work fluctuations associated with the protocol as a function of ϵ at the transition points between the cases (i)-(iv) above.
Here, we therefore choose a common rotation angle θ for all pairs of levels in P opt .If the target energy can be reached by a suitable choice of θ, then one selects this value, performs the operations, updates the probability distribution and the protocol is concluded.If the chosen rotation does not reach the target energy, one carries out the rotations with θ = π 2 for all pairs in P opt , updates the probability distribution and generates the corresponding new set P ′ opt .This procedure is continued until the target energy is reached.
Not only can we obtain the optimal-precision protocol in this way, we also observe that the final probability distributions change continuously at the transition points at the values ϵ = (n + 1) 1  2 for n ∈ N. As described, this approach for obtaining the optimal-precision protocol is independent of the degeneracy of the energy levels.Therefore, it can be applied for any system whose Hamiltonian can be written in the form of Eq. ( 9), for instance, N identical qudit systems with equally spaced Hamiltonians.
In panels (a) and (b) of Fig. 3, the results of the protocol are illustrated for a 5-dimensional system (d = 5, N = 1) and a 4-qubit system (d = 2, N = 4, same ω), respectively, showing the minimal unitarily achievable energy variance V (in units of ω 2 ) as a function of the energy input ∆ϵ.In both case, the systems are initially in thermal states with respect to their local Hamiltonians with temperatures (in units ̵ hω/k B ) from 0.1 to 1 in steps of 0.1.From Ref. [46], it can easily be seen that the fundamental variance limit for pure initial states (zero temperature, ground state) is From this formula we conclude that if the energy input is an integer multiple of ω, the variance vanishes and its maximum value is an integer multiple of ω 2 .In both plots, we see that for temperatures close to zero (in our case for T = 0.1), the minimum variance is well approximated by Eq. (27).When the initial temperature is raised, the exact periodic behaviour of the minimal variance for the harmonic oscillator discussed in Ref. [46] disappears.Instead, one now observes that the finite system dimension leads to a symmetric behaviour with respect to the point where the input energy ∆ϵ is exactly half-way between its minimal and maximal values as specified in (12).The positions of the local minima and local maxima of the minimal variance V (∆ϵ) for each fixed Hamiltonian depend on the initial temperature.
Here, it is interesting to observe two competing effects when comparing panels (a) and (b): For zero initial temperatures, the d-dimensional system and a (d−1)-qubit system appear as equivalently useful work-storage devices as far as the maximally achievable charge ∆ϵ max and the minimal achievable variance V (∆ϵ) are concerned.However, for higher temperatures, certain trade-offs become evident.On the one hand, the degeneracy in the energy-level structure of the multi-qubit system means that ∆ϵ max is decreased more strongly with increasing T than for a single d-level system.On the other hand, the degeneracy also means that the global minimum, min ∆ϵ V (∆ϵ), remains at smaller values as compared with the d-level system.

FUNDAMENTAL WORK FLUCTUATION LIMIT FOR ARBITRARY TEMPERATURE
We now turn our attention to the work fluctuations during the charging process of a general finite-dimensional system.For ease of presentation in this section we write the system Hamiltonian as H = ∑ d n=1 E n |n⟩⟨n|, where E n are sorted in increasing order.The work fluctuations ∆W are obtained from Eq. ( 6), and we again focus on two particular cases (cf.Sec.2): a qudit with equally-spaced Hamiltonian, and a system consisting of N identical qubits, each with the same local Hamiltonian.
Qudit protocol.For this scenario, we first make an ansatz for increasing the system's energy by formulating a protocol that consists of a sequence of permutations of populations of pairs of energy levels, as illustrated in Fig. 4. The sequence is specified by a free parameter, m.The parameter is a non-negative integer and the energy differences that can be achieved by all sequences labelled by m hence form a discrete set (for fixed initial temperature).As a consequence, this ansatz cannot reach arbitrary final energies, and so further finetuning is required in subsequent modifications of the protocol.
Qudit protocol-Phase One.To increase the system's energy starting from an initial thermal state we start with a sequence of permutations that move the populations of the m energy levels with the largest energies to the m levels with the lowest energy.The starting point for this sequence is to exchange the smallest population p d of any energy level (which initially corresponds to the largest eigenvalue E d of the Hamiltonian) with the population of the adjacent energy level with lower energy E d−1 .Taking the new population of this, now second largest energy level, as a new starting point, we exchange it with the population of its adjacent lower-energy neighbour E d−2 and repeat this procedure, step-by-step, until we reach the lowest energy level E 1 .In the process, all other populations are shifted upwards by one energy level.
Before we proceed with the remaining m − 1 levels, let us motivate this procedure by considering the limiting case where E d → ∞, i.e., a harmonic oscillator (since we assume equally spaced energy levels).In this limit, the result of the first sequence of pairwise permutations is a shift of the average energy by exactly one unit.At the same time, the original population p d of the only energy level experiencing a shift different from one unit vanishes in this limit, lim E d →∞ p d = 0.As a result the work fluctuations associated to this limiting-case process vanish.For a finite-dimensional system with finite energy gaps we have p d > 0 and so the overall shift in average energy will be less than one, and the fluctuations will not vanish but will be proportional to the smallest population p d .For a harmonic oscillator, the procedure can be repeated any number of times to raise the average energy by any non-negative integer number of units with vanishing work fluctuations [46].
For finite dimensions and energy gaps we can also repeat the procedure carried out above for p d , now for the secondsmallest population p d−1 but stop when it has reached the second-smallest energy level E 2 , and similarly for all of the remaining m − 2 populations among the smallest m populations.As a result, the smallest m populations end up as the populations, in increasing order, of the smallest m energy levels.It is clear that the resulting density matrix at the end of this step is diagonal with respect to the energy eigenbasis.The new probability weights with respect to this basis are given by Here, m remains as a free parameter that allows us to adjust the average energy of the final state.At this point, the total energy ε(m) of the system reads where m ∈ {0, 1, ⋯, d − 1}.Note that if m = 0, the first sum in Eq. ( 29) does not contribute to the average energy.
To better understand the proposed protocol, let us now again make a comparison with the optimal protocol for a harmonic oscillator from [46] by considering the simple case of a qudit system with equally spaced energy levels.In this case, by considering E n = (n − 1) ω, one can rewrite the corresponding energy increase of the system for given m in Eq. ( 29) as In general, we know that ∆ε(m) = ε(m) − ϵ 0 is a monotonically increasing function of m that takes discrete values.Therefore we are not able to cover all possible energy increases ∆ϵ using this approach.To remedy this, we first need to find a parameter m for a given ∆ϵ which minimizes ∆ε I (m) ∶= ∆ϵ−∆ε(m) under the constraint that ∆ε I (m) ≥ 0, such that That is, we find the parameter for the protocol described above for which the energy that is reached is as close as possible but still smaller than the desired target energy.
Qudit protocol-Phase Two.During the second phase of the protocol we then compensate for the missing energy ∆ε I ( m) by transforming the probability weights associated with the levels n = m + 1, m + 2, ..., d according to These pairwise rotations (of the largest and smallest, secondlargest and second smallest, etc.) of the probabilities in the considered subset about a common angle θ result in an energy change with respect to the first phase of the protocol given by Due to the minimization in Eq. ( 31) the rotation between these levels must be sufficient to reach the desired energy.Using Eqs. ( 31) and ( 33), the required angle for the rotation to compensate for the rest of the energy is obtained from Performance of the protocol.We now wish to examine how well, in particular, how close to the optimum, the proposed protocol performs.To do this, we once again consult the case of the harmonic oscillator.Here this corresponds to the limit d → ∞ of a d-dimensional system with equally spaced energy levels.For such a situation, the protocol described above reduces to the protocol from [46] which was shown to minimize the work fluctuations for fixed energy input: There, m is fully determined by ⌊∆ϵ⌋ which leads to a minimization of all terms except for the last sum in Eq. ( 35).One therefore only needs to investigate the contribution from last term in the work fluctuation on its own.The crucial step of the protocol from [46] is then to realize that the probability weights in the mentioned sum have all been shifted from some level with label n to another level with label m < n, whose energy gap E m − E n might diverge.In particular, (E m − E n − ω ∆ϵ) 2 can diverge, but the associated probabilities go to zero much faster (exponentially with E m − E n ), lim n→∞ p n = 0.So the last sum containing the probability weights p n for n = d − m + 1, . . ., d in Eq. ( 35) vanishes in the case of the harmonic oscillator.
In the finite-dimensional case we consider here, however, all energy gaps and all weights p n remain finite and give a nonzero contribution to the last sum in the work fluctuation in Eq. (35).In this case, it is generally complicated to confirm the optimality of the constructed protocol.However, in the regime of small temperatures (with respect to the maximum energy gap) it can be argued that the last sum in Eq. ( 35) is negligible for sufficiently small input energy and the protocol thus (at least) approximates the true optimal protocol.Furthermore, it is obvious that the optimal fluctuation protocol for pure states of any qudit or N -qubit system consist of partial shifts of the ground-state population to levels ⌊∆ϵ⌋ and ⌈∆ϵ⌉ such that the corresponding energy is equal to ∆ϵ, which is compatible with the optimal precision protocol.
To check our approach quantitatively we have numerically calculated the work fluctuations arising from the proposed protocol for a qudit system with varying input energy and for different initial temperatures, and compared them with a brute-force numerical search for the corresponding optimal ) obtained from the protocol proposed in Sec. 4 as a function of the input energy ∆ϵ (in units of ω) for (a) a 5-dimensional quantum system, and (b) for a 4-qubit system, for temperatures (in units of ̵ hω/kB) from T = 0.2 (blue) to T = 1 (red) in steps 0.2.For the chosen parameters (input energies, temperatures, and system dimension/qubit number) results from a numerical optimization of the work fluctuations are included but cannot be distinguished by the naked eye from the curves resulting from our proposed protocol, supporting the conclusion that the latter is at least a close approximation of the true optimum.
protocol, as illustrated in Fig. 5 (a) for d = 5.For the parameters we have considered the numerical differences between our protocol and the optimum are non-zero but invisible to the naked eye.We therefore now continue with an analysis of the properties of our proposed near-optimal protocol.
In the regime where the input energy is small compared to the inverse temperature of the initial state or compared with the system dimension, one can observe almost periodic behaviour of the fluctuations as functions of the input energy, approximating the periodic behaviour of the harmonic-oscillator case [46].In particular, if the system is infinite-dimensional or in a pure state, one may reach arbitrary input energies while keeping the fluctuations bounded from above by ∆W ≤ ω 2 , and for all input energies that are integer multiples of ω one can reach ∆W = 0.
However, in finite-dimensional systems with finite temperatures, the periodic behaviour and local minima and maxima gradually disappears with increasing input energy.The fluctuation becomes a monotonically increasing function of ∆ϵ as one approaches the maximum of the energy that can be transferred to the system unitarily.
Protocol for N qubits.So far, we have discussed fluctuations for a single-qudit system with equally spaced energy levels.But, as illustrated in Fig. 5 (b), one can also observe the same qualitative features when applying the proposed protocol to a system of N non-interacting qubits.To this end, we relabel the eigenvalues and eigenstates of the N -qubit Hamiltonian via a variable that we have already encountered in the description presented at the end of Sec.2.B, i.e.
. With this, we can easily apply the protocol derived for qudit systems above.
In Fig. 5 we showcase the performance of this protocol for a 5-dimensional system and for a 4-qubit system in terms of the work fluctuations as functions of the input energy for different temperatures.These plots illustrate that for small input energies one observes the discussed approximately periodic behaviour (as one encounters for the harmonic oscillator [46]) in finite-dimensional systems.It is clear that the periodic behavior is a result of transferring probability weights close to zero from high-energy levels to low-energy level with a negligible fluctuation cost.It tell us that the number of periodic cycles is given by the number of energy levels with negligible probability weights.However, in contrast to the harmonic-oscillator case, the work fluctuations observed in finite-dimensional systems strongly increase for larger energy inputs and any residual periodicity is gradually lost.

COMPARISON OF THE PROTOCOLS
After introducing the protocols that minimize the variance and work fluctuations in Secs. 3 and 4, respectively, we now wish to investigate the trade-offs between the two quantities by checking how well the protocols designed for minimizing one of them perform in terms of the respective other quantity.For pure initial states, it is known [46] that the protocols coincide, and so both the variance and the fluctuations can be minimized simultaneously.For finite temperatures, the protocols generally do not coincide.In Fig. 6 we therefore show the energy variance and the work fluctuations associated to the optimal-precision and the fluctuation protocol as functions of the invested energy ∆ϵ for both systems of interest and for different initial temperatures.However, we note that no further optimization of either protocol is carried out here in order to further adapt it to the respective second figure of merit.
As expected, the protocols coincide for small temperatures (e.g., as can be seen from the blue curves on the very bottom of all panels in Fig. 6), but the differences between the protocols become apparent for increasing temperatures.For the ex- ample shown in Fig. 6, the differences between the variances obtained from the fluctuation protocol and from the optimalprecision protocol [panels (a) and (b)] are smaller on average than the differences in the work fluctuations from using the optimal-precision protocol as opposed to the fluctuations protocol [panels (c) and (d)] when taking into account the different units on the vertical axes of panels (a) and (b) with respect to (c) and (d).
For example, for the highest shown temperature (T = 1, red curves), the ratio From the examples we have considered it thus appears that it is better on average to employ the optimal-precision protocol if one wishes to keep both the variance and fluctuations low with equal priority.At the same time, we observe that the differences between the protocols become less pronounced when the system under consideration offers more degeneracy in its energy levels: The ratios for the highest-temperature curves for the 4-qubit system in (b) and (d) evaluate to d V max /d ∆W max = 0.35 and A V /A ∆W = 0.34, respectively.
However, for both figures of merit we see that the differences between the two protocols become less pronounced, and partially even vanish altogether, when the input energy reaches sufficiently large values.Therefore, if we want to almost fully charge the quantum battery, there is no apparent priority to use one of the mentioned protocols rather than the other.We attribute the latter convergence of the protocols to the fact that the finite system dimension severely limits the possible options for protocols to differ when the input energy is large.
In our discussion up to this point, we have considered unitary operations that act globally on the entire Hilbert space for both types of considered systems.However, in particular for the N -qubit case, one might argue that non-local operations, i.e., operations that act jointly on several (or all) qubits and may entangle them, might be considerably more difficult to implement.We therefore briefly want to examine the role of local operations for the task at hand in the next section.

LOCAL VS NON-LOCAL OPERATIONS
Here we are interested in investigating the role of nonlocal operations in charging multipartite quantum batteries.In particular, we focus on the work fluctuations and the charging precision when unitarily transferring energy into a system comprising N non-interacting qubits with equal energy gaps.To do so, we compare the generically non-local processes that arise from applying the protocols considered in Secs. 3 and 4 to N -qubits with two alternatives: first, with a specific local process that we refer to as symmetric local charging, and second, with the numerically obtained optimal local processes for work fluctuations and the charging precision.
The motivation for this comparison are twofold: Apart from the observation that local operations might be easier to implement than non-local ones, the fact that local operations do not create any correlations between the qubits can help us to better understand the role of correlations in charging processes.Meanwhile, it is expected that, much like in the reverse process of (unitary) work extraction [65], the increase of the average energy via local unitaries performs worse than the corresponding global protocols.
Local charging processes.To study the charging process in terms of local operations, we consider the N -qubit system as a many-body quantum battery whose energy increase is brought about via local unitary operations U loc , i.e., U loc = ⊗ N i=1 U i , where U i acts unitarily on the Hilbert space of the ith qubit.If we assume that the initial state of the total system is an uncorrelated thermal state at inverse temperature β, then there also is no correlation in the final state.We then have where ϱ i = U i τ (β) U † i .Due to the identical local Hamiltonians, the work deposited in the single-qubit batteries through these operations can be obtained as In this case, since there are neither classical nor quantum correlations present, one can easily show that the variance with respect to the local Hamiltonians is given by the sum of the local variances, Similarly, it is straightforward to show in this case that the work fluctuation of the total system can be written as a sum of the local work fluctuations, Symmetric local charging.Now we specialize our discussion to a specific type of local charging process that we call symmetric local charging process (SLCP).In such a process, for a given total energy increase ∆ϵ, the energy of each of the N qubits increases by a fixed amount ∆ϵ N via some fixed local unitary, i.e., U i = U for all i.For general systems, e.g., harmonic oscillators, such an SLCP can result in variances and fluctuations that can be larger or smaller than those of other local charging processes at the same energy input, see, e.g., [46,Fig. 4].Here, however, we show that SLCPs for N qubits generate the maximal amount of work fluctuations and variance for a given energy increase among all local unitary charging processes.We thus find that the worst-case local scenarios for both quantities of interest are realized by the same protocol, see Appendix A.
However, this worst-case local process needs to be compared to two different optimal (local/global) protocols.As optimal global procedures we consider the protocols introduced in the previous sections of this article.For the optimal local strategies, we numerically determine the minimum variances and work fluctuations.As shown in Fig. 7 (a) and (b) for a 4qubit system, SLCPs generally do not coincide with either the global or local optimal protocols, even when the initial state is an energy eigenstate (the ground state, in the scenario we consider).The exceptions we observe are only the trivial cases of zero charge and maximal charge.In the considered temperature range, both the variance and work fluctuations obtained from SLCPs therefore maximize the difference to the respective quantities from both the optimal local and the optimal global processes.
For low initial temperatures, there does not appear to be a discernible difference between applying the optimal local or global operations, with perfect matching for the ground state.However, for increasing temperatures we observe that the gap between the optimal local strategy and the optimal global strategy increases, while the gap between the optimal local process and the SLCP decreases.
Regarding the correlation in the final state, we note that there can be many global unitaries generating the same probability distribution (on the diagonal of the density operator) for a given initial state.In particular, this distribution can arise from correlated or uncorrelated states, and which is the case depends both on the unitary and on the initial state.In addition, degeneracy in the energy basis results in a set of different optimal distributions for a given energy input.All of these factor suggest that relation between the specific correlations of the final state and achieved figure of merit are both difficult to determine and likely of no practical concern.Nevertheless it is clear that the advantage of the non-local operations comes from the fact that the set of probability distributions in the energy basis that is are reachable via such unitary operations is larger than the corresponding set for local operations.
Finally, let us remark that we do not know if the local worst case also represents the overall worst case among all (also global unitary) processes.For the case of harmonic oscillators this is trivially the case because the variances and fluctuations can diverge for both the local and global processes due to the infinite Hilbert-space dimensions of each individual os- (both in units of ω 2 ) obtained from the protocols (conjectured to) minimize the respective quantities are plotted as solid lines in (a) and (b), respectively, as functions of the energy input ∆ϵ (in units of ω) of a 4-qubit system for temperatures (in units of ̵ hω/kB) from T = 0.2 (blue) to T = 1 (red) in steps 0.2.The dot-dashed and dashed lines represent the respective quantities obtained from the local charging processes, and are the analytically worst-case and numerically best-case local processes for both of the quantities, respectively.cillators.At the same time we do not know the optimal local strategies for N qubits for either precision or fluctuations.Determining such strategies would require an optimization over all ways of splitting the energy contributions among the considered qubits.We leave both of these questions as open problems for future work.

CONCLUSION
In this article we have investigated the process of battery charging, i.e., depositing work, for finite-dimensional quantum systems.Starting from thermal states with no extractable work, we have considered processes that raise the average energy unitarily so that all deposited energy can be extracted again in principle, and we have focused our discussion on two figures of merit that characterise such a process: charging precision and work fluctuations.We have further centered our investigation on two exemplary systems of interest: ddimensional quantum systems with equally spaced Hamiltonians, and systems of N non-interacting qubits with identical energy gaps.
For these systems we designed protocols with the purpose of optimizing the charging precision, that is, minimizing the variance of the final average energy, and minimizing the work fluctuations during the charging process, respectively.While we show that our protocol for minimizing the precision is indeed optimal, a similar proof of optimality of the competing protocol for the work fluctuations remains elusive.Yet, all evidence we have gathered seems to suggest that the protocol is indeed optimal, in particular, it reduces to the optimal protocol known for harmonic oscillators [46] in the limit d → ∞.However, except for the notable case of initial zero temperature or the special case where the dimension of the qudit diverges, the two protocols generally differ and so optimizing with respect to one figure of merit comes at the expense of suboptimal performance in the other.We have therefore compared the performance of the two protocols with respect to both figures of merit.Based on the evidence we have currently avail-able, we see no fundamental reason to generally pick one of these protocols over the other.Although there are some isolated parameter regimes (see Fig. 6) where the two protocols give the same performance for one of the two quantities (but not the other), in general one of the protocols outperforms the other.Consequently, there is a trade-off between these quantities: Optimizing performance with respect to precision means sub-optimal fluctuations and vice versa.The choice of optimal protocol therefore strongly depends on the weighting one assigns to the two figures of merit, the specific system (dimension, composition, Hamiltonian) under consideration, and on the choice of input energy.
Another question that comes into play in the N -qubit scenario concerns the potential added complication of requiring non-local operations acting jointly on all qubits in order to achieve optimal performance.To illustrate this problem we have compared our protocols to a simple local charging scenario that requires only identical unitary operations to be performed individually on all N qubits.We show that such an approach results in the worst possible performance, which highlights that access and control over global unitary transformations is another resource that needs to be considered in this context (see, e.g., the discussion in [7]).We have also numerically compared both protocols with the optimal local protocol.We numerically show the advantage of non-local operations in the charging process when the system is initially in a state far from the ground state.
Meanwhile, the energetic correspondence between the two systems considered, N qubits versus a single qudit of dimension N + 1 with matching energy-level spacing, provides an opportunity to examine the role of the internal structure and the system dimension, 2 N versus N + 1, in determining advantageous charging strategies.We show that for given initial energy and energy input, accessing higher Hilbert-space dimensions (N qubits) allows us to achieve smaller fluctuations and higher precision as compared to lower Hilbert-space dimensions (N + 1-dimensional qudit).
For future work we envisage an even more allencompassing approach towards studying the trade-offs between the identified resources in achieving optimal or nearoptimal performance in terms of the figures of merit considered here, but also beyond, taking into account such aspects as the charging speed and stabilisation of the charge.
To determine the stationary points which are related to the best-case or worst-case scenario, we calculate the first partial derivatives of L F with respect to the variables p1i , The condition of vanishing first derivative at the extremal points yields p1i = ω p 0 (λ F − ω)(p 0 − p 1 ) ∀ i, (A17) which tells us that to reach the extremal point, the energy of the system should be increased via symmetric local unitary operations.By comparing the symmetric and asymmetric lo-cal transformations described in Appendix A.1, one then arrives at the same conclusion as before, that is, the obtained extremal point represents the local unitary process achieving maximal work fluctuations at fixed energy input.Thus, if we want to locally raise the energy of the N -qubit system by the amount ∆ϵ tot , the largest fluctuations are obtained when each of the qubits reaches a state with equal energy ε = ϵ 0 + ∆ϵ tot /N .We thus see that for a given energy ε, the obtained unique extremal point is given by p1i = ε/ω which results in the maximal possible work fluctuations for locally charging an N -qubit system.Consequently, we have obtained the local-unitary charging protocol with the largest fluctuations and hence with the worst performance.Moreover, the local unitary worst-case scenarios of work fluctuation and variance are compatible and match for any initial temperature in such a process.

d , where 1 d
denotes the ddimensional identity matrix, and H d = ∑ d−1 n=0 E n |n⟩⟨n| is the Hamiltonian of a d-dimensional system with equally spaced energy levels, i.e., E n = nω.

FIG. 3 .
FIG.3.Optimal precision protocol.The minimal energy variance V (in units of ω 2 ) obtained via the optimal-precision protocol is shown as a function of the energy input ∆ϵ for (a) a 5-dimensional quantum system with equally spaced energy levels, and (b) a system of 4 identical qubits for temperatures (in units of ̵ hω/kB) from T = 0.1 (blue, bottom) to T = 1 (red, top) in steps of 0.1.

FIG. 4 .
FIG. 4. Ansatz for minimizing work fluctuations.To illustrate the working principle of the proposed protocol to suppress work fluctuations during the charging process as much as possible, the initial and final probability distributions corresponding to the diagonal of the battery-system density operator with respect to the energy eigenbasis are shown for increasing the energy of the system by ∆ε(m, k).

FIG. 5 .
FIG.5.Proposed protocol for minimizing work fluctuation.The curves show the work fluctuations (∆W ) 2 (in units of ω 2 ) obtained from the protocol proposed in Sec. 4 as a function of the input energy ∆ϵ (in units of ω) for (a) a 5-dimensional quantum system, and (b) for a 4-qubit system, for temperatures (in units of ̵ hω/kB) from T = 0.2 (blue) to T = 1 (red) in steps 0.2.For the chosen parameters (input energies, temperatures, and system dimension/qubit number) results from a numerical optimization of the work fluctuations are included but cannot be distinguished by the naked eye from the curves resulting from our proposed protocol, supporting the conclusion that the latter is at least a close approximation of the true optimum.

2 FIG. 6 .
FIG. 6.Comparison of precision and work fluctuations in both protocols.We compare the protocol that minimizes the variance (dashed lines) and the proposed protocol aiming to minimize the work fluctuations (solid curves) by calculating both the variance [upper panels, (a) and (b)] and the work fluctuations [lower panels, (c) and (d)] in units of ω 2 , for both protocols, for both a 5-dimensional [(a) and (c)] and for a 4-qubit [(b) and (d)] system as functions of the energy input ∆ϵ (in units of ω) for temperatures (in units of ̵ hω/kB) for T = 0.1 (blue), T = 0.5 (purple), and T = 1 (red).
between the function values in (a) and the difference d ∆W max ∶= max ∆ϵ ((∆W opt.V ) 2 − (∆W opt.∆W ) 2 ) between the function values in (c) is d V max /d ∆W max = 0.48 < 1, and the ratio of the areas enclosed between the respective curves [shaded red areas in (a) and (c)] is A V /A ∆W = 0.51 < 1.

2 FIG. 7 .
FIG. 7.Local vs non-local operations.The variance V and work fluctuations (∆W ) 2 (both in units of ω 2 ) obtained from the protocols (conjectured to) minimize the respective quantities are plotted as solid lines in (a) and (b), respectively, as functions of the energy input ∆ϵ (in units of ω) of a 4-qubit system for temperatures (in units of ̵ hω/kB) from T = 0.2 (blue) to T = 1 (red) in steps 0.2.The dot-dashed and dashed lines represent the respective quantities obtained from the local charging processes, and are the analytically worst-case and numerically best-case local processes for both of the quantities, respectively.