Work Fluctuations in Ergotropic Heat Engines

We study the work fluctuations in ergotropic heat engines, namely two-stroke quantum Otto engines where the work stroke is designed to extract the ergotropy (the maximum amount of work by a cyclic unitary evolution) from a couple of quantum systems at canonical equilibrium at two different temperatures, whereas the heat stroke thermalizes back the systems to their respective reservoirs. We provide an exhaustive study for the case of two qutrits whose energy levels are equally spaced at two different frequencies by deriving the complete work statistics. By varying the values of temperatures and frequencies, only three kinds of optimal unitary strokes are found: the swap operator U1, an idle swap U2 (where one of the qutrits is regarded as an effective qubit), and a non-trivial permutation of energy eigenstates U3, which indeed corresponds to the composition of the two previous unitaries, namely U3=U2U1. While U1 and U2 are Hermitian (and hence involutions), U3 is not. This point has an impact on the thermodynamic uncertainty relations (TURs), which bound the signal-to-noise ratio of the extracted work in terms of the entropy production. In fact, we show that all TURs derived from a strong detailed fluctuation theorem are violated by the transformation U3.


Introduction
A quantum description of thermodynamic heat engines has lately become necessary to consider physical systems at the mesoscale and nanoscale [1][2][3], such as nanojunctions thermoelectrics [4], quantum dots [5], and biological [6,7] or chemical [8] systems.The optimal transport theory also has recently been embedded in a thermodynamic quantum framework [9].At the quantum level, the fluctuations of the thermodynamic variables play a fundamental role, due to the discrete spectral structure of quantum systems.
Figure 1.Scheme of a quantum thermodynamic engine based on the two-stroke Otto cycle with two qutrits as working fluid.In the first stage, the qutrits A and B with frequency ω A and ω B are at thermal equilibrium with the corresponding baths at temperature T A and T B , respectively, with T A > T B .In the second stage, the two systems are isolated and allowed to interact through a unitary evolution extracting work W. Finally, in the last stage, the systems A and B are allowed to relax to the corresponding thermal baths, implying that A absorbs the heat Q H and B releases the heat Q C , thus restoring the initial condition.
In the case of an engine based on a two-stroke Otto cycle the full probability distribution of work and heat has been retrieved for two qudits [35] and for two bosonic modes [36] as working fluids, where the transformation for the work extraction is the unitary partial-swap interaction.The two-stroke Otto engine is particularly interesting with respect to its well-known four-stroke version because it allows to extract the maximum amount of work in the adiabatic step of the cycle by a single unitary operation, the so-called ergotropy [49][50][51][52][53].Note that the extraction of the ergotropy necessarily depends also on the transformation that couples the systems.We will show here that, if the systems are qudits with dimension larger than two, unitary evolutions different from the swap interaction can increase the extracted work.
We will define a procedure for determining the unitary interaction that provides the maximum work from two multilevel systems A and B for a given choice of the relevant parameters, i.e. the frequency gaps ω A and ω B of the qudits and the temperatures T A and T B of the reservoirs.Then we take the specific case of a working fluid described by two qutrits and classify all the transformations that extract the ergotropy.Specifically, we find three different kinds of optimal unitary strokes: the swap operator U 1 , an idle swap U 2 (where one of the qutrits is regarded as an effective qubit), and a non trivial permutation U 3 given by a composition of the two previous unitaries, namely Each transformation extracts the ergotropy from a different regime defined by the frequency gaps ω A and ω B of the two qutrits and by the temperatures T A and T B of the baths.By deriving the characteristic function of work and heat, we evaluate the work statistics and the entropy production for every case.Note that a complete description of a quantum ergotropic heat engine and of the procedure for determining the work statistics is detailed in Appendix A. Then, we focus on the trade-off between ergotropy extraction and relative fluctuations var(W)/⟨W⟩ 2 , i.e. the inverse of the signal-to-noise ratio (SNR).The evaluation of the fluctuations allows us to establish the relation between the variance of the work and the mean entropy production in terms of the TURs.A standard reference TUR bounds the fluctuations with the inverse of the entropy production as follows [29] var(W) We show that all the three ergotropic transformations violate this TUR.Moreover, U 3 is proved to beat all the TURs derived by the strong fluctuation theorem where the forward and backward processes in Eq. ( 1) are related by the same condition The paper is structured as follows.In Section 1 we define the procedure for determining the transformations extracting the ergotropy in the case where the working fluid is described by two qudits with generic dimensions d A and d B .Then, in Section 2 we apply our procedure to the case of two qutrits.In particular, in Section 2.1 we classify all the transformations extracting the ergotropy and their properties.In Section 2.2 we evaluate the maximum work extracted by each transformation in terms of the frequency gaps and the temperatures.In Section 2.3 we study the mean entropy production related to each interaction.In Section 2.4 we derive the work distributions.Finally, in Section 2.5 we find the relative fluctuations of work, compare the corresponding SNR to the bounds provided by the most relevant TURs and discuss the assumptions required for these TURs to hold.In Section 3 we draw our conclusions.

