Precision and Work Fluctuations in Gaussian Battery Charging

One of the most fundamental tasks in quantum thermodynamics is extracting energy from one system and subsequently storing this energy in an appropriate battery. Both of these steps, work extraction and charging, can be viewed as cyclic Hamiltonian processes acting on individual quantum systems. Interestingly, so-called passive states exist, whose energy cannot be lowered by unitary operations, but it is safe to assume that the energy of any not fully charged battery may be increased unitarily. However, unitaries raising the average energy by the same amount may differ in qualities such as their precision, fluctuations, and charging power. Moreover, some unitaries may be extremely difficult to realize in practice. It is hence of crucial importance to understand the qualities that can be expected from practically implementable transformations. Here, we consider the limitations on charging batteries when restricting to the feasibly realizable family of Gaussian unitaries. We derive optimal protocols for general unitary operations as well as for the restriction to easier implementable Gaussian unitaries. We find that practical Gaussian battery charging, while performing significantly less well than is possible in principle, still offers asymptotically vanishing relative charge variances and fluctuations.


Introduction
Quantum thermodynamics (QT) deals with the manipulation and transfer of energy and entropy at the quantum scale. How well one can transfer energy depends greatly on the information one has about a system [1,2,3]. Consequently, the system entropy quantifying this information is rendered an important quantity for achievable state transformations [4]. At fixed energy, the entropy is maximized for thermal states, which allows for the definition of thermal equilibrium characterized by the emergent notion of temperature.
A system in such a thermal equilibrium with an environment at temperature T is thermodynamically useless in the sense that its energy cannot be extracted as work [5,6]. Therefore, much effort has been invested into understanding the emergence of equilibration and thermalization in quantum systems [7]. At the same time, quantifying extractable energy and identifying achievable transformations crucially depends on the control one assumes to have about microscopic degrees of freedom. For instance, acting only upon individual quantum systems from whom work is to be extracted gives rise to the notion of passive states [8] which cannot yield any work in cyclic Hamiltonian processes, even if the entropy at a given energy is far below the thermal entropy [9]. However, even non-passive states may still require complex operations and precise control over large Hilbert spaces that make them practically unfeasible sources of work. A recent focus of thermodynamic resource theories has thus been to investigate the role of precise control and practically implementable operations for achieving desired work extraction [10,11,12,13] and refrigeration [14].
The resource-theoretic view on quantum thermodynamics of course extends beyond the task of work extraction, and generally aims to identify the ultimate limitations of all single-shot processes [5,15,16]. More specifically, viewing quantum thermodynamics as a resource theory entails either the ability to perform any unitary operation induced by a Hamiltonian H(t) with control parameter t (i.e., the case of "driven" or controlled operations), or applying arbitrary "thermal operations", i.e., global energy conserving operations on the chosen system and arbitrary many auxiliary systems. When these auxiliary systems can feature coherence w.r.t. the energy eigenbasis (i.e., if coherent "batteries" are provided), these paradigms become equivalent [17,18]. However, neither paradigm limits the complexity of the allowed operations, requiring arbi-trary coherent energy shifts in subsystems from which energy is extracted or in which energy is stored [19,20]. This can lead to genuine quantum advantages, e.g., for the charging power of N -qubit batteries [21,22] and for small, finite-dimensional systems (e.g., few-qubit registers) such full control over the quantum systems may reasonably be expected. However, for larger systems such as registers of many qubits, arbitrary global operations may be difficult to realize and call for more specialized practical solutions [23]. In particular, this applies to infinite-dimensional quantum systems such as (ensembles of) harmonic oscillators. Besides the paradigmatic two-dimensional Hilbert spaces of qubits often favoured in information-theoretic approaches to quantum mechanics, harmonic oscillators play a crucial role for the description of physical systems in quantum optics and quantum field theory. Indeed, all current realistic proposals for and implementations of quantum machines involve at least one Hilbert space corresponding to a harmonic oscillator. Examples for such systems include superconducting resonators [24,25], modes of the electromagnetic field in a cavity [26], or vibrational modes of trapped ions [27,28]. It is hence of conceptual significance to understand the fundamental as well as the practical limitations for thermodynamic tasks in such infinite-dimensional continuous-variable (CV) systems. In particular, full control over such systems automatically implies the ability to create coherence between energy levels with arbitrarily large separation.
In contrast, a class of operations that can typically be realized comparatively simply in quantum optical realizations of CV systems is that of Gaussian unitaries [29]. In the context of driven quantum systems these operations naturally appear as most straightforwardly implementable in the hierarchy of driving Hamiltonians since they require H(t) to be at most quadratic in the system's creation and annihilation operators. Indeed, most natural interactions are appropriately described by such "bipartite terms" (usually resonant energy exchanges), whereas the creation of higher order terms is a challenge that is usually addressed only in a perturbative way. For the task of work extraction, the restriction of thermodynamic operations on CV systems to Gaussian transformations brings about the notion of Gaussian passivity [10], which encompasses states that are potentially non-passive, but are passive w.r.t. Gaussian transformations. Once work has been extracted, one would of course also like to put it to use, potentially at a later time. This requires the previously obtained energy to be stored and distributed. It is hence expected that practical limitations applying to work extraction -in particular, the restriction to Gaussian operations -will also be relevant for these tasks. sivity [10], we hence aim to quantify the limitations imposed by the restriction to Gaussian unitaries on the task of energy transfer to suitable quantum optical storage devices, i.e., charging batteries. More precisely, we consider ensembles of harmonic oscillators as batteries. These batteries are assumed to be initially uncharged in the sense that they contain no extractable work. That is, we consider the empty batteries to be in thermal equilibrium with the environment, and describe them by thermal states at the ambient temperature. We then study the task of unitarily increasing their energy by a fixed increment ∆E. Although such unitaries always exist in infinite-dimensional Hilbert spaces, unitaries achieving a given energy increase are not uniquely determined by ∆E, and may offer different charging precision, speed, and energy fluctuations during the charging process.
Here, we focus on two quantities characterizing the reliability of the charging process, as illustrated in Fig. 1. First, the charging precision, represented by the energy variance V of the battery, which is of interest since it is desirable that a charged battery is able to deliver the expected energy, not hazardously much more energy or disappointingly much less. Second, we consider the fluctuations during the charging process, captured by (∆W ) 2 the average square deviation of the energy transitions from the average energy supply. While the variance quantifies the usefulness of the battery in terms of potential fluctuations occurring when discharging the loaded battery, the variance loosely speaking only captures half of the problem. That is, taking into account the initial distribution of energies one may also be interested in the energy fluctuations during the charging process. The resulting distribution is often called "fluctuating work" [30,31,32] and characterizes the distribution of work if one were to measure the battery in the energy eigenbasis at the beginning and end of the charging protocol. For both of these characteristics we determine the ultimate limitations during arbitrary unitary charging processes by designing optimal protocols. We then specialize to Gaussian unitaries, for which we identify the optimal and worst charging protocols. In comparison, we find that Gaussian unitaries perform significantly less well than is possible in principle. Nonetheless, Gaussian battery charging can asymptotically achieve vanishing relative fluctuations V /∆E and (∆W ) 2 /∆E for large input energies by way of simple combinations of displacements and single-mode squeezing. Our results hence provide insights into both the fundamental and practical limitations of charging quantum optical batteries.
This article is structured as follows. In Sec. 2, we set the stage for the investigation and define the quantities of interest. We then present an investigation of the fundamental limitation of charging quantum batteries using arbitrary unitaries in Sec. 3, before we restrict to Gaussian transformations in Sec. 4. Finally, we draw conclusions in Sec. 5.

