Geometry of work fluctuations versus efficiency in microscopic thermal machines

When engineering microscopic machines, increasing efficiency can often come at a price of reduced reliability due to the impact of stochastic fluctuations. Here we develop a general method for performing multi-objective optimisation of efficiency and work fluctuations in thermal machines operating close to equilibrium. Our method utilises techniques from thermodynamic geometry, whereby we match optimal solutions to protocols parameterised by their thermodynamic length. We characterise the optimal protocols for continuous-variable Gaussian machines, which form a crucial class in the study of thermodynamics of microscopic systems.

When engineering microscopic machines, increasing efficiency can often come at a price of reduced reliability due to the impact of stochastic fluctuations. Here we develop a general method for performing multi-objective optimisation of efficiency and work fluctuations in thermal machines operating close to equilibrium. Our method utilises techniques from thermodynamic geometry, whereby we match optimal solutions to protocols parameterised by their thermodynamic length. We characterise the optimal protocols for continuous-variable Gaussian machines, which form a crucial class in the study of thermodynamics of microscopic systems.
Designing optimal protocols for heat-to-work conversion below the nanoscale remains an ongoing challenge in the fields of stochastic and quantum thermodynamics [1][2][3]. For microscopic machines, efficiency and power alone are not the only figures of merit due to the additional influence of stochastic fluctuations. While a machine may extract work efficiently on average, it may be subject to large work fluctuations which hampers any reliability. Considerable effort is now being devoted to study this interplay between efficiency, power and reliability in small scale systems [4][5][6][7][8][9][10][11][12][13].
Dissipation can be minimised for general far-fromequilibrium processes using methods from optimal control theory [14][15][16][17][18][19]. However finding solutions with such approaches is often limited to simple systems with few degrees of freedom, which makes it difficult to identify general design principles for efficient thermal machines beyond specific models. On the other hand, a general method for optimising the efficiencies of machines operating close to equilibrium was recently proposed by Brandner and Saito [20]. This method, applicable to both classical and quantum periodic heat engines, relies on expressing the engine's entropy production in terms of a metric over the Riemann manifold of equilibrium states of the working system. One can maximise the efficiency for any given protocol by reparameterising it in terms of the so-called thermodynamic length [21][22][23][24][25][26][27][28][29][30][31], which provides a measure of distance between configurations in the equilibrium manifold. The benefit of this approach is its simplicity; optimisation is achieved by a straightforward computation of the thermodynamic metric tensor which depends only on the equilibrium and relaxation properties of the machine [20].
While thermodynamic length provides a systematic way of determining efficient protocols, it remains to be seen how increasing efficiency impacts the work fluctuations. For systems connected to a single fixed-temperature reservoir, initial investigations have explored the simultaneous optimisation of the average dissipated work required to drive a system from one state to another and the associated minimal work fluctuations [9,32]. As this is a multi-objective optimisation problem, one must consider the boundary of allowed protocols where dissipation cannot be reduced any further without suffering an increase in fluctuations, or vice versa. This boundary, known as a Pareto front [33], only converges to a single point in regimes where the fluctuationdissipation relation holds true [34][35][36], in which case there exists a unique protocol with minimal average dissipated work and work variance. At present, Pareto optimisation has not been analysed in the context of periodic heat engines operating between different temperatures. In this situation the figures of merit are instead thermodynamic efficiency versus the resulting work fluctuations, whose optimal protocols are not expected to coincide.
In this paper, we outline a general method for finding Pareto-optimal protocols interpolating between maximum efficiency and minimal work fluctuations for engines operating close to equilibrium. Remarkably, we show that such protocols can be found by constructing a new form of thermodynamic metric tensor and corresponding length. By parameterising any given protocol in terms of this generalised thermodynamic length, one may identify regimes of optimal efficiency and reliability in a straightforward manner. We illustrate our approach for quantum heat engines operating along discrete step-equilibration cycles [37][38][39][40], though our results also apply to classical heat engines and open quantum systems undergoing Lindblad dynamics. As a core application, we derive analytic expressions for thermodynamic length in general continuous variable Gaussian quantum systems [41]. Such systems form a major platform for studying thermodynamic processes in the microscopic regime, and our approach can be used to optimise any Gaussian thermal machine by using the natural tools in the Gaussian formalism that operate on the steady-state covariance matrix. As an example we determine Pareto fronts of optimal efficiency and reliability for a system of coupled harmonic oscillators driven by periodic changes in temperature, frequency and coupling strength.
Let's begin by considering a quantum system weakly coupled to a thermal environment at inverse temperature β = 1/k B T . The system is subject to external control via a set of d mechanical parameters λ := (λ 1 , λ 2 , ...λ d ) and we denote its Hamiltonian by H( λ). In addition, we allow for external modulation of the environmental temperature. Collectively these set of variables define a cycle via a closed curve γ : t → Λ t in the parameter space containing the vectors and we label Λ 0 = β and Λ j = λ j for j ≥ 1. The parameter space defines a manifold of equilibrium states, defined by π( Λ) := exp (−βH( λ))/Tr exp (−βH( λ)) . Furthermore, we will introduce the following conjugate forces, During the cycle the inverse temperature undergoes a variation between a maximum and minimum, β h ≤ β(t) ≤ β c , which we will express in the form β(t) with α(t) ∈ [0, 1] a periodic function. For convenience we take t ∈ [0, 1] to be a dimensionless parameter, and denote a discretised set of N points along the curve evaluated at times t n = (n − 1)/(N − 1) for n ∈ [1, N ]. We have in mind thermodynamic cycles where each step can be approximated by a fast quench in the mechanical parameters λ tn → λ tn+1 , followed by relaxation to a new temperature β tn → β tn+1 . This means that at the beginning of each n'th step the system is in a thermal state π( Λ tn ), which is left unchanged during the quench step while work is performed. The state after the quench then relaxes to a new equilibrium state π( Λ tn+1 ) with no work done during this step.
A central quantity of interest is the average irreversible entropy production along the cycle γ, which can be expressed as where S(ρ||ρ ′ ) = Tr ρ ln (ρ) − Tr ρ ln (ρ ′ ) ≥ 0 the quantum relative entropy, and we identify the work done and supplied heat The second equality in (3) follows from using the periodicity Λ t1 = Λ tN (see Appendix B). This formula motivates a definition of efficiency for processes with a positive work output, W ≤ 0, given by the ratio [42,43] η where η C = 1 − β h /β c denotes the Carnot efficiency. Here we see consistency with Carnot's theorem, which follows as a consequence of the second law S irr ≥ 0.
In addition to efficiency, we will also be concerned with the amount of work fluctuations generated along the cycle. For quantum systems, stochastic work is determined from projective measurements of the system energy at the beginning and end of each unitary step [44]. By summing up the changes in energy across each step and computing the variance from the resulting work probability distribution, one can show [32,40] Var(W ) := At this stage, we restrict our attention to cycles composed of a large number of steps N 2 ≫ 1, which defines a regime that is close to quasi-static [37]. In this case we may replace the summation with an integral over the continuous path γ and obtain where summation is carried out over repeated indices, restricted to the mechanical variables only. Here m jk is a positive semi-definite and symmetric tensor and we set m j0 = 0 ∀j for later convenience. We also denote the shifted operators δX j ( Λ) = X j ( Λ) − Tr X j ( Λ)π( Λ) . Turning to the efficiency, we consider the fraction of the efficiency below the Carnot value, defined as δη := 1 − η ηC . Assuming a large number of steps and applying a Taylor expansion in 1/N to (3) gives where we now have a different metric tensor: and we have introduced the unitary channel U ν, λ [.] = e iνH( λ) [.]e −iνH( λ) . We have further defined the adiabatic work done [20]: which is a geometric quantity independent of the parameterisation and assumed negative W ≤ 0 to ensure a useful work extraction cycle. We provide a proof of (9) and (11) (16), corresponding to work fluctuations optimisation at ǫ = 1 (the blue curve) and efficiency optimisation with ǫ = 0 (the red curve). The parameters are set to ω0 = 2, Tc = 0.25ω0, T h = Tc + ∆T with ∆T = 2.5ω0, κ = 0.2ω0. Right-Work fluctuations vs. oscillator frequency ω0, for the linear protocol (red) and the optimal protocol (blue). We denote ∆W = Var(W ) as the standard deviation and ∆W * = kBTcL1/ √ N the optimal fluctuations given by (17). Both quantities are plotted relative to the adiabatic work extracted |W|, which is independent of parameterisation. The gray area is not accessible through any protocol. Here we set the parameters to N = 50, Tc = 0.5, T h = Tc + ∆T with ∆T = 5Tc, and κ = 0.2Tc.
Fisher information metric [45,46], which is a quantum analogue of the classical Fisher-Rao metric. At this stage we observe that the elements of this tensor typically differ from m jk due to possible non-commutativity between the conjugate forces, ie. if [X j , X k ] = 0. In Appendix A we highlight the different information-geometric interpretations of these two metrics.
To establish our optimisation problem, let us introduce the dimensionless multi-objective function where for convenience we have defined the work fluctuations in units of the cold temperature, Var(W ) = β 2 c Var(W ). The question we now address is: how long should the system spend at each point along the protocol γ in order to maximise efficiency while minimising fluctuations? This amounts to finding the best choice of parameterisa- to be determined so as to minimise the scalarized objective (13). Optimal parameterisations I * ǫ ≤ I ǫ lie along sections of the Pareto fronts [33]; these points form the boundary of protocols where it is not possible to increase efficiency (ie. reduce δη) without increasing work fluctuations, or conversely, reduce fluctuations without reducing efficiency. By combining (8) with (10) and applying the Cauchy-Schwarz inequality we arrive a geometric expression for the minimised objective function: and we define This follows from the fact that M ǫ jk ( Λ) gives another metric tensor, since it is formed by a positive-weighted linear combination of two metric tensors for ǫ ∈ [0, 1]. The function L ǫ may be interpreted as a generalised form of thermodynamic length [27], whose dependence on ǫ encodes information about the Pareto optimal solution. Crucially, these Pareto optimal solutions are determined by parameterising the protocol γ in terms of the modified thermodynamic length via the speed function t → φ ǫ t [47], which is obtained from the implicit equation: This means that for any given protocol, optimisation is achieved by changing the speed at which the curve is traversed by choosing a new parameterisation There are two limiting cases of the bound. For ǫ = 1 we obtain a geometric lower bound on the achievable work fluctuations: For ǫ = 0, we obtain a maximum upper bound on efficiency η ≤ η C (1−L 2 0 /N ). An analogous efficiency bound was previously obtained in [20] for continuous Lindblad dynamics.
The above construction gives a general recipe for finding Pareto optimal protocols for arbitrary quantum or classical systems, valid in regimes where the number of steps between equilibrium states is large. A particular class of systems that are frequently used to describe many relevant physical systems in thermodynamics, such as ion trap heat engines [48], are composed of Gaussian quantum states [41,49]. The Hamiltonian of a D-mode Gaussian system is quadratic in quadrature operators, namely with R = (x 1 , p 1 , . . . , x D , p D ) T being the quadrature vector. Here, the D × D dimensional symmetric matrix G λ contains all the quadratic couplings. For this general class we provide an analytic formula for the Pareto optimal solutions, which are given by (16) via computing the metric tensor (15). By defining X j := ∂ Λj G λ if j ≥ 1 and X 0 := β −1 G λ we can express the metric (15) as follows (see Appendix C): where a j0 = a 0j = 0 ∀j and a jk = 1 ∀j, k > 0, and we defineX with Ω being the symplectic form with Ω nm = i[R n , R m ], and σ Λ representing the steady state covariance matrix with elements [σ Λ ] nm = Tr π Λ {R n , R m } /2 − Tr π Λ R n Tr π Λ R m . Furthermore, the adiabatic work is found using W = 1/2 γ dλ j tr X j σ Λ . Notice that the 'tr' operation acts as a trace on the matrix space associated to the Gaussian covariance matrices, which should be distinguished from the trace 'Tr' which acts on the Hilbert space for density operators.
We have now derived the general form for the thermodynamic metric tensor for Gaussian heat engines. So long as one can compute this metric tensor, the speed function φ ǫ t can be approximately determined from (16) via point-wise inversion followed by numerical interpolation. We illustrate our method for the example of a pair of coupled harmonic oscillators with R = [x 1 p 1 x 2 p 2 ] T and Hamiltonian coefficient matrix were we chose equal frequencies ω 1 = ω 2 = ω and denote κ the coupling strength between the oscillators. As for the driving protocol, we consider control over the bath temperature alongside the joint frequency and coupling, ie.
The matrices X ω = ∂ ω G λ , X κ = ∂ κ G λ , and X 0 = β −1 G λ are easily found from Eq. (22). By substituting in Eqs. (20) and (21) one finds the correspondingX j and X j . Finally by plugging these into Eq. (19) we find the metric [50]. We choose a harmonic protocol path β(t) = β c + (β h − β c )sin 2 (πt), ω(t) = ω 0 1 + sin 2 (πt + π 4 ) , and κ(t) = κ 0 1 + sin 2 (πt + π 4 ) , with the parameters κ 0 , ω 0 , β c and β h being fixed during the cycle. In Figure 1 we compare the work fluctuations ∆W = Var(W ) for a linear parameterisation as a function of the optimal amount ∆W * = k B T c L 1 / √ N given by (17), both as a function of oscillator frequency ω 0 and expressed in units of the adiabatic work. We also plot the rate of change in the speed function φ ǫ t versus time for ǫ = 1, giving minimal fluctuations, compared with ǫ = 0 that gives maximum efficiency. One can clearly see that distinct protocols must be chosen in order to achieve either optimal efficiency or fluctuations. In Figure 2 we plot the Pareto fronts for the model for different choices of cold temperature T c and coupling constant κ. The points on each curve give the corresponding values of δη and ∆W for the optimal protocol Λ ′ t = Λ φ ǫ t determined from (16) for every ǫ ∈ [0, 1]. These curves form the boundary of achievable fluctuations and efficiency for the chosen protocol γ, with points below the curves inaccessible. In this case the Pareto fronts are strictly convex, and hence the entire front is determined by the minima of the scalarised objective function (13) [51].
To summarise, we have constructed a general method for performing multi-objective optimisation of efficiency and work fluctuations in microscopic heat engines operating close to equilibrium. This method relies on determining a thermodynamic metric tensor that encodes information about both the efficiency and fluctuations simultaneously.
Our formalism opens an avenue to further applications and generalisations. For example, in Appendix D, we show that an analogous metric providing Pareto optimal solutions for work fluctuations and efficiency can be derived for engine cycles described using Lindblad dynamics, and we derive the corresponding expressions for Gaussian Lindbladians. Future investigations could focus on extending our approach to strongly-coupled quantum heat engines [52], or regimes far away from equilibrium [19,53]. The Gaussian metric tensors we have derived can also be used to study other aspects of thermodynamic geometry, such as computing scalar curvature [26] and geodesics along the manifold of thermal Gaussian states [28].