Materials and Methods
In this work we fix the Planck and Boltzmann constants to natural units, i.e. h = k B = 1.We consider two qudits A and B in a product of Gibbs states, i.e. where n=0 n|n⟩⟨n| is the Hamiltonian of the system X = A, B, each one with equally-spaced energy levels, and Z X = Tr[e −β X H X ] denotes the corresponding partition function, and β X = T −1 X the inverse temperature.The number states |n⟩ in the expansion of the Hamiltonians are eigenstates of the occupation number n X ≡ H X /ω X .Without loss of generality, we fix T A > T B .We use the state in Eq. ( 3) as input to a two-stroke Otto engine.As depicted in Fig. 1, the process starts with the two qudits in thermal equilibrium with their baths, at temperature T A and T B .Afterwards, the two qudits are isolated from their baths and we make them interact through a unitary evolution in order to extract the ergotropy.The procedure for the ergotropy extraction will be detailed in the following.Once the work has been extracted through the interaction, we reset the two qudits to their equilibrium state by re-connecting them via weak coupling to their thermal baths.The sequential repetition of this process leads to our two-stroke cyclic engine.We fix the convention of positive work for the extraction from the system and positive heat for the absorption from the reservoirs.Then, in each cycle the average energy change in system A due to the unitary stroke corresponds to the average heat released by the hot reservoir A, namely ⟨Q H ⟩ = −⟨∆E A ⟩. Similarly, for the cold reservoir, ⟨Q C ⟩ = −⟨∆E B ⟩, and, for the first law of thermodynamics, the average work is given by ⟨W⟩ = ⟨Q H ⟩ + ⟨Q C ⟩ = −⟨∆E A ⟩ − ⟨∆E B ⟩. Correspondingly, the average entropy pro- Our goal is the investigation of an ergotropic heat engine based on the two-qudit system described above, i.e. an engine extracting the maximum work by exploiting the difference of frequency and temperature between the systems A and B. In other words, we are looking for the unitary transformations U mapping the input ρ 0 into a state ρ = Uρ 0 U † such that the average extracted work is maximized, i.e.

⟨W⟩ = max
where here H = H A ⊗ I B + I A ⊗ H B is the Hamiltonian of the system.The evolution that extracts the ergotropy was identified in Ref. [49] as the one minimizing the final energy Tr [ρH].In the present case, where the initial state ρ 0 has no coherence, namely it is diagonal in the energy basis, the ergotropic evolution is the transformation that permutes the eigenstates of the input state so that the magnitude order of the energy levels is reversed with respect to the corresponding occupation fractions.More explicitly, if we take the occupation fractions of the system e −(nβ A ω A +mβ B ω B ) /(Z A Z B ) in descending order, the transformation permutes the related eigenstates to set the corresponding energy levels in ascending order.
If the input state already displays this configuration, then the state is called passive and no transformation can extract work.In summary, since unitary transformations preserve the spectrum, the ergotropy is extracted by reversing all possible population inversion with respect to the energy levels.In the following, we provide a re-visited analysis of the first level maximization strategy developed in Ref. [50].
The procedure of ergotropy extraction can be formalized in a compact way for two subsystems A and B of dimension d A and d B as follows.We consider two different permutations P E and P ρ of the energy eigenstates with respect to their lexicographic order.The permutation P E sorts them so that the corresponding eigenvalues are set in ascending order, i.e.
where the vector of eigenvalues . Similarly, the permutation P ρ rearranges the occupation numbers of the initial state in descending order, namely, . Then we can straightforwardly find the transformation which minimizes the final energy from implying that the ergotropic transformation can be expressed as For instance, take two qubits in a Gibbs state Then the energies pertaining to the levels |10⟩⟨10| and |01⟩⟨01| are ω A and ω B , respectively, while the related occupation fractions are Z or the symmetric case where both the order relations are reversed, the transformation that swaps |10⟩⟨10| with |01⟩⟨01|, namely U = U † = |00⟩⟨00| + |11⟩⟨11| + |01⟩⟨10| + |01⟩⟨10|, extracts the ergotropy.This result appears immediately if we consider the permutation matrices P E and P ρ , which in this case read where θ(x) is the Heaviside function.The operator P † E P ρ promptly identifies the ergotropic transformations and the corresponding ergotropic regimes, since This simple example shows how the extraction of the ergotropy is entirely determined by the order relations between the parameters.In particular, the initial state of an equallyspaced two-qudit engine is described for any dimension of the qudits by a first partial order over the frequencies ω and a second one over the products βω.These order relations identify four basic partially-ordered sets (posets).In the two-qubit example, the ergotropy can only be extracted if the initial state belongs to where the bar denotes the same poset with A and B switched.The states belonging to the remaining two sets are passive.
The description in terms of posets becomes more complex in higher dimension.For a state as in Eq. ( 3), the ordering procedure for the ergotropy extraction needs to establish if kω A > jω B and if Even if the simplest non-trivial case would be a system made of a qubit and a qutrit, here, as mentioned above, we consider a two-qutrit system, so that we can use the results for the two-stroke swap Otto engine with two qudits with equal dimension studied in Ref. [35] as a benchmark.In this scenario, each of the four basic posets mentioned above is further partitioned in four subsets, defined by the order relations 0 < y X 1 < y X 2 /2 and y X 2 /2 < y X 1 < y X 2 , with y = ω or βω and X 1 ̸ = X 2 may be A or B. The total number of posets determining the regimes for the ergotropy extraction is then sixteen.We expect some of them to be passive regimes, i.e. the input state defined by those parameters is passive.
As for the others, we will show that a specific transformation can extract the ergotropy from different regimes, as we noted for the two-qubit case with the swap in the regimes Ω and Ω.