Charging a quantum battery
As battery systems to be charged we consider a number of bosonic modes (i.e., an ensemble of harmonic oscillators) initially in thermal states τ (β) = exp(−βH)/Z, where Z = Tr exp(−βH) is the partition function, β = 1/T is the inverse temperature of the battery (we use units where = k B = 1 throughout), and H = j ω j a † j a j is the system Hamiltonian. The mode operators a j and a † j satisfy the usual commutation relations [ a j , a † k ] = δ jk and [ a j , a k ] = 0. For such a noninteracting Hamiltonian, the initial state is a product state τ (β) = i τ i (β). The single-mode Gibbs states τ i (β) can be written as with respect to their respective Fock bases {| n i }, where a i | n i = √ n i | (n − 1) i and a † i | n i = (n + 1) i | (n + 1) i . This choice of initial state ensures that the batteries are truly empty at first, i.e., the initial state is passive for any number of such batteries because the Gibbs state is completely passive (uniquely at fixed energy).
We are then interested in applying a unitary transformation U to raise the average energy by ∆E, transforming the initial state τ (β) to a final state ρ = U τ U † , i.e., We quantify the charging precision via the increase of the standard deviation of the system Hamiltonian, that is, one of the quantities that we are interested in is where the variance w.r.t. H is given by Besides the precision of the final battery charge, one may also care about other quantities, for instance, the energy fluctuations 1 of the charging process. That is, we consider the average squared deviation from the average energy increase, given by where W m→n = E n − E m is the work relating the m-th and n-th energy levels, with H | n = E n | n , and is the probability of a transition from the m-th to the n-th energy eigenstate starting from the initial state τ (β) with diagonal elements p n = n | τ | n . To better understand the quantity ∆W it is useful to note that we can write the squared work fluctuation as the variance of the operator H ∆ =H −H in the thermal state, wherẽ H = U † HU , i.e., where the covariance is given by and {H, H} + =HH + HH denotes the anticommutator. In general, the operatorsH and H need not commute, but since the initial thermal state is diagonal in the energy eigenbasis, we can further simplify Eq. (7) and obtain The squared increase of the standard deviation, in comparison, can be written as Since one can write E(τ )E(ρ) = Tr(U † HU H τ τ ), it is easy to see the charging precision and fluctuations coincide when the initial state is an eigenstate of the Hamiltonian, because H τ τ = E n | n n | = Hτ and V (τ ) = 0. In this case (which, in our scenario only occurs for the ground state since our initial state is a thermal state), one has (∆W ) 2 = (∆σ) 2 = V (ρ).

Fundamental Limits for Battery Charging
In this section, we investigate the fundamental limits on the charging precision and fluctuations. As we shall see, optimal protocols can be constructed that minimize either the variance of the final energy or the fluctuations during the charging process, but these do not coincide for finite temperatures. However, the involved operations are often rather complicated in the sense that they require very specific interventions in particular subspaces of the infinite-dimensional Hilbert space, tailored to the initial temperature and energy supply. The results obtained in this section hence illustrate what is in principle possible and provide a benchmark for the precision and fluctuations achievable with Gaussian unitaries.

3.A Fundamental limits for zero temperature
Let us first consider a simple example to set the stage for a further, in-depth investigation. To this end, we consider a single-mode battery that is initially in the ground state, i.e., H = ωa † a and τ = | 0 0 |. In this case, the work fluctuations and charging precision coincide and are given by and ρ = U | 0 0 | U † = | ψ ψ | is a pure state. Since the Hilbert space in question is infinite-dimensional the energy variance of the final state is not bounded from above. This can be seen by choosing a super- In other words, for any chosen energy ∆E one can make k (and hence ∆σ = ∆W ) arbitrarily large by simultaneously choosing q sufficiently close to 1. So for arbitrarily  The maximal and minimal variances of the energy that are in principle possible for a battery that starts in its ground state and is being charged by ∆E are shown (in units of ω, with = 1) for a system of dimension d = 6. When the Hilbert space is infinite-dimensional, the lower bound periodically repeats, but the upper bound is no longer finite for any value of ∆E. Since the initial temperature vanishes, the bounds shown also apply to the charging fluctuation ∆W .
small energies, the energy variance and the fluctuations during the charging process may increase by an arbitrary amount. However, it is also clear that this is an artefact of the infinite-dimensional character of the system. If the dimension d of the system is finite (or there is some cutoff energy), then the maximal variance is obtained for a superposition of the eigenstates | 0 and | d − 1 , with minimal and maximal eigenvalues, respectively, resulting in The minimal achievable variance for any given energy is obtained by unitarily rotating to a superposition of the two energy eigenstates | n and | n + 1 that are closest to the available energy, i.e., such that n ≤ ∆ ≤ n + 1.
resulting, after some algebra, in a variance of Crucially, ∆σ min = 0 whenever ∆E is an integer multiple of the oscillator frequency, and the maximal value of ∆σ min is ω 2 , as illustrated in Fig. 2.

3.B Fundamental precision limits for arbitrary temperatures
Having understood the simple case of optimally charging a battery initially in the ground state, we now want to move on to the case of thermal battery states. On the one hand, the worst-case scenario immediately carries over from the situation discussed in the previous section. That is, in an infinite-dimensional system one may always find a unitary transformation that increases the energy by an arbitrarily small amount, while increasing the variance arbitrarily strongly. This can be seen by just noting that the two-level rotation used to rotate between the ground state and the level | k can also be applied to thermal states. The only difference is that the corresponding probability weights are now different from 1 and 0 initially. In contrast, the upper bound for the variance in a finite-dimensional system is always finite.
The optimally achievable charging precision, on the other hand, requires a more intricate analysis. The task at hand is to specify the state of minimal energy variance V (ρ) at a fixed average energy within the unitary orbit of a thermal state at a given temperature. In general, we cannot give a closed expression relating this minimal variance to the energy input and the initial temperature. However, one may formulate a protocol that provides (one of) these minimal variance states. Here, we will give a short, intuitive description of this protocol, and provide a detailed step-by-step account in Appendix A.1.
Let us now briefly explain the working principle of the optimal-precision charging protocol. First, recall that the initial thermal state has a density operator τ (β) that is diagonal in the energy eigenbasis with probability weights p n decreasing with increasing energies E n . The average energy E τ (β) is determined by the initial temperature and we hence know the target energy E(ρ) = E(τ ) + ∆E for any energy input. We can then naively apply two-level rotations to reorder the probability weights on the diagonal such that the largest weight p 0 is shifted to the eigenstate whose energy is closest to the target energy, the second largest weight is shifted to the second-closest eigenstate to E(ρ), and so on. This procedure results in the unique stateρ(β) within the unitary orbit of τ (β) whose average squared deviationṼ from the target energy is minimal.
Unfortunately, this state does not generally have the desired target average energy, i.e.,Ẽ = E(ρ) = E(ρ). Consequently, the average squared deviation from E(ρ) is generally not equal to the energy-variance,Ṽ = V . Moreover, both the casesẼ > E andẼ < E can occur and one therefore has to adjust the energy accordingly. This can be done by sequences of two-level rotations that change the energy by ∆Ẽ and increase the average squared deviation from E by ∆Ṽ . An ordering of these operations that is optimal is obtained when performing them in the order of increasing values of ∆Ṽ /|∆Ẽ|, starting with the smallest, i.e., when the increase ofṼ  per unit energy change is as small as possible. One carries on with this protocol until the desired target energy is reached, in which case the final value ofṼ becomes the variance of the energy V . The resulting variances for given energy input for a harmonic oscillator are illustrated in Fig. 3. It can be seen that for higher initial temperatures this optimal protocol can lead to decreasing variances in the battery state. Also note that the working principle of this optimal protocol is unchanged if one considers a finite-dimensional system instead and differences only arise because of the finite maximal energy input at any given temperature.