A. Thermodynamic metrics and Fisher information
In this Appendix we provide further background concerning the information-geometric interpretations of the metric tensors (9) and (11) for both classical and quantum mechanical systems.
Consider first a d-dimensional statistical manifold with coordinates Λ, with each coordinate defining a normalised probability distribution p n ( Λ) with discrete outcome space x n ∈ (x 1 , x 2 , ...). We denote the expectation value of an observable A with respect to p n ( Λ) as A = n p n ( Λ)A n . A particular divergence measure between two distributions on the manifold is the Kullback-Liebler divergence: which quantifies how distinguishable two distributions are from each other. If one considers two close points Λ and Λ + d Λ, the divergence between the distributions becomes where is known as the Fisher-Rao matrix. We may interpret this as a Riemann metric tensor on the manifold, since it is positive semi-definite, symmetric and smooth in Λ. This implies a notion of distance between points in the manifold, with squared line element ds 2 = F jk ( Λ)dΛ j dΛ k .
For quantum systems, one may instead consider a manifold of normalised density matrices ρ( Λ). The quantum relative entropy replaces (A2), defined as Considering again close points Λ and Λ + d Λ, we have where is known as the Kubo-Mori Fisher information matrix [45,46]. This provides a measure of distance between neighbouring density matrices on the manifold with squared line element ds 2 = F jk ( Λ)dΛ j dΛ k . We note however, that other information metrics may be obtained by replacing the relative entropy by a different choice of divergence. For example, consider instead the relative entropy variance [55] V Expanding this between neighbouring states gives whereF jk ( Λ) := defines another metric. This choice was first introduced in [56] as an alternative to the Kubo-Mori metric. Note that for quasi-classical states with a spectral decomposition of the form ρ( Λ) = n p n ( Λ)|n n|, with eigenstates {|n } that are independent of coordinates Λ, both metrics (A6) and (A9) become proportional to the Fisher-Rao metric (A3). The difference between these two metrics highlights the role of quantum coherence. If we restrict our attention to the manifold of thermal states ρ( Λ) = π( Λ) := exp (−βH( λ))/Tr exp (−βH( λ)) and introduce conjugate forces {X i } defined in (2), we see this difference between the metrics occurs when at least one pair of forces are non-commuting, [X i , X j ] = 0. By comparing the different metric expressions used in the main text, we find that the work fluctuations (9) are determined by components of the metric (A9), while the efficiency (11) is determined by the Kubo-Mori metric (A6). (8) and (10) In this section we provide details of the derivations of the metric expression for work fluctuations and efficiency. Let us begin by considering the protocol γ : t → Λ t with t ∈ [0, 1], evaluated at N discrete points t n = (n − 1)/(N − 1). We introduce the operator ∆Φ n /N = Φ( Λ tn+1 ) − Φ( Λ tn ), where Φ( Λ) = βH( λ) + lnZ( Λ) is the non-equilibrium potential with lnZ( Λ) = Tr exp (−βH( λ)) the partition function. Similarly we define ∆H n /N = H( λ tn+1 ) − H( λ tn ). It is then straightforward to see that the work fluctuations (7) to leading order in 1/N are