Ergotropic transformations
As mentioned above, we can jointly classify all the ergotropic transformations U and the corresponding ergotropic regimes by inspecting the permutations P E and P ρ .In the two-qutrit case we have four posets identified by ω A and ω B for P E , and four identified by β A ω A and β B ω B for P ρ .We find different permutations P E and P ρ for each of the corresponding four posets, i.e. four distinct transformations.We show them associated to the corresponding poset in Fig. 2. Note that, for what concerns P ρ , we have to distinguish three inequivalent cases identified by the relative position of points on the ω B axis according to the value of the ratio β A /β B .In summary, P E and P ρ are simply the identity I (i.e.no reordering is needed) for ω A < 2ω B and β A ω A < 2β B ω B , respectively.For ω B > 2ω A and β B ω B > 2β A ω A , both P E and P ρ are given by the swap U 1 , namely or, equivalently, U 1 = (24)(37)(68), using the cycle notation and the lexicographic ordering where the elements of the cycles are related to the kets as |nm⟩ → 3n + m + 1.For both P E and P ρ are given by both P E and P ρ are given by We notice that U 2 and U 3 are not invariant under swap symmetry.In particular, Ũ2 ≡ The unitary operators U 1 , U 2 and Ũ2 are also Hermitian and hence self-inverse.Notice also that and, similarly, Ũ3 The product P † E P ρ together with the composition rules for U 1 , U 2 and U 3 explored above allows to find the ergotropic transformations for each ergotropic regime identified by combining an ω poset with a βω poset.In particular, we remark that the ergotropic transformations resulting from the product P † E P ρ must be again U 1 , U 2 , Ũ2 and U 3 .They are five overall, considering the identity too, which pertains to initial passive states.We provide a direct visualization of the landscape of ergotropic transformations in Figs.3-5.Having set β A < β B , each figure is linked to a different regime for the ratio β A /β B .As outlined in Fig. 2, we can identify three distinct ranges of β A /β B with two critical values, namely 1/4 and 1/2.For each case, we show the ergotropic transformation related to each poset.In particular, we set β A /β B = 1/16 in Fig. 3, β A /β B = 5/16 in Fig. 4 and β A /β B = 9/16 in Fig. 5. Firstly, we observe that in the first two cases all the transformations found above appear (except Ũ3 , which pertains to the regime T A < T B ).In the case of Fig. 5, U 3 is never present and the number of passive regimes becomes four.Notice that in the region 0 < β A /β B < 1/4, it is possible to take the limits β A → 0 and β B → ∞.In this case, one of the passive regimes disappears and most of the parameter region is dominated by the swap.On the contrary, approaching the critical point β A /β B = 1/4 we see that the region where the swap extracts the ergotropy shrinks until it vanishes at the critical point.In the second case, in Fig. 4, the swap plays again a role, but the passive regimes grow as well until, at the critical point β A /β B = 1/2, the ergotropic region of U 3 vanishes and is replaced for Let us inspect more in detail the non-trivial ergotropic transformations U 1 , U 2 , Ũ2 and U 3 .The swap U 1 clearly commutes with the total number operator, namely On the other hand, the evolutions U 2 and Ũ2 act asymmetrically on the two systems, since they perform a permutation of the frequency levels of ρ 0 as if the system identified by the smallest frequency gap (B when the ergotropy is extracted by U 2 and A when it is extracted by Ũ2 ) were a two-level system, being its intermediate level |1⟩ left unaffected.Thus, we name U 2 as idle swap.In fact, for this asymmetry we have U 2 ̸ = Ũ2 .Differently from U 1 , the idle swaps U 2 and Ũ2 enjoy the conservation laws As for U 3 = U 2 U 1 , being the composition of the standard and the idle swap, we name it double swap.We noticed above that U 3 is not Hermitian.Indeed, one finds out that the double swap has multiplicative order six, namely U 6 3 = I, as it can be inferred from the cycle notation in Eq. ( 14).Furthermore, the double swap does not commute with any linear combination of n A and n B .In Appendix A, we prove that, if the transformation commutes with a linear combination of H A and H B , then all work and heat moments are proportional each another and hence the mean entropy production is proportional to the mean extracted work, as we will explicitly show for U 1 , U 2 and Ũ2 in the next Sections.