3.C Fundamental precision limits for multi-mode batteries
After obtaining the fundamental limits on the precision of charging a single-mode battery, it is of course natural to ask which possibilities arise when several such batteries are available. The worst case scenario for multiple modes trivially translates from our previous anal-ysis. Since the variance for any given energy input is not bounded from above for single-mode batteries, the same is also true for many modes.
To understand what can be achieved in the best case for multiple batteries, let us first consider the two-mode case, i.e., two batteries labelled A and B that are initially in a thermal state τ A ⊗ τ B . We are now interested in an increase of the average energy E(ρ) = Tr ρH w.r.t. E(τ ), where the bipartite Hamiltonian is H = The energy variance is then given by (16) where the covariance of Eq. (8) for the local (and hence commuting) operators H A and H B is For a local unitary charging protocol, i.e., where U = U A ⊗ U B , the initial thermal states remain uncorrelated and the covariance vanishes. That is, the final state ρ AB is a product state ρ AB = ρ A ⊗ρ B . In such a case not only the average energies but also the variances are additive. Inspection of Fig. 3 then shows that having two or more batteries available can be beneficial even when they are charged independently. For instance, when the supplied energy would lead to a local maximum of the variance if all energy is stored in one battery, it may be prudent to reduce the energy supply to this battery to reach a (local) minimum instead. The remnant energy can then be stored in a second battery. When the two modes have the same frequency and the initial temperature is nonzero the resulting overall variance is then smaller than or equal to that of charging only one battery, as we can see from Fig. 4. In short, the availability of several battery modes at potentially different frequencies hence provides a certain flexibility to reach local minima of the variances of the individual batteries, but the exact performance for a given set of battery modes requires to be worked out on a case-by-case basis.
For unitaries that are not local and can correlate the two batteries, the situation is even more involved but in principle such unitaries may help to achieve an even better performance. To see this, let us return to the optimal protocol of the last section. In the first step of this protocol, the probability weights of the initial thermal state are reordered to create a distribution that is as narrow as possible around the target energy. The resulting state is diagonal in the energy eigenbasis. Since this is a product basis w.r.t the tensor product structure of different modes, the state is still uncorrelated. However, in the second step, where the energy of the distribution is adjusted to the target energy, two-level rotations with optimal ratios ∆Ṽ /|∆Ẽ| may occur between states | m, n and | m , n with m = m and n = n and hence correlate the systems. For batteries at different frequencies there can thus be an advantage in introducing (specific) correlations, whereas a situation as just described can always be avoided for batteries with equal frequencies. In the general case of arbitrary frequencies it is interesting to note though, that the creation of correlations may be marginally helpful but is not the key ingredient. This is in contrast to recent results on the charging power, where the ability to create quantum correlations, i.e., access to entangling operations (albeit not necessarily the actual creation of entanglement) can be extremely useful [21,22].
To reach optimality it nonetheless remains to be determined how the energy can be optimally split between the oscillators, or invested in correlations. Unfortunately, this is difficult to answer in general, and is even rather complicated for uncorrelated charging due to the non-monotonic behaviour of the optimal single-mode charging protocol illustrated in Fig. 3, which is illustrated in Fig. 4. There, the specific optimal splitting depends on the initial temperature, the specific energy input, and the (number and) frequencies of the battery modes involved. The optimal performance hence has to be determined on a case-by-case basis. However, one can state quite generally that the optimal final variance of the joint system is never larger than the optimal variance when all the energy is stored only in one of the modes. In other words, having several battery modes available is never detrimental. Indeed, having more empty batteries at different frequencies at one's disposal can be considered a nontrivial resource for precise charging.
Having discussed which charging precisions can be achieved in principle, let us briefly turn to the fundamental limitations arising for the charging fluctuations.

3.D Fundamental fluctuation limits for arbitrary temperatures
To complete the investigation of the fundamental restrictions of charging a quantum battery, let us consider a protocol that minimizes the fluctuations ∆W . For simplicity, let us start with the case where the input energy is exactly one unit, ∆ = 1. Then the infinitedimensional Hilbert space allows keeping the fluctuations arbitrarily small. To achieve this, we perform a unitary permutation operation on the first N energy levels that shifts the weight p n = (1 − e −βω )e −nβω from the level n to the level n + 1 for n = 0, . . . , N − 1, while the last weight p N is shifted to the ground state level. In the limit N → ∞, the energy is increased by ∆E = ω and since m | U | n = δ m,n+1 , the fluctuations vanish.
When the input energy is less than one unit, i.e., when 0 < ∆ < 1, the fluctuations do not vanish, but can be minimized in a simple way. Suppose that we perform the same permutation as before, but start shifting weights upwards at some finite n = k rather than at n = 0, such that a vanishingly small weight is placed on the k-th level. The corresponding final state energy (in units of ω) would bẽ Now, generally, (βω) −1 ln(1/∆ ) is not an integer and, consequently, the energy shift upwards starting from k =k := (βω) −1 ln(1/∆ ) is not enough, ∆ I := e −kβω ≤ ∆ . However, if we perform the shift from k onwards nonetheless, the difference ∆ II = ∆ − ∆ I can be obtained by continuously rotating between the levelk − 1 and the (now effectively unoccupied) levelk, i.e., by a mapping The corresponding rotation angle θ is given by This protocol is optimal since each (finite size) weight is shifted by either 0 or 1 units of energy, i.e., the shifts closest to ∆ since 0 ≤ ∆ ≤ 1. Explicitly, we can calculate the corresponding fluctuations by splitting the contributions for the differently shifted weights, i.e., For where we have used (20). The remaining shifts fromk upwards give rise to When summing up these contributions, substituting sin 2 θ = ∆ II /pk −1 from Eq. (21), and noting that ∆ = ∆ I + ∆ II , we find ∆Wmin ω for 0 ≤ ∆ ≤ 1. Finally, consider the case where ∆ > 1. Then we perform the protocol just described, but replace ∆ with the difference ∆ − ∆ to the lower integer value. The remaining energy is now an integer multiple of ω and can be gained by shifting the entire distribution upwards by ∆ units, whilst filling the gaps with vanishing contributions from arbitrarily high levels (as described for ∆ = 1 at the beginning of this section). Since the last integer shift does not add any fluctuations, we arrive at the optimal value ∆Wmin ω Note that this expression for the minimal fluctuations at arbitrary temperatures coincides with the expression for the minimal variance (∆σ min /ω) 2 achievable at zero temperature, as given in Eq. (15) and illustrated (by the lower curve) in Fig. 3, but for finite temperatures the protocol minimizing the fluctuations does not minimize the variance, and vice versa. As a remark, note that in contrast to the optimal precision protocol, the protocol for minimal fluctuations does not translate directly to the finite-dimensional case.
As for the case of the variance, let us now turn to the case of several modes, starting with two. Here it is first important to note that a second battery can be added without increasing the fluctuations since for local unitaries one finds where ∆E = ∆E A + ∆E B . Second, one may note that the protocol described above can now achieve vanishing fluctuations also for energies ∆E = mω A + nω B for m, n ∈ N 0 , not just for integer multiples of a single frequency. In addition, the optimization of the energy splitting between the two modes can lead to lower fluctuations as compared to only charging one of the batteries also for energy values that lie in between two choices of m and n. All of this can be done using only local unitary charging. Correlating unitaries again only play a minor role in the sense that they may be employed in optimizing the second part of the protocol, where the missing energy ∆ II is added. The presence of multiple modes as batteries to be charged can hence be considered to be helpful. However, as before, the exact optimal protocols for multiple modes depend on the respective frequencies, temperatures, and on the input energy, and hence require case-by-case analyses. It thus becomes ever more clear that the operations to optimize either the variance or fluctuations are generally complicated and require extreme levels of control over the infinite-dimensional systems we consider here. It is hence of great interest to turn to practical operations such as Gaussian unitaries, and investigate their limitations for realistic battery charging.

