Thermodynamics of precision in quantum non equilibrium steady states

Autonomous engines operating at the nano-scale can be prone to deleterious fluctuations in the heat and particle currents which increase, for fixed power output, the more reversible the operation regime is. This fundamental trade-off between current fluctuations and entropy production forms the basis of the recently formulated thermodynamic uncertainty relations (TURs). However, these relations have so far only been derived for classical Markovian systems and can be violated in the quantum regime. In this paper we show that the geometry of quantum non-equilibrium steady-states alone, already directly implies the existence of a TUR, but with a looser bound. The geometrical nature of this result makes it extremely general, establishing a fundamental limit for the thermodynamics of precision. Our proof is based on the McLennan-Zubarev ensemble, which provides an exact description of non-equilibrium steady-states. We first prove that the entropy production of this ensemble can be expressed as a quantum relative entropy. The TURs are then shown to be a direct consequence of the quantum Cramer-Rao bound, a fundamental result from parameter estimation theory. By combining techniques from many-body physics and information sciences, our approach also helps to shed light on the delicate relationship between quantum effects and current fluctuations in autonomous machines.


INTRODUCTION
Autonomous machines, whether classical or quantum, generically operate in non-equilibrium conditions. They harvest work by consuming resources such as heat or fuel and, in order to maintain functionality, dissipate into the environment [1]. As such they operate in a regime known as non-equilibrium steady state (NESS), characterized by a nonzero entropy production rate [2] and by an ability to maintain non-zero average currents across the system. Accurately describing the physical properties of a NESS is central to the development of mesoscopics [3,4] as well as fundamental to our understanding of nano-scale autonomous machines such as molecular electronics [5,6], nano-junction thermoelectrics [7], single electron circuits [8], quantum dots [9] and even ultra-cold atoms [10].
An important advance in this respect has been made recently with the discovery that classical Markovian systems obey the so-called thermodynamic uncertainty relations (TURs) [11,12]. The basic idea is that in meso-and microscopic systems, fluctuations of the currents around their mean values become significant. The TUR provides a bound on these fluctuations by relating them to the NESS entropy production rate according to (taking k B = 1 and = 1 throughout) where Ĵ α is the average current (of particles, charge or heat), ∆Ĵ α ≡ lim T→∞ T Ĵ 2 α − Ĵ α 2 denotes its normalized variance and σ is the average entropy production rate in the NESS. Since the original inception of this result, there has been a flurry of activity aimed at further exploring their con-sequences in various settings [13][14][15][16][17][18]. Not only is the TUR expected to have implications for the functioning of biological clocks [19] and control techniques [20], it was demonstrated recently that it has significant consequences for the operation of autonomous machines. For instance, it follows from Eq. (1) that the fluctuations in the output power of an engine operating between two reservoirs at temperatures T C and T H > T C is bounded by [21] P ≤ where P denotes the average power, ∆ P its normalized fluctuations, η the engine's efficiency and η C = 1 − T C /T H the corresponding Carnot efficiency. This result implies that, in order to have an autonomous steady-state machine which operates at finite power as η → η C , one must incur fluctuations that diverge at least as ∼ (η c − η) −1 . Recently, it has been shown [22,23] that the classical TUR Eq. (1) can be violated in the quantum regime. While the precise mechanisms responsible for these violations are still not fully understood, this opens up an interesting perspective, as it would in principle allow one to use quantum effects to reduce the deleterious current fluctuations without compromising the engine's efficiency and output power [24][25][26]. Moreover, these violations also naturally lead one to ask, to what extent, are TURs really a universal feature of non-equilibrium steady-states.
In this work we show that the geometry of quantum NESS, by itself, already implies the existence of a TUR of the form where ∆ is the normalized covariance matrix between different steady-state currents. By restricting to the diagonal ele-ments, one immediately in particular obtains which, compared to the classical result in Eq. (1), shows that our bound involving the variance of currents is four times looser. This result shows the universal nature of TURs, generalizing the concept of thermodynamics of precision to a much larger class of physical processes. The key to achieve this result is, rather than adopting a dynamical approach, to exploit the generalization of the idea of Gibbs distributions to the set of NESSs known as the non-equilibrium statistical operator approach, or the McLennan-Zubarev form [27][28][29]. This description allows us to write the average entropy production in the NESS as a quantum relative entropy and investigate the geometrical structure of the manifold defined by this family of states. Making use of concepts borrowed from the geometry of quantum states and quantum estimation theory, this allows us to derive the geometric bound Eq. (3). We finally illustrate the implications of our findings to the output power of an autonomous mesocopic heat engines, in the same spirit as Ref. [20]. This leads to the following new general independent upper bounds where B P S is given by Eq. (2), ∆P ,ĴH denotes the normalized correlation between the power and the heat current from the hot reservoir, while ∆Ĵ H the variance of the latter. On the one hand, one can deduce from the first inequality in Eq. (5) that the maximum reachable power for a nanoscale steady-state engine is four times larger than any classical Markovian counterpart. Furthermore, our generalized quantum TUR allows us to see and quantify how much the fluctuations of the incoming heat current from the hot reservoir affects the achievable output power. A systematic study of these two bounds in several physical platforms and models will be pursued in a forthcoming work.