Ergotropy
Now we are ready to provide the mean work of Eq. ( 4) extracted by each ergotropic transformation.In the case of the swap U 1 , the ergotropy can be expressed in terms of ω A − ω B units and reads In the case of the idle swaps U 2 and Ũ2 we obtain and Here we recognize the action described above: the lower frequency qutrit is taken as a qubit whose gap is 2ω B for ⟨W 2 ⟩ and 2ω A for ⟨ W2 ⟩, so that the extracted work is proportional to ω A − 2ω B and ω B − 2ω A , respectively.As expected, the work extracted from U 2 is obtained from the one extracted by Ũ2 just by swapping A with B. From Eqs. (20)(21)(22) one also verifies that where the ratio x ≡ ω B /ω A is a relevant parameter, as we will find in the following.
In the case of the double swap, we have Here we see the effects of the atypical behavior of U 3 : the extracted work is not proportional to any frequency gap.On the contrary the frequencies ω A and ω B appear multiplied with different weights.Notice that for x = 1/2 one has ⟨W 2 ⟩ = 0 and from the second line of Eq. ( 24) the double swap U 3 extracts the same work as U 1 , i.e. ⟨W 3 ⟩ = ⟨W 1 ⟩.Instead, for x = 1, namely ω A = ω B , one has ⟨W 1 ⟩ = 0 and ⟨W 3 ⟩ = ⟨ W2 ⟩.Finally, for x = 2 we have ⟨ W2 ⟩ = 0 and again ⟨W 3 ⟩ = ⟨W 1 ⟩.In Fig. 6 we represent the ergotropy extraction in the case β A /β B ∈ (0, 1/4).In particular, we set the ratio β A /β B = 1/16, as in Fig. 3, with β B = 10.Note that the pretended discontinuities in the transitions between different ergotropic regions are just cusps, as it can be recognized in Figs.7 -11.In these figures, we show specific examples of ergotropy extraction as a function of ω B , by fixing all the other parameters.Fig. 7 displays the case β A /β B < 1/4, with β A /β B = 1/8.Therefore, this is not a critical point and for varying ω B we span all the non-equivalent ergotropic transformations.The black dot line displays the work extracted from the standard swap U 1 , so that we can see how it is outperformed by the other unitaries outside its own ergotropic regime.Moreover, the solid lines, corresponding to U 2 and Ũ2 , show that the regime of operation of an ergotropic heat engine is enlarged with respect to the swap Otto engine.In Fig. 8 we consider the critical point β A /β B = 1/4, which represents the transition between the cases in Fig. 3 and Fig. 4, where the ergotropic regime of the standard swap vanishes.Indeed, here we do not have any ergotropic contribution from U 1 , except for the limiting case ω A = 2ω B , where the work extracted coincides with the one provided by U 3 , identified by the red point.In Fig. 9 we show the ergotropy as a function of ω B for the critical point β A /β B = 1/2, which is the transition point between the cases of Fig. 4 and Fig. 5.As expected, the double swap U 3 is never required to extract the ergotropy.The maximum work is extracted by the idle swap U 2 for ω B < ω A /2, by the standard swap U 1 for ω A /2 < ω B < ω A and by Ũ2 for ω A < ω B < 2ω A .For the case β A /β B > 1/2 of Fig. 5, we fix in Fig. 10 β A /β B = 3/4.As in the previous case, U 3 is not needed and, furthermore, there are two more passive regions.Finally, in the last example in Fig. 11, we plot the ergotropy for the ideal case β A /β B = 0, by setting β A to 0 and finite large values for ω A and β B .In particular, the high value of ω A allows to see that the extracted work is large when ω A − ω B is large, except for the limiting case ω B → 0 (in such case indeed we would have     In summary, in the regime of operation of the standard swap Otto engine, i.e.
the work extraction may be improved by replacing the swap U 1 with the permutation U 3 .Moreover, the idle swaps U 2 and Ũ2 even allow to enlarge the range of operation of the heat engine.

Entropy production
Let us now evaluate the mean entropy production of the quantum heat engine in order to study its relation with the work fluctuations and to explore the validity or violation of TURs.As mentioned in Section 1, the mean entropy production is given by ⟨Σ⟩ We can evaluate the moments of W and ∆E A through the derivatives of the characteristic function, according to Eqs. (A27) and (A28) of Appendix A. Due to the conservation laws for U 1 , U 2 and Ũ2 as in Eqs. ( 18) and ( 19), according to Eq. (A34), we have where for Ũ2 .Hence, the entropy production of U 1 , U 2 and Ũ2 is proportional to their pertaining work, and one has where ⟨W 1 ⟩, ⟨W 2 ⟩ and ⟨ W2 ⟩ are given in Eqs. ( 20), ( 21) and (22), respectively.Eq. ( 26) does not hold for U 3 , and the entropy production explicitly is given by Note that in all cases the mean entropy production is positive and depends only on the ratios between frequency and temperature and not on the bare frequencies.

Work distribution
We can now provide the explicit expression for the distribution of work p(W) pertaining to each ergotropic transformation.As shown in Appendix A [see Eq. (A15)], we have where p n,m is the energy distribution of the input state, namely while q(l, s|n, m) is the energy conditional distribution after the evolution U, given the input energy levels n and m, i.e.
In the case of the standard swap U 1 , the conditional distribution reads q 1 (l, s|n, m) = δ l,m δ n,s , and hence which is a 5-point distribution.Explicitly, upon naming k ≡ n − m, one has A specific example is plotted in Fig. 12. Eq. ( 33) is consistent with the general result given in Ref. [35] for the work distribution in swap engines based on two qudits.Now we focus on the idle swap U 2 .Due to its asymmetric action on systems A and B, the conditional distribution is slightly more complicated and reads where denotes the sum mod 3. Hence, one retrieves the following 3-point distribution An example is depicted in Fig. 13.
Finally, since U 3 = U 2 U 1 , for the double swap we readily find and then one obtains the following 7-point distribution A specific example is provided in Fig. (14).

Work fluctuations and TURs
Here we evaluate the relative fluctuations of the work extracted by the ergotropic transformations and compare them to the lower bounds identified by different thermodynamic uncertainty relations (TURs).We can find the relative fluctuations as the ratio between the variance of the extracted work and the square of its mean value, namely var(W)/⟨W⟩ 2 = ⟨W 2 ⟩/⟨W⟩ 2 − 1, with var(W) = ⟨W 2 ⟩ − ⟨W⟩ 2 .The second moment of the extracted work can be obtained from the characteristic function as in Eqs.(A28) and (A29), and one has For the standard swap U 1 one obtains var(W 1 ) which is in agreement with the general result of the swap engine with two qudits of Ref. [35].As expected, the fluctuations of the standard swap are invariant under swapping A and B. For the idle swap U 2 , we have (41) As for the ergotropy in Eqs. ( 21) and ( 22), the expression for var( W2 )/⟨ W2 ⟩ 2 is simply obtained by exchanging A with B in Eq. (41).Note that the fluctuations of both the standard and the idle swap depend only on the products βω.This is not the case for the double swap U 3 , which depends also on the frequency ratio x = ω B /ω A as follows var(W 3 ) For all the ergotropic transformations the fluctuations are minimized in the limiting case where βω → 0 for one qutrit and βω → ∞ for the other one.In the case of the swap, being naturally invariant under swap symmetry, we can either set β A ω A to zero and β B ω B to infinity or the other way around.On the contrary, the case of the idle and the double swap is asymmetric and we achieve the minimum of the fluctuations for Here we mainly focus on the transformations that extract the ergotropy in the same poset identified by the products βω.In particular, we choose the poset defined by β A ω A < β B ω B , where the optimal evolutions are U 1 , U 2 and U 3 .In the case of the double swap U 3 , the minimization has to be performed also on the frequency ratio and the infimum is obtained for x → 0. The optimization of the fluctuations over the whole span of the parameters readily provides Therefore, it turns out that U 2 and U 3 achieve smaller fluctuations than U 1 .
We will now investigate if damping the noise comes together with the extraction of the ergotropy.While for U 1 this is always the case, the same is not true for U 2 and U 3 .The idle swap U 2 extracts the ergotropy for β B ω B < β A ω A , where the condition for the minimization of fluctuations corresponding to U 2 does not hold.Interestingly, in that region it is Ũ2 the ergotropy extractor.Within the ergotropic region of U 2 , we need to take For U 3 , on the contrary, the condition on the ratios βω for optimal fluctuations is compatible with the extraction of ergotropy, but with the additional constraint x ≥ 1/2.The minimization over x then sets it to 1/2, and, as discussed after Eq. ( 24), for that frequency ratio ⟨W 3 ⟩ = ⟨W 1 ⟩.To sum up, if we aim to optimize the noise inside the ergotropic regimes of each ergotropic transformation, we find that the best performance is achieved by the standard swap, since 2 = inf In this last regime where ergotropy extraction and minimal noise coexist, we finally note that the standard swap extracts more work than the idle and double swap.In fact, one has We remark that the results found so far do imply that the standard swap is the best operation in terms of fluctuations and extracted work in the optimal limiting case β A ω A → 0 ∧ β B ω B → ∞, but the same does not hold for intermediate values of βω, as we shall see in the following.Now we compare the relative fluctuations of the ergotropic engine in asymptotic and non-asymptotic cases with the bounds derived from the most significant TURs.We recall that the double swap U 3 is not Hermitian.Therefore, as remarked in the Appendix after Eq. (A26), U 3 could violate all the TURs based on the equivalence between forward and backward processes.On the other hand, we already know from previous works [35] that the swap itself breaks the standard TUR in Eq. ( 2).We study in Fig. 15 the violation of the standard TUR as a function of βω in conditions of minimal fluctuations, independently from the ergotropic regime.Namely, in the case of U 1 (red dashed line), U 2 (purple solid line) and U 3 (blue dot-dashed line) the free variable is β B ω B with β A ω A ≪ 1.Just for U 3 , we also need ω B /ω A ≪ 1.We remark that here we are not focusing on the ergotropy extraction, but only on the properties of the evolutions U 1 , U 2 and U 3 in terms of work fluctuations.We find that all the three ergotropic transformations break the standard thermodynamic uncertainty relation.In particular, the violation due to U 3 is impressive.As found in [35], when the evolution is the standard swap the relative fluctuations for the extracted work satisfies var(W) Figure 15.Product of the relative fluctuations with the mean entropy production, which is lower bounded by 2 in the standard TUR of Eq. ( 2) and by a function of the mean entropy in Eq. ( 50).We set The variation of Eq. ( 46) from the standard TUR explains the slight violation found in Fig. 15, where the lower bound from the standard TUR is displayed as a black dotted line.Similarly to U 1 , also U 2 and Ũ2 satisfy Eq. ( 46).In fact, where The fluctuations originated from U 3 , instead, can break the TUR in Eq. ( 46).Such violation stems from the asymmetry of the process described by U 3 , which is not Hermitian.Indeed, we note that a necessary condition for the TURs in Eqs. ( 2) and ( 46) to hold is the equivalence between forward and backward process, i.e. p B (W, ∆E A ) = p(W, ∆E A ).Moreover, note that the double swap is the only transformation whose fluctuations depend also on the frequency ratio, while leaving the mean entropy as a function of just β A ω A and β B ω B .Therefore, in this case we can optimize over a third parameter without changing the lower bound of the TUR.The violation of the TUR in Eq. ( 46) by U 3 can be found also in realistic cases, i.e. even if we do not set the parameters to the values minimizing the fluctuations.Actually, these cases are the most relevant to be considered, not only because closer to experimental applications but especially because they keep into account the ergotropy extraction provided by the different evolutions.For instance, consider the case of Fig. 7, where ω A = 1, β A = 0.5, β B = 4 and ω B is left free.Correspondingly, in Fig. 16 we plot the signal-to-noise ratio (SNR) of the extracted work for each transformation in its ergotropic regime, together with the lower bound of Eq. ( 2) (dotted lines).Firstly, note that the double swap violates the TUR even if we are far from the optimal conditions on the parameters maximizing the SNR.Second, the TUR is violated in both regimes where U 3 extracts the ergotropy (ω B ∈ [1/8, 1/4] ∪ [1/2, 1]).Third, differently from what we found in the case of optimal conditions, U 2 and U 3 can achieve better SNRs than the standard swap U 1 where the ergotropy is extracted.2) and the tight TUR in Eq. ( 50), respectively.The dashed horizontal line highlights the asymptotic limit of the SNR.
The standard TUR is not the only relevant lower bound which we show in Fig. 16.The tightest TUR that cannot be violated by any time-symmetric process was found in Ref. [33] and, applied to the extracted work, reads var(W) where g(x) is the inverse of the function x tanh(x).Therefore, we expect that nor U 1 neither U 2 can violate this TUR, while U 3 in principle can.This is what we see in Fig. 16, where the dot-dashed lines corresponds to the lower bound determined by Eq. ( 50): the SNR identified by the double swap U 3 is the only one that can violate the tight TUR, also within its ergotropic regime.
We focus more in detail on the violation of the TURs above in Figs.17-19, where we plot the SNRs for the three evolutions both for optimal values of the parameters independently from the ergotropy extraction and within the corresponding ergotropic regime.In particular, Fig. 17 displays the performance of the standard swap U 1 .Here we set β A ω A ≪ 1, which implies that the fluctuations are minimized for large β B ω B .As β B ω B increases, the signalto-noise ratio approaches the inverse of the minimal fluctuations, namely 3/2, in agreement with Eq. ( 43).Again, we find a slight violation of the standard TUR (dotted line).
In Fig. 18 we show the performance of the idle swap U 2 where it maximizes the SNR (first panel) and where it extracts the ergotropy (second panel).Therefore, in the former case we set β A ω A ≪ 1 and retrieve the optimization of the SNR for large values of β B ω B , as in Eq. ( 43).In the regime where U 2 extracts the ergotropy, as in Eq. ( 44), we find an optimal SNR approaching 1/2 for β A ω A ≫ β B ω B ∼ 0 and an almost negligible violation of the standard TUR.Neither the standard nor the idle swap violate the tight TUR in Eq. ( 50), displayed as a dashed-dotted line.
The case of the double swap, displayed in Fig. 19, is radically different.If we neglect the conditions for the ergotropy extraction, here we can optimize also over the frequency ratio ω B /ω A and we can set it to zero, while β A ω A ∼ 0, implying that we expect to find the optimal SNR for large β B ω B , as in Eq. (43).Again, the standard TUR is violated, but, compared with the previous cases, the corresponding bound is saturated for larger values of β B ω B , where the SNR approaches its maximum.Most importantly, the tight TUR of Eq. ( 50) is also violated, both when the SNR is optimized (first panel) and when the ergotropy is extracted (second panel).
We also compare the SNR of U 3 with the loosest bound that always holds for timesymmetric processes [21,31,34] given by var(W) The bound from Eq. ( 51) is displayed as a brown line in Fig. 19.The violation that we find is a consequence of the fact that U 3 is not Hermitian.
In the second panel of Fig. 19, as mentioned above, we explore the performance of the double swap U 3 in its ergotropic regime, where ω B /ω A ∈ [1/2, 1].The best performance is obtained for ω B /ω A = 1/2, where the amount of work extracted by U 3 is the same as the one extracted by U 1 (red line in Fig. 19).We also plot the case ω B /ω A = 3/4, in blue.We get a worse SNR, but still can observe a violation of all the TURs derived for time-symmetric processes.
The only TURs that can set a bound which cannot be violated by U 3 are those obtained without posing the symmetry between forward and backward process.In fact, the TURs in Eqs. ( 50) and ( 51) have been generalized respectively in Refs.[39] and [32] by releasing the assumption that forward and backward processes share the same distribution of the stochastic variables.These new bounds are given by var(W) + var(W) B (⟨W⟩ and var(W) + var(W) B (⟨W⟩ where the quantities with subscript B are referred to the backward process and a = (⟨Σ⟩ + ⟨Σ⟩ B )/2.In the case of U 3 , the statistics of W for the backward process is easily found since U −1 3 = Ũ3 .Hence, U −1 3 outputs the same work statistics as U 3 provided that systems A and B are swapped.Then, ⟨W 3 ⟩ B , var(W 3 ) B and ⟨Σ 3 ⟩ B can be obtained from Eqs. ( 24), ( 28) and ( 42) simply swapping labels A and B. Note also that the bounds (right-hand sides) given by the TURs in Eqs. ( 52) and ( 53) depend only on the products βω, while the corresponding bounded quantities depend also on the frequency ratio ω B /ω A .In Fig. 20, we compare the reciprocal of the left-hand sides of Eqs. ( 52) and ( 53) for U 3 with the corresponding bounds as a function of β B ω B with fixed β A ω A ≪ 1.In this regime U 3 maximizes the SNR.We show the two limiting cases ω B /ω A ≪ 1 (thick dark-blue curve) and ω B /ω A ≫ 1 (thin light-blue curve) together with the bounds obtained from the TURs in Eqs. ( 52) and (53), identified by the dot-dashed brown curve and the dashed green curve, respectively.Note that these TURs are never violated and, as expected, the first is tighter than the second.Having set β A ω A ∼ 0, the maximum is asymptotically reached for β B ω B ≫ 1 and ω B ≫ ω A , and amounts to 8/9 (dashed horizontal line in Fig. 20).52) and Eq. ( 53), respectively.The dashed black horizontal line identifies the asymptotic value, which amounts to 8/9 and is achieved for both β B ω B , and ω B /ω A → ∞.

Conclusions
We devised a consistent description of ergotropic heat engines for the optimal work extraction from a couple of quantum systems which are cyclically restored to the canonical equilibrium at two different temperatures.We provided an exhaustive study for the case of two qutrits with equally-spaced energy levels by deriving the optimal ergotropic transformations, the statistics of the extracted work and the mean entropy production.We focused on the relative fluctuations of the work extracted by each ergotropic transformation and showed that one of them, the double swap U 3 , violates many common TURs, specifically the ones based on the assumption that the distributions of the extracted work for the forward and backward process are the same.
The application of our procedure to systems with higher dimension is promising because it will lead to the generalization of the ergotropic transformations found for the qutrit case and will allow to find new transformations which, as shown in this work, may possibly extract more work on average with lower fluctuations with respect to the two-stroke Otto engines based on qudits.
where Z X = Tr[e −β X H X ].We perform a two-strokes cyclic heat engine by i) isolating the two quantum systems from the reservoir at t = 0 + ; ii) extracting work by a unitary transformation U acting on the two systems up to time t = t; iii) reconnecting the two quantum systems to their respective reservoirs by weak coupling to achieve complete thermalization at t = t ′ ≫ t.We remark that the unitary U incorporates the free evolution of the two systems and their interaction obtained by external (possibly time-dependent) driving protocols, with the condition of being cyclic, namely such that initial and final Hamiltonian coincide, i.e.H X = H X (0) = H X ( t) for both A and B.
The average value ⟨W⟩ of the work extracted on a cycle corresponds to the opposite of the variation of the internal energy of A and B during the unitary stroke U, i.e.