Battery Charging Using Gaussian Unitaries 4.A Preliminaries: Phase space and Gaussian states
In the following, we want to study the restrictions imposed on the battery charging scenario when only Gaussian unitaries are used, i.e., unitary operations that map Gaussian states to Gaussian states. To examine this class of states, note that any quantum state ρ in the Hilbert space L 2 (R N , dx), i.e., the space of squareintegrable (with respect to the Lebesque measure dx) functions over R N , can be assigned a Wigner function W(x, p) given by are appropriate position and momentum coordinates andx = (x 1 ,x 2 , . . . ,x N ) T and p = (p 1 ,p 2 , . . . ,p N ) T are the corresponding position and momentum operators, respectively. The eigenstates | x and | p of these operators, respectively, satisfŷ It is convenient to collect x and p into a single phase The commutation relation [ a m , a † n ] = δ mn then implies the canonical commutator [x m ,p n ] = iδ mn , and vice versa. For the Wigner function, the normalization of the density operator translates to the condition Expectation values of Hilbert space operatorsĜ can be computed from the Wigner function via where the Wigner transform g(x, p) of the operatorĜ is given by With these basic definitions at hand, we can now return to Gaussian states and operations. Gaussian states are defined as those states in H whose Wigner function is a multivariate Gaussian, i.e., of the form (35) for some vector X ∈ R 2N and a real, symmetric 2N × 2N matrix Γ. These quantities are called the first and second statistical moments of a quantum state, where the vector of first moments is simply X = X ρ and the components of the covariance matrix are given by Note that we have included a conventional factor of 2 in the definition of the covariance matrix w.r.t. the actual covariances of the operators, compare, e.g., Eq. (17). Via Eq. (35) Gaussian states are hence fully determined by X and Γ. Gaussian unitaries, which map the set of Gaussian states onto itself, are represented by affine maps (S, ξ) : Here ξ ∈ R 2N are displacements in phase space represented by the unitary Weyl operators D(ξ) = exp iX T Ωξ , which can shift the first moments, but leave the covariance matrix unchanged. The objects S are real, symplectic 2N ×2N matrices which leave the symplectic form Ω invariant, i.e., The components of Ω are given by . For more information on Gaussian operations and states see, e.g., Refs. [34,29].

4.B Charging precision for single-mode Gaussian unitaries
We now want to study the previous situation of precisely charging quantum batteries based on harmonic oscillators under the restriction to Gaussian unitaries. To this end, first note that any initial thermal state τ (β) is Gaussian for all temperatures (for the usual Hamiltonian H = n ω nNn withN n = a † n a n ). The corresponding first moments vanish, X = 0 and the covariance matrix is diagonal, where the single-mode covariance matrices are given by Γ n (β) = coth βω n /2 1 2 . In particular, when the temperature is zero, we have the ground state with Γ vac = 1, and the corresponding Wigner function To determine the energy of any Gaussian state for the Hamiltonian H = n ω n a † n a n we do not need to use Eq. (33). Inspection of the first moments X (n) = (X 2n−1 , X 2n ) T = (x n ,p n ) T of each mode and local covariances simply reveals 2 that 2 Note that there is a typographical error in the prefactor of ||X (n) || 2 in Ref. [10,Eq. (12)].
However, to compute the variance V (ρ) we require also the expectation value of H 2 . For our single-mode example (for notational convenience we drop the mode label n on all quantities from now on) we have H = ωN witĥ N = a † a. We hence need to find the Wigner transform ofN 2 . With some straightforward calculations which are shown in detail in Appendix A.2, one obtains the expression With this, we can compute the expectation value N 2 for arbitrary single-mode Gaussian states in terms of the corresponding first and second moments. Since we are dealing with a mode-local operator, we can use the single-mode version of Eq. (35) to do so, i.e., using only the vector X ∈ R 2 and the 2 × 2 covariance matrix Γ.
After some lengthy but straightforward algebra we find that for any single-mode Gaussian state Since the first term on the right-hand side of Eq. (41) With this knowledge at hand, we can now return to our problem of raising the energy of a single-mode battery using Gaussian unitaries. For single-mode batteries that are initially at a finite temperature, the initial energy and corresponding variance can be calculated from Eqs. (39) and (42) by noting that the corresponding first moments vanish, X = 0 and the covariance matrix is Γ(β) = coth βω/2 1 2 . With this we have We can then apply Gaussian unitaries to these states. For instance, we may consider single-mode displacements to raise the energy of initial thermal states. covariance matrix, the latter remains that of a singlemode thermal state, while the first moments are transformed to X i = ξ i . We hence have For displaced thermal states we consequently find i.e., an asymptotic increase of the energy standard deviation with the square-root of the energy increase. As we shall see, pure displacements are neither optimal (minimal ∆σ for given ∆E), nor the worst possible Gaussian operations for battery charging, but nonetheless, make for an interesting comparison. This is illustrated in Fig. 5, where we have also included results for the optimal and worst operations, which we shall derive next.

4.C Optimal and worst-case Gaussian charging precision
Let us now investigate these best-case and worst-case Gaussian operations. The action of an arbitrary local Gaussian unitary U G results in some (generally nonzero) first moments ξ = X(U G τ U † G ), while Γ(β) is mapped to the covariance matrixΓ of an arbitrary single-mode Gaussian state with the same mixedness via a local symplectic operation S loc , i.e.,Γ = S loc ΓS T loc . Any singlemode symplectic operation can be decomposed [35,29] into (phase) rotations R and single-mode squeezing transformations S(r) as where θ, φ are real rotation angles, r ∈ R is the squeezing parameter, and With this, the transformed covariance matrix can be written asΓ and the corresponding average energy of the state ρ = U G τ U † G evaluates to Combining this with Eq. (43), we then have the energy input For the variance of the energy of the final state we first inspect the term ξ TΓ ξ and note that for our intents we can absorb the rotation R(θ) into the choice of the first moments, since ||ξ|| 2 = ||R T (θ)ξ|| 2 . We hence find We then proceed in the following way. First, we note that ∆E is a function of r and ||ξ|| 2 , whereas V (ρ) depends on ξ only via the term ξ T S(2r)ξ = ξ 2 1 e −2r +ξ 2 2 e 2r . Since e −2r ≤ e 2r for r ≥ 0, the maximal and minimal values of V (ρ) for fixed ∆E and T = 1/β must be attained for combinations of r ≥ 0 and ξ with ξ 1 = 0 and ξ 2 = 0, respectively. Conversely, this means that the remaining quantities ξ 2 2 and ξ 2 1 in Eq. (52) can be identified with ||ξ|| 2 , i.e., from Eq. (51) we have We can thus define two functions V + (r) and V − (r), given by V±(r) which correspond to the respective restrictions of the final state variance V (ρ). Moreover, when the initial temperature and input energy are fixed these functions depend only on r, allowing one to straightforwardly determine the maxima max r V + (r) = max r V and minima min r V − (r) = min r V , respectively. As we show in detail in Appendix A.3, for every fixed ∆E ≥ 0 and T ≥ 0, there exist unique values r ± ≥ 0, such that max r V + (r) = V + (r + ) = max r V = V (r + ) and

