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.

The probability distributions of a set of thermodynamic variables $\left\{X_i\right\}$ are related to the entropy production $\Sigma$ through the so-called fluctuation theorems, which in general can be expressed as [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28] where $p_B$ refers to the backward process, i.e. to the time-reversed process identified by p [for a self-contained derivation of Eq. (1) see Appendix A and Eq. (A22)]. Through this relation, the symmetries of the processes set relevant constraints on the statistics of the variables $\left\{X_i\right\}$. Another class of relations that connects the statistical properties of mesoscopic and nanoscopic systems to the entropy production is given by the so-called thermodynamic uncertainty relations (TURs) [9,26,29,30,31,32,33,34,35,36,37,38,39]. It has been shown that there is a strong connection between fluctuation theorems and TURs, i.e. every fluctuation theorem implies a specific TUR [32]. Note that the converse does not hold: it was recently found in Ref. [38] a TUR which does not stem from any fluctuation theorem.

Thermodynamic engines which admit a straightforward quantum description are the ones based on the Otto cycle [28,35,36,40,41,42,43,44,45,46,47,48] since the work and heat exchanged are unambiguously identified by their respective distinct strokes. The case considered in this paper, namely a two-stroke Otto cycle, is outlined in Figure 1, where the working fluid is represented by two qutrits.

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 ${\omega }_A$ and ${\omega }_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 $U_3=U_2{U}_1$. Each transformation extracts the ergotropy from a different regime defined by the frequency gaps ${\omega }_A$ and ${\omega }_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 $\mathrm{var}\left(W\right)/{\langle W\rangle }^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] 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 $p_B(\left\{X_i\right\},\Sigma )=p(\left\{X_i\right\},\Sigma )$.

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.

In this work we fix the Planck and Boltzmann constants to natural units, i.e. $\hslash =k_B=1$. We consider two qudits A and B in a product of Gibbs states, i.e. where $H_X={\omega }_X{\sum }_{n=0}^{d_X-1}n|n\rangle \langle n|$ is the Hamiltonian of the system $X=A,B$, each one with equally-spaced energy levels, and $Z_X=\mathrm{Tr}\left[e^{-{\beta }_X{H}_X}\right]$ denotes the corresponding partition function, and ${\beta }_X=T_X^{-1}$ the inverse temperature. The number states $|n\rangle$ in the expansion of the Hamiltonians are eigenstates of the occupation number $n_X\equiv H_X/{\omega }_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 Figure 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 $\langle Q_H\rangle =-\langle \mathsf{\Delta }{E}_A\rangle$. Similarly, for the cold reservoir, $\langle Q_C\rangle =-\langle \mathsf{\Delta }{E}_B\rangle$, and, for the first law of thermodynamics, the average work is given by $\langle W\rangle =\langle Q_H\rangle +\langle Q_C\rangle =-\langle \mathsf{\Delta }{E}_A\rangle -\langle \mathsf{\Delta }{E}_B\rangle$. Correspondingly, the average entropy production reads $\langle \Sigma \rangle =-{\beta }_A\langle Q_H\rangle -{\beta }_B\langle Q_C\rangle =({\beta }_A-{\beta }_B)\langle \mathsf{\Delta }{E}_A\rangle -{\beta }_B\langle W\rangle$. 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 ${\rho }_0$ into a state $\rho =U{\rho }_0{U}^{†}$ such that the average extracted work is maximized, i.e. where here $H=H_A\otimes {\mathbb{I}}_B+{\mathbb{I}}_A\otimes 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 $\mathrm{Tr}\left[\rho H\right]$. In the present case, where the initial state ${\rho }_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{\beta }_A{\omega }_A+m{\beta }_B{\omega }_B)}/\left(Z_A{Z}_B\right)$ 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_{\rho }$ 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 $\tilde{\mathit{E}}={\left\{{\tilde{E}}_l\right\}}_{l=0}^{d_A{d}_B-1}$ satisfies ${\tilde{E}}_l<{\tilde{E}}_{l+1}\forall l\in [0,d_A{d}_B-1)$. Similarly, the permutation $P_{\rho }$ rearranges the occupation numbers of the initial state in descending order, namely, and $\tilde{\mathit{r}}={\left\{{\tilde{r}}_l\right\}}_{l=0}^{d_A{d}_B-1}$ is such that ${\tilde{r}}_{l+1}<{\tilde{r}}_l\forall l\in [0,d_A{d}_B-1)$. 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\rangle \langle 10|$ and $|01\rangle \langle 01|$ are ${\omega }_A$ and ${\omega }_B$, respectively, while the related occupation fractions are $Z_A^{-1}{Z}_B^{-1}{e}^{-{\beta }_A{\omega }_A}$ and $Z_A^{-1}{Z}_B^{-1}{e}^{-{\beta }_B{\omega }_B}$. If we have ${\omega }_A>{\omega }_B$ and ${\beta }_A{\omega }_A<{\beta }_B{\omega }_B$ or the symmetric case where both the order relations are reversed, the transformation that swaps $|10\rangle \langle 10|$ with $|01\rangle \langle 01|$, namely $U=U^{†}=|00\rangle \langle 00|+|11\rangle \langle 11|+|01\rangle \langle 10|+|01\rangle \langle 10|$, extracts the ergotropy. This result appears immediately if we consider the permutation matrices $P_E$ and $P_{\rho }$, which in this case read where $\theta \left(x\right)$ is the Heaviside function. The operator $P_E^{†}{P}_{\rho }$ 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 equally-spaced two-qudit engine is described for any dimension of the qudits by a first partial order over the frequencies $\omega$ and a second one over the products $\beta \omega$. 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 $\mathsf{\Omega }\equiv \{{\omega }_A>{\omega }_B\wedge {\beta }_A{\omega }_A<{\beta }_B{\omega }_B\}$ or $\overline{\mathsf{\Omega }}$, 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{\omega }_A>j{\omega }_B$ and if $k{\beta }_A{\omega }_A>j{\beta }_B{\omega }_B$ for every pair of natural numbers $k\in [0,d_A)$ and $j\in [0,d_B)$.

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=\omega$ or $\beta \omega$ and $X_1\ne 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 $\mathsf{\Omega }$ and $\overline{\mathsf{\Omega }}$.