⟨W⟩ = −⟨∆E
Let us now describe the above thermodynamical cycle by a set of stochastic trajectories which correctly reproduce the mean values of all thermodynamic variables by an average of stochastic variables over all possible trajectories.We adopt an operational approach based on complete energy measurements at different times [2,16,23] as in the typical derivation of Jarzynski equality [12].This approach will allow us to study the complete statistics of work extraction and heat exchanges, and in particular to evaluate the fluctuations of work and their relation with the entropy production.As depicted in Fig. A1, we identify a single stochastic trajectory by the outcomes of the sequential fine-grained energy measurements of H A , H B , H R A , H R B at the beginning of the cycle t = 0 + ; of H A , H B at the end of the work stroke t = t operated by U; and finally of H A , H B , H R A , H R B at the end of the thermalization stroke t = t ′ .We denote by V X , with X = A, B, the unitary operator representing the joint evolution of system X and reservoir R X by weak coupling and energy-preserving interaction up to complete thermal equilibrium at t = t ′ ≫ t.By respecting the order of the above measurements, we denote by γ the stochastic trajectory corresponding to the sequence of outcomes, namely γ = {n, m, u, v; l, s; n ′ , m ′ , u ′ , v ′ } (A3) for arbitrary density matrix σ of system X.This fact can be used to simplify Eq. (A11) by summing on all reservoir indexes {u, u ′ , v, v ′ }.Hence, p(W, Q H , Q C ) can be rewritten in terms of measurements outcomes only on systems A and B, namely Eq. (A13) allows one to study the complete statistics of an ergotropic heat engine.The first principle of thermodynamics for the cycle is recovered since the average variation of the internal energy ⟨∆U⟩ = ⟨Q H + Q C − W⟩ correctly gives zero, as shown as follows One also has ⟨Q H ⟩ = −⟨∆E A ⟩, where ⟨∆E A ⟩ is the average variation of the internal energy of A during the unitary stroke, whose expectation can be obtained by averaging E A l − E A n over all trajectories.Similarly, ⟨Q C ⟩ = −⟨∆E B ⟩, where ⟨∆E B ⟩ has corresponding stochastic values given by E B s − E B m .By introducing the trajectories for the thermalization stroke, we remark that the present result allows us to refine the approach of Refs.[33,35,36], where the stochastic values of Q H were identified with −∆E A (and analogously for Q C with −∆E B ).In fact, notice that the relation W = Q H + Q C does not generally hold at the trajectory level.Anyway, since ⟨Q H + ∆E A ⟩ corresponds to the average of E A n ′ − E A n over the trajectories, for increasing number of cycles the discrepancy between the stochastic variables Q H and −∆E A remains bounded by the finite energy of system A, whereas both Q H and ∆E A increase linearly with the number of cycles, thus providing ⟨Q H ⟩ + ⟨∆E A ⟩ = 0. Analogous point applies for the variables Q C and ∆E B .
All results about the stochastic efficiency η of the heat engines given in Refs.[35,36] rigorously hold for η defined in terms of ∆E A , namely η = −W/∆E A .One could refine and compare the results for η defined as η = W/Q H by means of the probability distribution presented here in Eq. (A13).The subtle difference between these two definitions of stochastic efficiency was already discussed in Ref. [45], where the thermalization in a two-stroke Otto engine was modeled by the quantum jump method.
The probability of work extraction is the marginal of p(W, Q H , Q C ) in Eq. (A13) with respect to the heat exchanges, and one has Let us now consider a backward protocol where the measurements on the quantum systems and the reservoirs are performed in the reverse ordering, along with the time-reversal evolution of all interactions.The initial state for the backward protocol is again taken as the product of canonical density matrices, namely as in Eq. (A1).By assuming that all Hamiltonians are invariant under time-reversal at all times [17,18,23], the backward protocol is then equivalent to follow Fig. A1 from the right to the left, along with the replacement of V A , V B and U with V † A , V † B and U † , respectively.We can compare the forward and the backward protocols by the probability P[γ] of a trajectory γ = {n, m, u, v; l, s; n ′ , m ′ , u ′ , v ′ } and the probability P B [γ B ] for the occurrence of the specular reverse trajectory with the same measurement outcomes, namely γ B = {n ′ , m ′ , u ′ , v ′ ; l, s; n, m, u, v}.The comparison between the forward and the backward