IV.C.1 Worst-case Gaussian charging precision
In particular, we find (see Appendix A.3.I for details) that for any temperature and energy input, the maximal values of ∆σ are obtained for ξ = 0, that is, when all energy is transferred to the battery via single-mode squeezing. In this case both ∆E and ∆σ are functions of r only, and we can use Eq. (51)to relate the two quantities directly. In other words, we find along with the maximal variance increase As we see, in this case the energy standard deviation increases linearly with the energy input in the asymptotic regime (as ∆E → ∞).

IV.C.2 Optimal Gaussian charging precision
While the worst-case Gaussian unitary transformation has thus been identified as pure single-mode squeezing, the Gaussian unitary transformation that minimizes the variance of the energy can be identified as a combination of squeezing and displacement that depends on the input energy and temperature, see Appendix A.3.II. That is, in our conventions, the optimal performance is achieved for ξ 2 = 0 and generally nonzero values of ξ 1 and r = r − , where the latter is determined by the condition ∂V − /∂r| r=r − = 0, which implies Inserting Eq. (57) into V − from Eq. (54) then permits us to write Although a closed expression for r − (∆E, β) cannot be given, we show in Appendix A.3.II that r − (∆E, β) ≥ 0 exists and is unique and can thus easily be determined numerically for any given ∆E and T = 1/β via the implicit formula in Eq. (57). The value r − obtained in this way can then be inserted into Eq. (58) to arrive at the minimal variance. Results for the bounds on ∆σ for a initial thermal states are shown in Fig. 6 for a range of temperatures and input energies.
Here it is interesting to note that, while the upper bound is achieved for pure squeezing transformations, the lower bound is a combination of squeezing and displacements. In the optimal case, the energies ∆E sq and ∆E d invested into squeezing and displacement, respectively, are expressed via the optimal squeezing parameter r − as ∆Esq ω In the limit of large energy supplies, ∆E → ∞ (at fixed temperature this implies r − → ∞), we hence have That is, the energy invested into squeezing grows much less strongly than that invested into displacements. Moreover, note that while pure displacement asymptotically (i.e., for ∆E → ∞) leads to a linear scaling of the final variance with ∆E, that is, V (ρ)/∆E → ω coth( βω 2 ) as ∆E → ∞, see Eq. (46), the optimal local strategy provides a more favourable scaling behaviour even though most of the energy is invested into displacement. More specifically, considering the relative variance V − /∆E by combining Eqs. (57) and (58), and taking the limit ∆E → ∞ (corresponding to r − → ∞), one finds that V − /∆E → 0, i.e., the optimal variance scales sub-linearly with the input energy.

4.D Charging precision for multi-mode Gaussian unitaries
Let us now finally turn to the case of bounding the charging precision of a multi-mode battery under the restriction to Gaussian unitaries. As we have discussed in Section 3.C, in a general optimal protocol for the charging precision, correlations between the individual battery systems can be helpful in principle. However, this appears to be the case only if one can selectively rotate between specific energetically desirable levels. For Gaussian unitaries, such specialized operations with, in a manner of speaking, surgical precision are out of the question. In particular, one may view any multi-mode Gaussian operation as a combination of local operations and beam splitting [35]. The latter may shift the average excitation numbers between the modes but may do little more. Another aspect of introducing correlations during the charging process is that any energy stored in this way also has to be extracted globally from the joint system if optimality is to be preserved. In other words, introducing correlations raises the effective local temperatures (and hence the local entropies), reducing the local free energy. In the spirit of restricting to practical operations, we shall hence consider only local Gaussian operations from now on.
Nonetheless, it may sometimes be useful to split the energy supply between different modes in specific ways, depending on the initial temperature and energy supply. To understand why it is useful, consider the (nonoptimal) case of charging two modes labelled A and B (with frequencies ω A and ω B , respectively) via pure displacements. For such local operations, no correlations are introduced. If the energy is split in such a way that for some real p ∈ [0, 1] the energy p∆E is stored in the mode A and (1 − p)∆E in the mode B, inspection of Eq. (46) reveals that the variance of the final state ρ AB behaves as with ν i = coth βωi 2 for i = A, B. When the two modes have the same frequencies, ω A = ω B , the variance becomes independent of p, i.e., V (ρ AB ) = ν A ω A ∆E + 2V (τ A ) = ν B ω B ∆E + 2V (τ B ), and it does not matter how the energy is split. Otherwise, that is, when ω A = ω B , it becomes beneficial to store all the energy in the lower frequency mode. Now, recall that this is the case for pure displacements, which are not optimal. The optimal strategy, in contrast, provides an increase of the variance that is sub-linear with the input energy. In this case it matters how the energy is split for all frequency combinations. E.g., for ω A = ω B it becomes optimal to evenly divide the energy between the two batteries. In general, the optimal energy per battery is determined by the number and temperature of the batteries and their respective frequencies.
The worst case local scenario is obtained when the energy is used only for single-mode squeezing, where the splitting of the energy between the modes again depends on the specific situation. For instance, when the frequencies of both modes are the same, investing the same energy in both batteries via squeezing will lead to the largest variance. This is because the local variances increase stronger than linearly with the input energy for single-mode squeezing.