C. Gaussian formalism
In this section we provide a derivation of the analytic expressions for the metrics for Gaussian states, namely (19). Let us start by finding the elements X i in Eq. (2) for the Gaussian scenario. By defining the matrix Σ whose elements are the second order symmetric quadratures Σ ij = 1/2{R i , R j } and using the fact that the matrix G λ and its derivatives are symmetric, we have Tr π( Λ)tr (X j Σ) Where we defined X j := ∂ Λj G λ . Notice that "Tr . " is different from "tr ()", which represents the expectation value. As for the temperature element X 0 we have where we define X 0 := β −1 G λ .
In the Gaussian formalism we know how the displacement vector and the covariance matrix-and hence all higher order moments-evolve under unitary transformations. In particular, under a unitary that is generated by the Hamiltonian where we use the BCH lemma, and define S ν G λ = e −iνΩG λ . With these at hand we can work on the metric (11) where we used the Wick's theorem in order to expand the fourth order correlations in terms of second moments. Recall that we also have the first moments, and hence all odd moments, vanish. In the same manner, we can find the work fluctuations metric: which completes the proof of (19). Finally, If we substitute Eq. (C1) in (12) of the main text, we find the adiabatic work in the Gaussian formalism In evaluation of (C10) recall that the temperature element is not included in the integration, i.e., j ≥ 1.