NON EQUILIBRIUM STEADY-STATE STATISTICAL OPERATOR
We consider here the typical NESS scenario depicted in Fig. 1(a). The main idea of the statistical operator approach is to use a density matrix ensemble description for the NESS, taking into account the additional conserved quantities which are responsible for the currents. This is based on propagating the entire composite system from the infinite past to the present using a generalized Gibbs ensemble first derived by McLennan [28,30] and Zubarev [27,31,32], and is hence known as the McLennan-Zubarev form. An alternative form for the NESS statistical operator was derived in the 90's by Hershfield [29] using a scattering-theory based approach. Not only are the approaches known to be equivalent [33][34][35] but they can also can be obtained from a max-entropy approach [36,37] in analogy with the well known case for usual equilibrium ensembles [38,39].  Figure 1. Diagrammatic illustration of the scenario considered in this paper. A central system (C) is coupled to two infinite thermal reservoirs prepared at temperatures TL,R and chemical potentials µL,R. After a sufficiently long time, this system will tend to a global nonequilibrium steady-state (NESS) (a) characterized by the existence of a finite current of particles (ĴQ)and energy (ĴE) through the central system. When the two biases in temperature and chemical potential is set to zero, i.e TL = TR = T and µL = µR = µ, the asymptotic state becomes an equilibrium state known as local equilibrium state (LES) (b). The two states are closely related to each other as shown in Eq. (12) For clarity, we will motivate this approach from the perspective of a central system connected to two fermionic baths (for a generalization of the NESS to the case of multiple baths, see e.g. Ref. [40]; for a derivation from a purely classical perspective see Ref. [28]). In the following we will consider a schematic model where a central conductor C, which itself could be a complex many-body system, is attached to two semi-infinite fermionic leads, L and R, acting both as thermal and particle reservoirs such that a global steady state is reached [see Fig. 1(a)]. The overall system is therefore described by the total HamiltonianĤ =Ĥ 0 +Ĥ int , witĥ H 0 =Ĥ C +Ĥ L +Ĥ R representing the sum of all the autonomous terms for each part, andĤ int =V LC +V RC incorporating all the couplings between the three parts.
At the infinite past, t 0 = −∞, the three components are considered decoupled, with the two baths being at equilibrium characterized by two different inverse temperatures β L and β R and two chemical potentials µ L and µ R . The total system is therefore taken to be at the stateρ(t 0 ) =ρ L ⊗ρ C ⊗ρ R , wherê in which Z α = Tr a [e ··· ] andN a is the total particle number operator for reservoir a. We are interested in the steady-state and so any observable of interest in this limit will be independent from the initial state of the central systemρ C . Immediately after the initial time t 0 the coupling between the central system and the two leads is switched-on adiabatically according toĤ ǫ (t) =Ĥ 0 + e −ǫ|t|Ĥ int , with ǫ being an arbitrary small positive constant [41], soĤ ǫ (t 0 = −∞) =Ĥ 0 ,Ĥ ǫ (0) =Ĥ. Making use of the interaction picture evolution operatorÛ I,ǫ (0, −∞) generated by this Hamiltonian, the NESS statistical operatorρ ness ≡ lim ǫ→0 +Û I,ǫ (0, −∞)ρ 0Û † I,ǫ (0, −∞) is given by the exact and intuitive form (see SI and Refs. [27,31,40,42]) where Z ness is the NESS partition function, Crucially, Eq. (7) includes the entropy production operatorΣ which is defined aŝ and for X = E, Q, we introduce the time-integrated operatorŝ withX H (t) = e iĤtX e −iĤt . The last equality in Eq. (10) is a direct consequence of Abel's theorem [42,43].Q + andÊ + therefore represent the operators connected to time-integrated particle and energy exchanged between the L and R leads. It is important to emphasize that in this framework, the NESS state Eq. (7) refers to the global state of the composite LCR system and not just the reduced state of C. The structure of Eq. (7) implies that the NESS is a generalized Gibbs ensemble where in addition to the usual conserved quantitiesĤ andN , there are the additional conserved quan-titiesÊ + andQ + , whose Lagrange multipliers are the well known thermodynamic affinities δ β and δ βµ that drive the energy and particle currents respectively [2]. Another observation is that Eq. (7) also illustrates that the cumulants of the steady-state entropy production, here connected to the current fluctuations, are generated by taking derivatives of the corresponding Massieu potential ψ ness ≡ − ln Z ness [44]. This is in direct analogy to the case for equilibrium ensembles where they are connected to the equilibrium fluctuations of the energy. In the SI we explicitly show that the expectation value of the entropy production operatorΣ for the NESS, calculated through the Massieu potential, recovers the familiar result for the average entropy production rate where Ĵ Q,E = Tr Ĵ Q,Eρness are the usual asymptotic steady-state particle and energy currents [34,45], with the current operator ofX given asĴ X (t) ≡ dXH (t) dt for X = E, Q.