4.E Charging fluctuations for single-mode Gaussian unitaries
At last, let us turn to the question of bounding the possible fluctuations that may appear in Gaussian battery charging. From Eq. (9) we already know how to express the energies and variances of ρ and τ in terms of the temperature and final first and second moments. However, we still need to calculate Tr(HHτ ) for arbitrary Gaussian unitaries, where we restrict to local operations, as before. For a single mode with frequency ω we may write the term in question as 1 ω 2 Tr(HHτ ) = n p n n n | U † a † a U | n . (62) In principle, an arbitrary single-mode Gaussian unitary U G may be decomposed into a combination of singlemode squeezing operations, local rotations, and displacements, none of which commute with each other. In spite of this, for any unitary and any initial state ρ o we can find a single-mode Gaussian unitaryŨ G = D(−ξ)U G such that the first moments ofŨ G ρ oŨ † G vanish. Conversely, this means that for the purpose of calculating n | U † a † a U | n we may assume that U G may be written as U G = D(ξ)Ũ G , where D(ξ) is a pure displacement andŨ G leaves the origin of the phase space invariant. Consequently, we can use the Bloch-Messiah decomposition [35] to writeŨ G as a combination of single-mode squeezing U S (r) and local rotations R(θ), i.e., Inserting into Eq. (62) we then find that the rotations either act on the rotationally invariant Fock states (the phases cancel), or can be absorbed into the direction of the displacement (using the same symbol ξ in a slight abuse of notation). We hence find We then use the simple relations to obtain the desired matrix element We can then reinsert this result into Eq. (62) and evaluate the sum over n. Further inserting into the squared work fluctuations of Eq. (9), and combining this with the expressions for the variances and average energies from Eqs. (43), (44), and (51) we obtain where V (ρ) is given by Eq. (52). Here, we note that, apart from V (ρ), no dependency on the displacement ξ appears in Eq. (66). Consequently, we may argue as in Section 4.B, that for any given value of r, the maximal and minimal function values (here of (∆W/ω) 2 ) for fixed initial temperature and fixed ∆E are attained for combinations of (single-mode) squeezing r ≥ 0 and displacements ξ with ξ 1 = 0 and ξ 2 = 0, respectively. In other words, we are interested in determining the maximum of ∆W + (r) and the minimum of ∆W − (r), where and V ± (r) is given by Eq. (54). As we show in detail in Appendices A.3.III and A.3.IV, both extremal values exist and are unique for any given initial temperature T = 1/β and input energy ∆E ≥ 0. Once again, the corresponding extremal squeezing parametersr ± (which are in general different from the optimal squeezing parameters r ± for the charging precision) are only given implicitly, i.e., by the conditions ∂∆W ± /∂r| r=r ± = 0, which can be expressed as ∆E ω = 1 2 coth(βω/2) e ∓2r ± cosh(4r ± ) − 1 (68) respectively. Nonetheless,r ± and hence the exact optimal and worst single-mode Gaussian fluctuations can easily be obtained numerically, which we have done for some sample values shown in Fig. 7. Interestingly, the additional terms appearing in Eq. (66) besides the variances lead not only to a different optimum in terms of the relative strengths of squeezing and displacements, but also mean that the worst case is now also attained for nonzero displacements. In particular, we find that the energy ∆E ± d invested into displacement in the extremal cases is given by while the energy invested into squeezing is ∆E ± sq = ∆E − ∆E ± d . In the limit of large energy supplies, i.e., ∆E → ∞ (corresponding tor ± → ∞ at fixed temperature), we have while ∆E + sq → ∞. One thus finds that the worst-case Gaussian strategy invests almost all energy into squeezing asymptotically. Conversely, for the minimal fluctuations we have ∆E − d /∆E − sq → e −4r − asr − → ∞. In the limit of large input energies it is thus optimal to invest almost all energy into displacement to minimize the fluctuations.
The crucial feature to note is that the optimal Gaussian strategy results in a sub-linear increase of the fluctuations with the input energy. Similarly as before for the Gaussian strategy optimizing the charging precision, we can consider the limit ∆E → ∞ of the relative fluctuations ∆W 2 − /∆E. For the optimal strategy, Eq. (68) tells us that this corresponds to the limit r − → ∞, for which ∆W 2 − /∆E(r − ) → 0, since ∆W 2 − and ∆E grow as e 4r ± and e 6r ± , respectively, in this limit.
Finally, let us briefly comment on the Gaussian multimode scenario. Much like before for the charging precision, using Gaussian operations that generate correlations seems to be practically irrelevant since any energy stored in such global correlations could not be accessed locally. Nonetheless it can again be useful to split the energy in specific ways (depending on the respective frequencies) between two (or more) batteries, since the optimal protocol brings about a sublinear increase of (∆W ) 2 with the input energy.

Conclusion
In this work, we have investigated fundamental and practical limitations on the precision of charging quantum batteries and on the work fluctuations occurring during the charging process. The battery systems we consider are infinite-dimensional bosonic systems, i.e., collections of harmonic oscillators, which are paradigmatic in the theoretical description of physical systems in quantum optics and quantum field theory and are hence of high conceptual significance. We assume these systems to be initially thermalized at the ambient temperature. That is, from the point of view of a resource theory of extractable work, empty batteries are considered to be for free, as no work can be extracted from them. We find that, on the one hand, neither the fluctuations nor the precision of the charge for any finite energy input are bounded from above in principle when increasing the average energy of such batteries. On the other hand, we are able to provide lower bounds for both quantities, presenting the respective optimal protocols minimizing the energy variance or fluctuations at given energies and temperatures.
In general, these optimal protocols, though theoretically easily describable, are practically difficult to realize, since they require sequences of precise interventions in particular subspaces of the corresponding infinite-dimensional Hilbert spaces. Therefore, it is interesting to understand which limitations apply in scenarios where the energy storage is performed using practically realizable transformation. A set of operations that can usually be implemented comparably simply in such systems is the family of Gaussian unitaries. Here, we have determined the optimal and worst-case Gaussian operations for charging quantum batteries. We find that energy increase via pure single-mode squeezing is the least favourable operation if one wishes to obtain a precise charge for single-mode batteries, whereas the optimal precision, as well as the smallest and largest fluctuations within the restricted set of Gaussian operations are obtained for combinations of squeezing and displacements. For multiple modes, the situation becomes more complicated in principle, but it can be said that it is in general useful to have access to multiple batteries, and that correlations between them are not necessarily detrimental, but also only helpful indirectly.
Overall, we conclude that while the optimal Gaussian operations do not achieve results comparable in quality with optimal non-Gaussian protocols, the worst performance achieved with Gaussian operations still produces finite variances and fluctuations, whereas this is not guaranteed in general. In particular, the relative variance and relative fluctuations w.r.t. the energy input asymptotically vanish for large energy supply for the optimal Gaussian charging operations. Gaussian unitaries are hence nonetheless practically useful for battery charging. In a sense, Gaussian operations hence represent a trade-off between performance versus reliability and practicality. This is reminiscent of similar contrasts between usefulness and severe limitations of Gaussian operations for tasks in quantum information, e.g., the non-universality of Gaussian operations for quantum computation [36]. Another observation of this kind can also be made in a different quantum thermodynamical context, where Gaussian operations achieve optimal scaling for the entanglement creation for large energy inputs, but fail to create entanglement in thermal states of finite temperatures when the energy input is too small [37,38].
This work hence adds to recent efforts [10] of understanding the usefulness of Gaussian operations for quantum thermodynamical tasks, providing investigations of Gaussian unitary work extraction and energy increase. Nonetheless, future work may expand on a number of open questions. For instance, we have here mostly focused on individual batteries since any work stored in joint systems would also require joint extraction. In other words, the role of correlations for work fluctuations and charging precision may be of interest, in particular, in relation to recent results on the workcost of creating correlations [37,39,40,41,42]. In addition, it would also be of interest to consider the consequences of restricting to Gaussian operations for the charging speed (or charging power) as considered in Ref. [21,22]. Finally, we note that, while some of the results presented here (e.g., the optimal precision charging protocol) directly translate to finite-dimensional systems, other aspects of this work are applicable only to the infinite-dimensional case. An in-depth investigation of the fundamental and practical limitations of precision and fluctuations in charging finite-dimensional systems, although certainly of interest [23], goes beyond the scope of this paper.
In this appendix, we give a detailed description of a unitary battery charging protocol that raises the average energy of an initial single-mode thermal state whilst keeping the energy variance of the final state minimal. The initial thermal state with density operator τ (β) is diagonal in the energy basis. The corresponding diagonal elements are the probability weights p n = (1 − e −βω )e −nβω , which are decreasing with increasing energies E n = nω. The initial average energy E τ (β) = o ω, where o = e −βω (1 − e −βω ) −1 and the initial energy variance V τ (β) = ω 2 e −βω (1 − e −βω ) −2 are determined by the initial temperature T = 1/β. We are then interested in increasing the average energy by an amount ∆E = ω∆ to reach a state ρ with E(ρ) = ω = E(τ ) + ∆E. In particular, we aim to achieve this increase unitarily, i.e., such that ρ = U τ U † . Moreover, we want to keep the energy variance of ρ minimal. In other words, we would like to determine the (non-unique) minimal energy-variance state ρ with given average energy ω in the unitary orbit of τ (β).
The protocol that we present here to obtain such a state consists of two parts (I & II). Each of these parts can be described as a series of (unitary) two-level rotations, ensuring that the final state is within the unitary orbit of the initial state. The two-level rotations between pairs of energy levels are used to appropriately shift and reorder the probability weights p n of the initial state. As we shall explain, part I of the protocol reaches the unique stateρ in the unitary orbit of τ (β) that minimizes the average squared distance to the target energy. The stateρ is diagonal in the energy eigenbasis, and arises from a permutation of the weights p n that assigns positions with increasing distance to the target energy to weights with decreasing size. However, the state obtained in this way does not have the desired target energy, i.e., in general E(ρ) = E(ρ). During part II of the protocol, this deviation of the average energy is corrected in such a way that the target energy is reached whilst only minimally increasing the average squared distance to it.