Figure 2 .
Figure 2. Scheme of the transformations realizing the permutations P E and P ρ in the different regimes identified by ω in the former case and by βω in the latter.We show these regimes by fixing ω A and the three inequivalent cases for the temperature ratio β A /β B and studying P E and P ρ for increasing ω B .As ω B increases, we find that both the permutations are given by the identity (black thin line), U 2 (purple thick line), U 3 (blue dotdashed line), and the swap U 1 (red dashed line).

Figure 11 .
Figure 11.Ergotropy ⟨W⟩ as a function of ω B in the limiting case β A /β B = 0, with ω A = 100, β A = 0, β B = 10.Blue dot-dashed line: double swap U 3 .Red dashed line: standard swap U 1 inside the corresponding ergotropic regime.Purple solid line: idle swap Ũ2 .Black dotted line: standard swap for any ω B such that the extracted work is positive.

Figure 12 .
Figure 12.Distribution p 1 (W = k(ω A − ω B )) of the work extracted by the standard swap U 1 in ω A − ω B units.We set β A ω A = 0.5 and β B ω B = 2.

Figure 13 . 2 ∑
Figure 13.Distribution p 2 (W = k(2ω B − ω A )) of the work extracted by the idle swap U 2 in 2ω B − ω A units.We set β A ω A = 0.5 and β B ω B = 2.Similarly, in the case of Ũ2 one has