D. Thermodynamic geometry for open quantum systems
In this section, we demonstrate how our formalism in the main text can be extended beyond step-equilibration processes to continuous Markovian processes. Rather than modelling the system evolution by a sequence of Hamiltonian quenches followed by relaxation, we instead consider a weakly coupled open system ρ t whose evolution over time interval t ∈ [0, τ ] is given by a time-dependent Lindbladian:ρ This Lindbladian depends on both temperature and the mechanical variables of the corresponding Hamiltonian H( λ), collectively labelled by Λ = {β, λ} as before. We assume the evolution obeys quantum detailed balance [57] with a unique thermal fixed point for each parameter Λ The process is cyclical and described by a closed curve γ : t → Λ t . In this case, the average work done and irreversible entropy production As we saw with (B4), the corresponding efficiency of the process can be related to the ratio between entropy production and work [43]: We now restrict to a regime where the system remains close to equilibrium at all times. This occurs whenever the characteristic timescale t eq of the system is always small compared to the total duration, namely (t eq /τ ) 2 ≪ 1 [58]. This slow driving regime can be thought of as an analogue of the large step approximation used in the main text. By expanding up to linear order in t eq /τ , it can be shown (see [20,59]) that the fraction of efficiency below Carnot is where W is the adiabatic work as defined in (12), and we introduce a new metric tensor whereg with L † the adjoint.
To quantify the work fluctuations Var(W ), one needs to unravel the system evolution in terms of quantum jump trajectories [60][61][62]. This amounts to monitoring the energy exchanges with the environment, with jump statistics determined by Born's rule. The precise formalism for computing Var(W ) can be found in [59]. There it is shown that under slow driving the work fluctuations can be expressed as: The corresponding metric is given by andm j0 = 0 ∀j > 0. With these two metrics, one may minimise the objective function (13) by constructing the analogous metric (15). Going further, we now derive a set of closed expressions for evaluation of the thermodynamical quantities in the Gaussian formalism. The main difference with the non-dissipative case-that was presented in the main text and proved in the previous section-is the presence of a Gaussian Lindbladian master equation-which is determined by the term e νL † Λ . Since the dissipative dynamics is Gaussian too, it can always be characterized efficiently. Generally speaking for an arbitrary observable O we have with H λ = 1 2 R T G λ R being the Hamiltonian, and the Lindbladian operators L k, Λ = c T k, Λ R are linear in quadratures. Here, c k, Λ are 2D dimensional complex vectors.
Like any other Gaussian quantuman channel one only needs to identify how they transform the first and second order moments. In our case, we only focus on vanishing first order moments. Then the dissipative Gaussian channel can be characterized in the Heisenberg picture as follows where the matrices F ν, Λ and Y ν, Λ can be found from the basic elements of the Gaussian Lindbladian master equation i.e, G λ and c k, Λ . Specifically, we have F ν, Λ = e νA Λ and Y ν, Λ = ν 0 dν ′ e ν ′ A Λ D Λ e ν ′ A T Λ with A Λ = −iΩ(G λ − Im(C Λ C † Λ )), and D Λ = Ω Re(C Λ C † Λ )Ω, where we define C := (c T 1 ; c T 2 ; . . . c T m ) T -see e.g., Section 5 and Appendix C of [49] for derivation. The application of the channel to σ Λ should be understood through its application on the identity operator, because in fact by [σ Λ ] jk we mean [σ Λ ] jk I. Since, the map is unital, it leaves σ Λ unchanged. Thus, when applied to Σ − σ Λ we have where we use the fact that σ Λ is the fixed point of the dissipative dynamics i.e., σ Λ = F ν, Λ σ Λ F T ν, Λ + Y ν, Λ for ∀ν. Putting everything together we can find the metrics. (D10) and (D7). Firstly, we havẽ m jk ( Λ) = 1 2 ∞ 0 dν Tr π Λ e νL † Λ δX j ( Λ) , δX k ( Λ) where we used the Wick's theorem in order to expand the fourth order correlations in terms of second moments. We also extend the definition of X j -from Eq. (21)-to the open dynamic scenario Moreover, the matrixg jk ( Λ) can be found in a similar manner: whereX j is still given by Eq. (20) whereas X k is given by (D17) above.