A.1.I Part I of the protocol
In part I, we first identify the energy level (labelled k) that is closest to the desired target energy, i.e., we define where we distinguish between two cases, depending on whether is closer to the energy level above or below its value. The probability weights for the case where k = are illustrated in Fig. A.1 (a). Part I of the protocol then consists of a reordering of the weights p n such that the largest weight p 0 is moved to the energy level k, the second largest weight p 1 to the second closest level to , and so forth. After part I, the density operator is still diagonal, but the probability weights on the diagonal are now either given bỹ for n = 0, . . . , k p 2(n−k)−1 for n = k + 1, . . . , max{1, 2k} The resulting probability distribution, illustrated in Fig. A.1 (b) for the case k = , has an average energỹ I = np n n and its average squared distance from the target is minimal, i.e., we have arrived at the unique stateρ in the unitary orbit of the initial state τ that minimizesṼ I = np n (n − ) 2 . However, as˜ I in general does not match , which implies thatṼ I also is not equal to the energy variance, we are not yet done. Interestingly, for both k = and k = one may encounter combinations of T and ∆ such that˜ I < or˜ I > .

A.1.II Part II of the protocol
In part II of the protocol we hence have to appropriately adjust the average energy. This can again be done by a sequence of two-level rotations. Each of this transformations will bring the average energy closer to , but since we start from a minimum of the average squared deviation from , the value of the latter will increase. We are hence interested in selecting the optimal sequence of these two-level rotations. To start, consider a rotation between the levels m and n with weightsp m andp n by an angle θ. This corresponds to the map (p m ,p n ) → (cos 2 θp m + sin 2 θp n , cos 2 θp n + sin 2 θp m ), (A.4) and leads to a change in energy given by Meanwhile, the increase of the mean squared deviation from can be written as Now, let us pick two such values on either side of , i.e., m = k−l and n = k+l+j, where l ∈ N 0 determines the distance from the level k, and j ∈ Z quantifies the difference in distances to k (or equivalently to ) between the levels m and n. More specifically, we have Moreover, this implies that the energy change from Eq. (A.5) is given by With this we further find that the changes of the energy and ofṼ have the relation With this knowledge, we come to a more detailed description of the protocol. First, we set˜ =˜ I and V =Ṽ I , and we distinguish between the two situations < and˜ > . On the one hand, when˜ < , we need to increase the energy, which means picking levels m and n such thatp m >p n and 2l + j > 0, whilst choosing j as small as possible to minimize ∆Ṽ ∆˜ . On the other hand, when˜ > , we need to decrease the energy, suggesting that one should select levels m and n such thatp m <p n and 2l+j > 0, whilst choosing j as large as possible to minimize ∆Ṽ ∆˜ . When such levels are chosen and the potential energy change exceeds what is needed to reach the target, one appropriately fixes the rotation angle θ such that˜ + ∆˜ = . If the energy change achievable with a specific such rotation is not sufficient to reach the target, one rotates by θ = π 2 , up-dates˜ ,Ṽ , and {p n } and continues with the next viable pair of levels minimizing ∆Ṽ ∆˜ . Inspection of all cases then reveals that the (first phase) of part II consists of k or k + 1 two-level rotations labelled by l = l min , . . . , k, where and for each of these rotations we choose Since ∆Ṽ ∆˜ does not depend on l and the rotations all commute (they pertain to different subspaces), the order of these operations within the first phase is irrelevant. However, not even all k (or k + 1) rotations may The probability weights pn of the initial state decrease with increasing energy. The initial average energy o (we use the dimensionless variables here for simplicity) is to be raised by ∆ to a value , that is closer to the energy level k = = 4 rather than = 5. (b) Part I: After rearranging the probability weights to place the largest weights closest to k, one obtains a distribution {pn} with an average energy of˜ I > . The numbers m = 0, . . . , 10 above the vertical lines at horizontal position n indicate that the corresponding weightpn corresponds to the value of the original weight pm. (c) Part II: Additional two-level rotations adjust the energy to the target . The first two of these rotations (corresponding to ϕ = 1 with j = 0 and l = 1, 2, see Sec. A.1.II) have angles π 2 and hence completely exchange the populations of the levels (m, n) = (3,5) and (2,6). The third rotation between levels 1 and 7 requires only a smaller angle 0 < θ < π 2 until reaching the target energy.
generally be enough to sufficiently adjust the average energy. Consequently, part II may consist of an arbitrary number of phases labelled by ϕ = 1, 2, . . ., where Let us now give a more compact description of part II. After part I, set {p} as the initial distribution, and further set˜ =˜ I ,Ṽ =Ṽ I , and ϕ = 1. Then perform the following steps: (i) Set j = j(ϕ), and l = l min (ϕ).
(iii) If˜ < , then ∆˜ II max > 0 and θ l is set to the value If˜ > , then ∆˜ II max < 0 and θ l is set to the value (A.14) Then continue with step (iv).
After the conclusion of part II, the target energy has been reached,˜ = , and the average squared deviation from hence becomes the energy variance. The second part of the protocol is illustrated in Fig. A.1 (c) and the variances resulting from the protocol for different temperatures and input energies are shown in Fig. 3 of the main text.

A.2 Wigner representation of squared number operator
We do this by using the formulas of Eqs. (33) and (35). To this end, we start by rewritingN 2 in terms of the local position and momentum operators aŝ We then insert term by term into Eq. (34) and calculate for functions f of the operatorsx andp, where we have used that After some algebra we then find the phase space representation of the operatorN 2 to be given by The second term on the right-hand side can be understood in the distributional sense. That is, defining the distribution γ[g(p)] via the function γ(p) given by Finally, we can compute the integral of W(x, p) with the second term on the right-hand side of Eq. (A.21) and find 1 16π dx dp W(x, p) dq dpq 2p2 e i(p−p)q where we have integrate over p using Eq. (A. 19), followed by an integral over the delta function δ(q − 2y), and finally made use of Eqs. (A.22) and (A.23). It is then easy to see that only the term 2h(x, 0) remains after taking the derivatives and evaluating at y = 0. Using the normalization of the wave function in Eq. (A.24) and arrive at Due to the linearity of the Wigner function in ρ, this computation extends from | ψ to arbitrary single-mode mixed states, and since the number operator of each mode is a local observable also to arbitrary N -mode states. We can hence rewrite Eq. (A.21) as

A.3 Extremal Gaussian precision and fluctuations
In this appendix, we give detailed proofs for the existence and uniqueness of the extremal values of the charging precision and fluctuations when restricting to single-mode Gaussian unitaries at fixed initial temperature and energy input. The corresponding results are presented and discussed in Sections 4.C and 4.E.