ENTROPY PRODUCTION AS A RELATIVE ENTROPY
Armed with the NESS statistical operator we are now ready to introduce our first main result. Using the NESS ensemble Eq. (7) we may write the average entropy production as a quantum relative entropy [46], which plays a central role in non-equilibrium quantum thermodynamics [47][48][49][50][51][52]. In the formalism presented here, the entropy production operator Σ defined in Eq. (8) represents a conserved quantity in the NESS, i.e. it commutes with the total HamiltonianĤ. Consequently the exponential in Eq. (7) can be factorized leading to the following insightful re-expression whereρ les represents the local equilibrium state (LES) of the total system, i.e. the equilibrium condition that would be reached if both leads had an inverse temperatureβ and chemical potentialμ [See Fig. 1 where the second term on the r.h.s. corresponds to the difference in the Massieu potentials of the LES and of the NESS We stress that the expectation value of the even powers of the entropy production operator in the last equality are calculated over the local equilibrium state (see Supplementary Information -SI -for derivation and details). One may physically interpret D (ρ ness ||ρ les ) as the available work that could be extracted by letting the leads equilibrate to β and µ [53].
Two relevant considerations can be made. First, all the expectation values Σ 2n les are positive quantities, therefore leading to the conclusion that ∆ψ ≥ 0. Combining this with the fact that by Klein's inequality the relative entropy is snon-negative, allows us to obtain the non-negativity of the steady-state mean entropy production Σ . Along the same lines, upon introducing the entropy production rate σ from Eq. (11) into Eq. (13), one obtains an analogous reformulation of this quantity in terms of a relative entropy and also prove its positivity. Second, in the absence of a temperature gradient, the first term in the series, i.e. Σ 2 les , is proportional to the Johnson-Nyquist noise [42], which therefore becomes the leading term whenever an expansion over the affinity β(µ L − µ R ) is performed.

THERMODYNAMICS OF PRECISION
Using our derived expression for the average entropy production defined by Eq. (13) we are now in a position to exploit the mathematical properties of the relative entropy to derive our TUR. We start by rewriting the NESS operator in Eq. (7) in the generic form where Einstein's summation notation has been adopted, where X = (H, −N, E + , −Q + ) T and where the vector of affinities (or generalized forces) λ = β, βµ, δ β , δ βµ T physically represents the set of experimentally controllable parameters of the system, thus defining the manifold of thermodynamically accessible states. We will henceforth refer to this as the manifold of steady-states (SSM), see Fig. 2. It is immediate to see that equilibrium states of the formρ les will belong to this manifold, as they correspond to vectors λ having the last two affinities δ β and δ βµ equal to zero. The local curvature of the manifold is then given by the Fisher information, which quantifies the sensitivity of the system to small variations of the control parameter λ. Explicitly where the Fisher information I is in this multidimensional case a matrix with elements with {ρ k } denoting the set of populations, i.e. the projections of the density matrix on the the energy eigenbasis ofĤ. The SSM can be shown to be a Riemannian manifold over the set of parameters λ, whose metric (called Fubini-Study [54]) induces a notion of statistical distance between two generic statesρ 1,2 . In fact looking at the SSM in this geometrical form can be used to generalised the concept of thermodynamic length, introduced by Crooks for the class of equilibrium ensembles [55], to the class of NESS (see also extension to nonunitary dynamics [56]).
In order to shed light on the geometry of the thermodynamics of precision we invoke the Cramer-Rao bound [57] which puts a fundamental lower bound on the precision of estimation of a parameter λ, or a function g(λ) of that parameter, labeling a statistical ensemble. Concretely, the latter reads where K is the Jacobian matrix of transformation with elements and where Cov denotes the covariance matrix with elements Equation (18) can be interpreted as stating the positive semidefiniteness of the matrix Cov − KI −1 K T . Given this, the NESS defined in Eq. (7) is diagonal in the total energy eigenbasis so we can use the classical version of the bound. Motivated by what happens in actual experimental platforms, rather than estimating the vector of affinities λ we choose to estimate the average steady-state currents Ĵ α λ ≡ Tr Ĵ αρ (λ) , with α potentially labeling particle currents (α = Q), energy (α = E), heat (α = H), left/right (α = L, R) and/or even difference between these two which defines the work current (α = W ). Let us start from a state T being a small increment in the inverse temperature and chemical potential imbalances. These two states represent respectivelyρ les andρ ness . By combining Eq. (16), Eq. (18) and Eq. (13), one can finally arrive to the following inequality (see SI for details), where, following the standard procedure [58], we have defined the so-called normalized covariance matrix with elements and where we used the entropy production rate defined in Eq. (11) as σ = lim T→∞ T −1 Σ . We remind the reader that the thermodynamic limit lim T→∞ has to be performed at the very end of the calculation of interest, and therefore the ratio in Eq. (21) is well-defined and finite. From Eq. (21) it immediately follows in particular that which represents a TUR of the same form as Eq. (1) but with a constant which is four times looser, dictated by the geometry of quantum steady-states. It is however worth remarking that our result Eq. (21) generalizes the classical TUR in that it involves in general the full covariance matrix ∆.

IMPLICATIONS FOR MESO-AND NANOSCOPIC HEAT ENGINES
We will now discuss the consequences that our new bound Eq. (21) on precision has on autonomous quantum steadystate machines operating at small biases in temperature and chemical potentials. Let us assume, with reference to the color scheme used in Fig. 1 (a) and without loss of generality, that T L > T R and that the thermal gradient is exploited to drive the current against the chemical potential difference µ L < µ R (Seebeck effect) [45,[59][60][61].
In their seminal work, Pietzonka and Seifert (PS) showed that, for steady-state engines described by classical Markovian stochastic processes, the application of TUR to the work current, i.e. the power P ≡ Ĵ w = Ĵ L − Ĵ R (witĥ J a =Ĵ E,a − µ aĴQ,a , a = L, R, being the heat current from the a reservoir) leads straightforwardly to the upper bound Eq. (2), with η = P / Ĵ L being the efficiency of the engine and η C = 1 − T R /T L denoting the corresponding Carnot efficiency (and having set T L = T H and T R = T C ). A crucial identity that is used to derive Eq. (2) is the expression of the entropy production rate in terms of the output power and of the efficiency, i.e.
The bound Eq. (2), of paramount importance due to its wide applicability to many systems ranging from colloidal systems [62] to biological clocks [19], is however expected to fail whenever quantum systems, such as nanoscale heat engines, cannot be suitably described by effective Markovian processes. Recent results have in fact witnessed violations of the classical TUR Eq. (1) [22] and of the bound on the output power Eq. (2) in paradigmatic toy models, such as resonant single-dots and serial (or side-coupled) double-dot junctions [23].
A straightforward application of our new bound Eq. (21), when restricted to Eq. (23) withĴ α =Ĵ w , leads to This result extends the validity of the conclusions obtained in Ref. [21] that were summarized in the introduction of this work. What is more remarkable, however is that Eq. (25) indicates that the allowed output power for given engine efficiency and constancy (i.e. power fluctuations) can potentially be four times larger than any counterpart described as a classical Markov stochastic process. It can moreover be easily checked that all the violations observed in the above mentioned toy models analyzed in Refs. [22,23] are well within our new bound, even in the presence of Coulomb interaction between the quantum dots [23] (see also SI, Subsection B.4). On top of this, an additional bound can be derived by exploiting the full covariance matrix ∆Ĵ α ,Ĵ β of Eq. (21). If we consider in particularĴ α =Ĵ w andĴ β =Ĵ L (the latter being by convention the heat current from the hot reservoir), we have that (see SI, Section B.4) where ∆P ,ĴL is the normalized covariance betweenP and J L and ∆Ĵ L the normalized variance ofĴ L (see Eq. (22)). This new upper bound complements the one of Eq. (25) and shows an unexpected relation between the maximum amount of power output and the incoming heat current from the hot (left) reservoir. Eq. (26) implies that when these two quantities becomes highly correlated, i.e. the denominator goes to zero). This is the case, for example, in tight coupling regime, where the heat current becomes proportional to the particle current [23,63]. Since, by definition, η = P / Ĵ L , this means that, for strongly correlated systems, P ∝ √ ∆ P and likewise for the heat current from the hot reservoir. Such relation between the mean values and their variances is typical, e.g., of Gaussian distributions and is expected to hold for ergodic systems.

CONCLUSIONS AND DISCUSSIONS
In this paper we have explored the thermodynamics of precision for NESS. We exploited a statistical ensemble description and an expression of the entropy production in terms of relative entropy in order to bound the dissipation from below by the covariance matrix of currents. Our result differs from the standard approaches in the literature for the thermodynamics of precision -not only is it derived in a fully quantum mechanical way, it also is geometrical in nature, reflecting the underlying universality of the concept. It is important to stress that the derivation of our result in Eq. (21) does not assume any Markov approximation and it is valid at second order in δ β and δ βµ , thus beyond the linear response regime. Moreover, as it is the case for the classical TUR, it holds true for any current in the steady-state system. However, it also goes further in that it contains information on the covariance between different currents, whereas the usual TUR that instead concern only the variances of currents, i.e. the diagonal elements of the above covariance matrix. Employing our bound in the context of mesoscopic engines, allows us to demonstrate that a machine not modelled with a classical Markovian description can be more powerful. We speculate, from our example that this is due to the presence of quantum coherence. A detailed exploration of this speculation is a study which we are currently undertaking. In addition in future work we plan to investigate the repercussions of our bound on the precision of quantum clocks [64]. In this section we present some technical details concerning the derivation of the McLennan-Zubarev generalized statistical ensemble [Eq. (7) of the main text], which will be necessary for the derivations of our main results. We begin by considering the unitary evolution of a system described subject to a Hamiltonian of the formĤ =Ĥ 0 +Ĥ int . The Schrödinger picture density matrixρ S (t) will then evolve from some initial time t 0 according toρ S (t) =Û (t, t 0 )ρ S (t 0 )Û † (t, t 0 ), whereÛ (t, t 0 ) = e −iĤ(t−t0) . We move to the interaction picture with respect toĤ 0 by definingρ I (t) =Û † 0 (t, 0)ρ S (t)Û 0 (t, 0), whereÛ 0 (t, 0) = e −iĤ0t . Note that here we have chosen the time t = 0, and not t 0 , as the coincidence time between operators and states in the two pictures.
The time evolution ofρ I (t) will then be given byρ The expectation value of an arbitrary observableÂ not explicitly dependent on time is given by Â (t) = Tr Â I (t)ρ I (t) , whereÂ I (t) =Û † 0 (t, 0)ÂÛ 0 (t, 0). In the particular case whereρ S (t 0 ) commutes withĤ 0 , thenρ S (t 0 ) =ρ I (t 0 ) =ρ 0 and Â (t) simplifies to Following standard literature, the steady-state expectation value ofÂ is then defined as the asymptotic limit Â ness = lim |t−t0|→∞ Â (t) . We will in particular consider the case where t = 0 and t 0 → −∞. On the one hand, it can be proven that this choice corresponds to taking the correct causal constraint rather than the steady-state corresponding to the advanced solution [40,42]. On a more intuitive basis however, one can also motivate this choice on physical grounds by arguing that, when dealing with a steady-state system, one wants to calculate thermodynamic quantities in the steady-state, i.e. once the latter is established. This choice is then also consistent with choosing the coincidence time for the quantum mechanical pictures for the evolution to coincide at time t = 0, where the steady-state is assumed to be reached. One therefore obtains that where we have introduce the Møller operatorΩ We will now specialize this to the McLennan-Zubarev NESS operator [Eq. (7) of the main text], and show how it can be constructed starting withρ 0 as given by Eq. (6) of the main text. To derive the explicit form ofρ ness used in Eq. (7) of the main text, we next introduce the notion of adiabatic switching of the interaction [41].
To that end, we distort the original Hamiltonian to read with ǫ being a positive infinitesimal number and g is a dimensionless bookkeeping parameter that can be formally set to unity at the end. This new Hamiltonian smoothly interpolates between the free HamiltonianĤ 0 at |t| → ∞ and the total HamiltonianĤ at t = 0. The adiabatic limit corresponds to ǫ → 0 + , which should be taken only in the end of all calculations. One may now directly verify that for 0 ≥ t ≥ t 0 , the new interaction picture evolution operatorÛ ǫ,I (t, t 0 ) satisfies the differential equation [65] iǫg whereĤ ǫ,I =Û † 0 (t, 0)Ĥ ǫÛ0 (t, 0). Specializing to the case where t = 0 and t 0 = −∞ and using the fact thatĤ ǫ (0) =Ĥ and H ǫ (−∞) =Ĥ 0 , we get Taking the limit ǫ → 0, so the l.h.s. vanishes, and identifying the Møller operator (S4) aŝ then finally leads to the so called intertwining propertyĤ This somewhat unintuitive result, which is equivalent to the Gell-Mann and Low theorem, shows that the Møller operators can actually be used to formally connect the free and full Hamiltonians.
We now rearrange the different terms as follows. First, we defineβ = (β which are related to the asymptotic values of the energy and particle imbalances. We may then write, for instance, where we used the fact thatĤ L,+ +Ĥ R,+ +Ĥ C,+ =Ĥ is the full Hamiltonian. A similar result holds for the particle number operators. With this rearrangement, Eq. (S9) may now be written aŝ ρ ness = 1 Z L Z R e −β(Ĥ−µN )−δ βÊ+ +δ βµQ+ e β(ĤC,+−µNC,+)ρ C,+ . (S11) As shown in Refs. [34,42] the stateρ ness is independent of the initial stateρ C of the central system. This proof, however, is not trivial. Instead, we shall assume for simplicity that the initial stateρ C iŝ which results in the last two terms in Eq. (S11) cancelling out since e β(ĤC,+−µNC,+)ρ C, Consequently, Eq. (S11) then becomes Eq. (7) of the main text, specificallŷ where Z ness is the normalization constant. It is finally worth mentioning that, by exploiting the Dyson's expansion ofÛ ǫ,I (0, −∞) and Abel's theorem [40], it is possible to expressÊ + andQ + as time-averaged Heisenberg-picture operators aŝ withX H (t) = e iĤtX e −iĤt . We refer the interested reader to Ref. [42] for all the details In light of the formalism illustrated in detail in the previous Subsection, the NESS statistical operatorρ ness is reached at time t = 0 through an adiabatic switching on of the coupling at initial time t 0 = −∞, where the initial state wasρ 0 . What we defined to be "entropy production operator"Σ = δ βµQ+ − δ βµÊ+ appearing at the exponent of the NESS state Eq. (7) is therefore, by construction, a quantity which expresses the dissipated work necessary to create the steady-state starting from the factorized stateρ 0 . It is therefore far from evident why its average over the NESS should correspond to the usual expression considered in the steady-state, where Σ = lim T→∞ T σ with σ being given by [66] where the affinities δ β and δ βµ have been defined previously and Ĵ Q,E are the steady-state values of the particle and energy currents, respectively. The latter are in fact known to be constant in time and all the above mean values are assumed to be taken with respect to the NESS state, i.e. · ≡ Tr [ ·ρ ness ]. Put in another way, can one prove that the mean value of that operator, i.e. the average dissipated work to create the NESS, is actually equivalent to the average entropy production generated in the steady-state (i.e. once the NESS is obtained)? The answer to this question is affirmative, and to prove this important fact the first step is to notice that the average currents (particle and energy) in the NESS are time independent. This was explicitly shown in Ref. [42] but, in light of its importance, we will repeat here some of the main steps using our notation for convenience. Let us consider first the expectation value of the particle current at a generic time t > 0, i.e. once the NESS is established, see Fig. S1 whereÛ ǫ (t, 0) is the evolution operator corresponding to the HamiltonianĤ ǫ (t) and where we have used the notation Note that the steady-state solution for the statistical operatorρ ness Eq. (7) is obtained fromρ ǫ by taking the adiabatic limit, i.e. ρ ness = lim ǫ→0 +ρ ǫ . Finally, in what follows, we will also use the alternative notationρ T −1 to equivalently denoteρ ǫ after we have switched the limits lim ǫ→0 + into lim T→∞ using Abel's theorem (as in Eq. (10)). We recall the very important fact that the thermodynamic limit lim ǫ→0 + , which was taken before to get the final closed expression forρ ness , must always be performed only at the end of the calculations (i.e. in this case after the average is taken). By making use of the identity Eq. (10), one has that where we have defined the current operator (in Heisenberg picture) It is then possible to show that (see Eqs. (31) to (34) of Ref. [42] for the details), upon definingĴ ′ Q (τ, t) ≡Û ǫ (τ, 0)Ĵ Q (t)Û † ǫ (τ, 0), one can compute the second term in the above expression which, substituted back into Eq. (S18), gives Analogous calculations hold for the mean steady-state energy current Ĵ E . Equipped with these results, we can now consider the entropy production operatorΣ defined in Eq. (8). We have that This concludes the proof, as we can evaluate the mean entropy production rate in the NESS as the latter being the average of the entropy production operator appearing in the NESS statistical operator.
In order to prove Eq. (14) it is useful to keep in mind that, as stated in the main text (and proven e.g. in Refs. [36,67]), the timeaveraged entropy production operatorΣ is a conserved quantity and therefore it commutes with the total HamiltonianĤ (and consequently also with the total number operator N ). This allows to 'disentangle' the exponential exp −β Ĥ − µN +Σ = exp −β Ĥ − µN exp Σ . We have therefore that where the first term comes from the m = 0 term, i.e. the identity 1, which therefore gives back Tr [ρ les ] = 1, and the remaining terms of the series have been sorted into odd and even powers ofΣ. Since the expectation value is taken with respect to the LES statistical operator, the former (i.e. the odd powers of the entropy production operator) vanish and only the even powers survive, thus leading to Eq. (14).

Subsection B.3 -Proof of our bound on thermodynamic precision
In this Subsection we will provide the explicit derivation of Eq. (21). As explained in the main text, the starting point is to perform the following transformation on the manifold of steady-states (SSM) Finally, one can express the above bound in terms of the average steady-state entropy production rate, given the relation Reminding the reader of the very important fact that the thermodynamic limit lim T→∞ must be performed only at the very end of calculations, one then introduces, retracing the exact same steps as done in standard literature of TUR, the normalized covariance matrix and simplifies the time T from the expression before taking the thermodynamic limit. It is straightforward to show that this immediately leads to Eq. (21). Finally, it is worth stressing that, from the above result which generalizes the TUR as it involves the full covariance matrix, it follows in particular that the diagonal elements must be positive as well, from which one obtains which has the same structure as the classical TUR but a bound four times looser than the Markovian counterpart.

Subsection B.4 -The double serial quantum dots steady-state engine
We will devote this Subsection to briefly show the application of the TUR-derived upper bound on power in a toy model considered in [22,23]. In particular, we will choose the serial double quantum dots junction model since it has been shown to manifest violations of the TUR even at arbitrary small biases δ β and δ βµ when a second order expansion is considered. We remind in fact that this corresponds in fact to the regime of validity of our geometrical new TUR.
Let us therefore consider a 1D junction system made of two quantum dots, with energies E L,R , coupled to each other coherently through a tunnelling amplitude Ω. The two dots are then respectively hybridized with their corresponding lead with a tunnelling amplitude t α and chemical potential µ α . The total Hamiltonian is then where {ĉ a ,ĉ † a } are the annihiliation and creation operators for the quantum dots, while {b a,k ,b † a,k } are those of the fermionic leads for an eigenstate k with energy ǫ k . We also assume, without any loss of generality, that T L > T R .
In the main text we showed how the TUR implies an upper bound to the power, dictated by its fluctuations ∆ P and by the efficiency ǫ, according to Eq. (25). The quantities entering this bound are usually calculated using the formalism of non equilibrium Greens function [45,[59][60][61] and concretely are given, in the wideband limit, by the Landauer-Büttiker formulas with f (E) = e β(E−µ) + 1 −1 being the Fermi-Dirac distribution and τ (E) being the transmission function, which we will assume of the form where Γ a = 2π|t a | 2 d a is the real part of the dot's self-energy quantifying the coupling of lead a to dot a, with d a denoting the density of states of the lead. Finally, the efficiency of the steady-state engine is given by ǫ = P / Ĵ L , with denoting the heat current from the left reservoir, and In Fig. S2 the power (blue solid curve) is displayed as a function of the two quantum dots' detuning E L − E R = ∆E. In the resonant case ∆E = 0 a Markovian description fails and a violation of the classical TUR arises because of the degeneracy in the system's Hamiltonian. For increasing values of ∆E the validity of B P S (red curve) is recovered as sequential hopping becomes dominant. These results are in line with those obtained in Refs. [22,23] for the same model. Since the violation of the classical TUR in this setup is less than 1% of B P S it well below our new upper bound B GG demonstrating that it holds true even in known quantum regimes. Comparisons with the results obtained in Ref. [23] for the violations of B P S in presence of a Coulomb interaction term in the Hamiltonian Uĉ † Lĉ Lĉ † Rĉ R show that also in this case they are well within the predictions of our bound B GG .

Subsection B.5 -Proof of the new lower bound on power
Let us start from the re-expression of the entropy production rate in terms of the power and of the efficiency Let us chooseĴ = Ĵ w ≡P ,Ĵ L T (the latter component being by convention the heat current from the hot reservoir). The inverse of the normalized covariance matrix can be calculated explicitly using the following relation, true for any square n × n matrix M with C being the square matrix of cofactors of M, i.e. C ij = (−1) i+j m ij with m ij being the minor of M obtained deleting the ith row and jth column. The result is given by where we have defined, in conformity of notation with the main text, the normalized correlation function between the power and the heat current from the hot (left) reservoir and the normalized variance ofĴ α . The direct application of Eq. (21), making also use of the definition of the efficiency η ≡ P / Ĵ H , leads straightforwardly to Eq. (26)