Figure 14 .
Figure 14.Distribution p 3 (W = kω A )) of the work extracted by the double swap U 3 in ω A units in the case ω B /ω A = 3/4.We set β A ω A = 0.5 and β B ω B = 2.

10 − 3 .
The plot shows the violations due to the standard swap U 1 (red dashed line), the idle swap U 2 (solid purple line) and the double swap U 3 (blue dot-dashed line) as a function of β B ω B .

Figure 17 .
Figure 17.SNR of the work extracted by the standard swap U 1 (solid line) in ideal optimal conditions, with β A ω A = 10 −3 .The dotted and dot-dashed lines display the upper bound from the standard TUR in Eq. (2) and the tight TUR in Eq. (50), respectively.The dashed horizontal line highlights the asymptotic limit of the SNR.

Figure 18 .
Figure 18.SNR of the work extracted by the idle swap U 2 (solid lines).The dotted lines display the upper bound from the standard TUR in Eq. (2), while the dot-dashed lines display the upper bound from the tight TUR in Eq. (50).The dashed horizontal lines highlight the limit of the SNR.Up panel: conditions for the maximum SNR independently from the ergotropy extraction, namely β B ω B > β A ω A ∼ 0.Here we set β A ω A = 10 −3 .Bottom panel: conditions for the maximum SNR within the ergotropic regime of U 2 , namely β A ω A > β B ω B ∼ 0.Here we set β B ω B = 10 −3 .