A.3.I Maximal variance
We begin with the maximally possible variance that single-mode Gaussian unitaries allow for, i.e., the worstcase scenario. Here, one is interested in determining the maximum of the function V + (r) from Eq. (54) over all r ≥ 0. For brevity, we will (again) use the notation ∆ = ∆E/ω and ν = coth βω 2 , as well as define V + := 4V + /ω 2 , such that the function that we wish to maximize can be written as To determine the extremal points, we calculate the first and second partial derivatives of V + w.r.t. r, i.e., The condition ∂V+ ∂r | r=rextr. = 0 at the extremal points r = r extr. yields Now note that the left-hand side is greater or equal than 1, while the function f (r) = e −2r cosh(4r) appearing on the right-hand side satisfies f (r = 0) = 1 and has a minimum at r = ln(3)/8 > 0 (as can be seen by setting ∂f (r)/∂r = 0), and f (r) r → ∞ −→ ∞. Consequently, Eq. (A.30) has two solutions r ± extr. , with r + extr. > ln(3)/8 > 0 and r − extr. < 0. Inserting (2 ∆ ν + 1) from Eq. (A.30) back into the second partial derivative gives which is negative for r + extr. (which is > ln(3)/8) and positive for r − extr. < 0. We thus have a local maximum at r = r + extr. , while r − extr. is a local minimum. Then, recall that we are interested in non-negative solutions (a negative squeezing parameter would reverse the roles of ξ 1 and ξ 2 in our treatment, see Section 4.C), and can thus discard r − extr. . Moreover, this eliminates all values of V + (r) for r < r − extr. , which can become larger than V + (r + extr. ). In other words, in the relevant range r ≥ 0, a (global) maximum of V + (r) can be found at r + extr. . Nevertheless, this is not the sought-after maximum for the variance, as we shall explain next.
If r + extr. > 0, implicitly determined by Eq. (A.30), were the correct solution, we could express the total input energy as < 0.
In other words, the (local) maximum of the function V + (r) at r = r + extr. is not physically realizable, since we must require ∆ d ≥ 0. Put simply, the maximum at r = r + extr. would require more energy to be invested into squeezing than is available overall, ∆ sq > ∆ . Since ∆ sq is strictly increasing with increasing squeezing parameter, we are hence looking for a solution for some r = r + < r + extr. that maximizes V + within the physically allowed range. Our previous analysis informs us that such a solution exists uniquely. The function V + (r) has one local minimum at a negative argument, and one local maximum for r > 0, and must hence be strictly increasing between r and r + extr. . The solution we are looking for is thus unique and found for the maximal value r = r + allowed by the global energy constraint, that is Expressing r + as a function of ∆ and ν then yields the result presented in Eq. (55), i.e., the worst precision for Gaussian single-mode unitaries is achieved when all energy is invested into single-mode squeezing.

A.3.II Minimal variance
Next we are interested in determining the optimal strategy using Gaussian single-mode unitaries. To this end, we similarly define V − := 4V − /ω 2 with V − as in Eq. (54), that is, we have to minimize The partial derivatives w.r.t. to r yield The extremal condition ∂V− ∂r | r=rextr. = 0 then yields (2 ∆ ν + 1) = e 2rextr. cosh(4r extr. ) , (A. 38) which has a unique solution r extr. = r − for r extr. ≥ 0 since the function e 2r cosh(4r) is greater or equal than 1 and is strictly increasing for r ≥ 0. Therefore, we have one and only one solution r − ≥ 0 for every value of (2 ∆ ν +1). Moreover, inserting into the second derivative gives Moreover, for the minimum we find that the energy input splits into squeezing and displacement according to such that, unlike the maximum at r + extr. discussed before, the desired minimum can be physically realized for all ∆ and ν.

A.3.III Maximal fluctuations
Let us now determine the maximal fluctuations that are possible during a charging process at fixed input energy via single-mode Gaussian unitaries. That is, we are interested in finding the maximum value of (∆W + (r)/ω) 2 from Eq. (67) over all r for fixed ∆ and ν. To simplify this task, we note that this is equivalent to the maximization problem for the function W + (r), given by where we have inserted for the partial derivatives of V + from Eqs. (A.28) and (A.29). For the purpose of solving the maximization problem, we introduce the notation χ := 2 ∆ ν + 1 ≥ 1 and λ := (ν − 1) 2 /ν 2 , with 0 ≤ λ ≤ 1 since ν ≥ 1, such that the extremal condition ∂W + /∂r r=r + = 0 at the extremal point r =r + reads λ 1 2 (1 − e −4r + ) + 1 2 e 2r + + e −6r + = χ. The maximization problem for W + (r) can thus be formulated as the question: Does there exist a u = u χ , with 0 < u χ ≤ 1, for every pair of λ and χ, such that f λ (u = u χ ) = χ? To answer this question, we first determine the extremal point u λ of f λ , i.e., such that which implies g(u λ ) := 1 2 3u λ − 1 u 3 λ = λ for u λ > 0. Since g(u λ ) is a continuous, strictly increasing function of u λ that can take the values g(u λ = 3 −1/4 ) = 0 and g(u λ = 1) = 1, there is exactly one u λ that satisfies g(u λ ) = λ for any λ between 0 and 1, suggesting that f λ (u) always has a unique local extremal point within the interval ]0, 1[. Inspection of the second partial derivative, i.e., then reveals that the local extremum is a local minimum. Moreover, since f λ (u = 1) = 1, this means that the minimal value is below one, f λ (u λ ) < 1. In contrast, we have f λ (u) u → 0 −→ ∞, suggesting that f λ (u) is strictly decreasing on the interval [0, u λ [ and can there take any value between f λ (u λ ) < 1 and ∞ (in particular, any value χ). We have thus found that there is a unique u χ < u λ for every λ and χ such that f λ (u χ ) = χ.
In other words, Eq. (A.46) has a unique solutionr + for every valid λ and χ. To check that this solution is a maximum, we calculate which is negative since u χ < u λ is below the minimum u λ of f λ (u). Consequently, the extremal value of W + is a maximum, which we have thus shown exists and is unique for any fixed λ and χ corresponding to fixed values of ∆E and T .

A.3.IV Minimal fluctuations
Similarly, we now wish to determine the minimal fluctuations that are possible during a single-mode Gaussian unitary charging process at fixed input energy, i.e., we want to minimize (∆W − (r)/ω) 2 from Eq. (67) over all r for fixed ∆ and ν. As in the previous section, we simplify the problem by considering the equivalent minimization of the function W − (r), given by To verify, that this condition can be met for all λ and χ, we again define a new variable v := e −2r with 0 ≤ v ≤ 1 for r ≥ 0, and define a family of functions h λ (v) via The minimization problem for W − (r) can thus be formulated as: Does there exist a v = v χ , with 0 < v χ ≤ 1, for every pair of λ and χ, such that h λ (v = v χ ) = χ?
To provide an answer, we start again by determining if h λ has any extremal points in the allowed range of v.
Finally, we check that we have indeed found a minimum of W − by evaluating the second partial derivative at r =r − (corresponding to v = v χ ), i.e., where we have substituted e −2r − = v χ . Inserting for χ = h λ (v χ ) from Eq. (A.56), and comparing with Eq. (A.57), we obtain which is nonnegative since h λ (v) is strictly decreasing for 0 ≤ v ≤ 1, and hence has a negative first derivative on this interval. We can therefore finally conclude that the extremal value of W − is a minimum that exists and is unique for any fixed λ and χ corresponding to fixed values of ∆E and T .