Figure 19 .
Figure 19.SNR of the work extracted by the double swap U 3 .The dotted and dot-dashed lines display the upper bound from the standard TUR in Eq. (2) and the tight TUR in Eq. (50).The dashed brown lines display the bound from the loosest TUR for time-symmetric processes in Eq. (51).The dashed horizontal lines highlight the limit of the SNR.Up panel: conditions for the maximum SNR independently from the ergotropy extraction, namely β B ω B > β A ω A ∼ 0 and ω B /ω A ∼ 0.Here we set β A ω A = ω B /ω A = 10 −3 .The solid blue line displays the SNR.Bottom panel: conditions for the maximum SNR within the ergotropic regime of U 3 , namely β B ω B > β A ω A ∼ 0 and ω B /ω A ∈ [1/2, 1).Here we set β A ω A = 10 −3 and show the cases obtained from two different choices of the frequency ratio.The red line displays the choice optimizing the SNR, i.e. ω B /ω A = 1/2, which reduces the statistics of the work extracted by the double swap to the one extracted by the standard swap.The blue solid line displays the case ω B /ω A = 3/4.

Figure 20 .
Figure 20.Ratio between the squared sum of the mean works extracted in the forward (⟨W 3 ⟩) and backward (⟨ W3 ⟩) processes and the sum of the corresponding variances as a function of β B ω B (solid lines).We set β A ω A = 10 −3 .We display the cases ω B /ω A = 10 −2 (dark-blue thick line) and ω B /ω A = 10 2 (light-blue thin line).The dot-dashed brown and the dashed green curve represent the upper bounds given by Eqs.(52) and Eq.(53), respectively.The dashed black horizontal line identifies the asymptotic value, which amounts to 8/9 and is achieved for both β B ω B , and ω B /ω A → ∞.

)
During the thermalization stroke each system comes back to equilibrium, namely system A absorbs the average heat⟨Q H ⟩ = Tr[H A ρ 0 ] − Tr[H A Uρ 0 U † ] = −⟨∆E A ⟩ from the hot reservoir, and system B dumps ⟨Q C ⟩ = Tr[H B ρ 0 ] − Tr[H B Uρ 0 U † ] = −⟨∆E B ⟩ on the cold one.In our convention the cycle operates as a heat engine when ⟨W⟩ > 0, ⟨Q H ⟩ > 0, and ⟨Q C ⟩ < 0. Clearly, the first law is obtained as ⟨W⟩ = ⟨Q H ⟩ + ⟨Q C ⟩.At each cycle the two quantum systems come back to their respective equilibrium states, and hence the average entropy production per cycle simply corresponds to ⟨Σ⟩

Figure A1 .
Figure A1.Scheme for the stochastic description of a cycle of the Otto two-stroke engine.Labels A and B identify the two systems operated by the engine as a working fluid, R A and R B are the corresponding reservoirs.The unitary U is the transformation extracting work, while V A and V B are the energy-preserving unitaries for the thermal relaxation of each system with its pertaining reservoir.