# Quantum machines powered by correlated baths

Gabriele de Chiara, Mauro Antezza

► To cite this version:

Gabriele de Chiara, Mauro Antezza. Quantum machines powered by correlated baths. Physical Review Special Topics: Physics Education Research, 2020, 2, pp.033315. 10.1103/PhysRevResearch.2.033315. hal-02923218

HAL Id: hal-02923218

https://hal.science/hal-02923218

Submitted on 14 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantum machines powered by correlated baths

Gabriele De Chiara<sup>1,2</sup> and Mauro Antezza<sup>2,3</sup>

<sup>1</sup>Centre for Theoretical Atomic, Molecular and Optical Physics, Queen's University Belfast, Belfast BT7 1NN, United Kingdom

<sup>2</sup>Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France

<sup>3</sup>Institut Universitaire de France, 1 rue Descartes, F-75231 Paris Cedex 05, France

<figure><img src="image_2.png" /></figure>

(Received 26 June 2020; accepted 7 August 2020; published 26 August 2020)

We consider thermal machines powered by locally equilibrium reservoirs that share classical or quantum correlations. The reservoirs are modeled by the so-called collisional model or repeated interactions model. In our framework, two reservoir particles, initially prepared in a thermal state, are correlated through a unitary transformation and afterward interact locally with the two quantum subsystems which form the working fluid. For a particular class of unitaries, we show how the transformation applied to the reservoir particles affects the amount of heat transferred and the work produced. We then compute the distribution of heat and work when the unitary is chosen randomly, proving that the total swap transformation is the optimal one. Finally, we analyze the performance of the machines in terms of classical and quantum correlations established among the microscopic constituents of the machine.

DOI: 10.1103/PhysRevResearch.2.033315

I. INTRODUCTION

The interest in thermal machines powered by quantum working media has recently surged thanks to the technological advancement in the realization and control of individual quanta [1–6]. This tremendous progress has led to the first realizations of quantum engines and thermal devices [7–13].

The theoretical modeling of such devices usually involves the system in contact with equilibrium uncorrelated baths at different temperatures. However, some papers have generalized this picture to nonequilibrium reservoirs [14–21], including the case of the Otto engine in contact with squeezed reservoirs [22–25], which can lead to efficiencies and performances beyond the Otto and Carnot limit. This conclusion, obviously, does not take into account the cost of maintaining a nonequilibrium reservoir which is then considered as a free resource but shows how to best employ these resources (see also Ref. [26]). Other works have considered thermal devices coupled to spatially separated reservoirs which share correlations, classical or quantum [27–33].

Here we propose a general framework based on collisional models [19,23,34–58] which allows us to analyze, in a consistent thermodynamic sense as recently proven in Ref. [59], the effect of classical and quantum correlations between reservoirs in the functioning of quantum thermal machines.

Our setup, depicted in Fig. 1, consists of a working medium made of two quantum systems $S_1$ and $S_2$. Each of these is in contact with a reservoir modeled by the repeated interaction of flying auxiliary particles. These particles are first prepared in

a thermal state at $T_1$ and $T_2$, respectively, and then undergo a unitary operation $U$, which correlates them before their collision with the systems $S_1$ and $S_2$. We study the steady state of the system after many collisions with the flying particles.

Such a microscopic model has the advantage that all thermodynamic contributions, e.g., energy, heat, work and entropy, are accountable and that it is consistent with the laws of thermodynamics [42,43,59]. From the point of view of open quantum systems, our model leads to a nontrivial evolution of a quantum system in contact with correlated reservoirs. As we will show, under the assumption of continuous evolution in the limit of collisions lasting an infinitesimal amount of time, the system’s evolution can be cast in the form of a Markovian Lindblad master equation with collective jump operators acting nontrivially on both systems (see also Ref. [60]).

In this paper, we showcase the functioning of our model assuming the system and environment’s particles to be qubits. However our framework is general and could be equally applied to higher-dimensional systems including infinite-dimensional ones, for instance, quantum harmonic oscillators.

After introducing the setup in more detail in Sec. II and discussing its thermodynamics in Sec. III, we assume in Sec. IV the unitary operation $U$ to be a partial swap. We will show how the amount of swapping between the flying qubits controls the amount and direction of the heat flow among the system’s qubits. We then consider in Sec. V the most general two-qubit unitary operations and study the distribution of work produced and heat exchanged by the system when the unitaries are chosen randomly. We find that the extremal points of the distribution correspond to eight noncorrelating unitaries with the optimal one corresponding to the complete two-qubit swap. In this context, we analyze the quantum and classical correlations established among the quantum constituents of our setup. We find that, while correlations among the flying qubits are not necessary for the machine to work, quantum and classical correlations among the system’s particle must be

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
<figure><img src="image_3.png" /><figcaption>FIG. 1. Setup of our scheme: thermal particles emerging from the cold and hot baths and prepared in the states $\tilde{\rho}_{th}(n_1)$ and $\tilde{\rho}_{th}(n_2)$, respectively, are made to collide to each other under a unitary operation $U$. The emerging correlated particles in the global state $\rho'_B$ collide with the system's particles $S_1$ and $S_2$.</figcaption></figure>

nonzero for the optimal performance. Finally, in Sec. VI, we summarize and conclude.

## II. SETUP AND PRELIMINARIES ON REPEATED INTERACTIONS

We assume the system to be composed of two coupled qubits described by the XXZ Hamiltonian:

$$H_S = J(\sigma_{x1}\sigma_{x2} + \sigma_{y1}\sigma_{y2} + \Delta\sigma_{z1}\sigma_{z2}) + B_1\sigma_{z1} + B_2\sigma_{z2}, \quad (1)$$

where $J$, $\Delta$, and $B_i$ are the interaction strength, anisotropy and local magnetic field, respectively. Here and throughout the paper, we assume $\hbar = k_B = J = 1$ expressing all physical quantities in these units. The operators $\sigma_{xi}$, $\sigma_{yi}$, $\sigma_{zi}$ are the Pauli matrices for the qubit $i$. Notice that the total magnetization $S_z = \sigma_{z1} + \sigma_{z2}$ is a conserved quantity as it commutes with the system Hamiltonian.

We assume the system to be affected by a reservoir modeled by the so-called repeated interactions [19,23,34–58]. In this model, represented in Fig. 1, each qubit of the system interacts with a stream of uncoupled environment qubits (or flying qubits). The interaction between the system qubit and a flying qubit only lasts for a short time $\tau$ during which the interaction Hamiltonian is constant and given by

$$H_{SB} = \sum_{i=1,2} \sqrt{\frac{\gamma(2n_i + 1)}{2\tau}} (\sigma_{xi}\tilde{\sigma}_{xi} + \sigma_{yi}\tilde{\sigma}_{yi}), \quad (2)$$

where the operators $\tilde{\sigma}_{xi}$, $\tilde{\sigma}_{yi}$, and $\tilde{\sigma}_{zi}$ are the Pauli matrices for a flying qubit interacting with the system qubit $i = 1, 2$. The coefficient $\gamma$ determines the strength of the interaction while $n_i = (e^{2B_i\beta_i} - 1)^{-1}$ models the thermal occupation of the flying qubit and is related to the inverse temperature $\beta_i = 1/T_i$ of each bath. Additionally each flying qubit is subject to

the local Hamiltonian:

$$H_{Bi} = B_i \tilde{\sigma}_{zi}, \quad i = 1, 2, \quad (3)$$

and we define $H_B = H_{B1} + H_{B2}$.

Usually, in the literature, the state of the environment qubits for different reservoirs has been assumed to be uncorrelated. In this paper, we challenge this assumption and introduce some correlations, not necessarily quantum ones, between the flying qubits. These are initially uncorrelated and prepared in a thermal state corresponding to a thermal occupation $n_i$:

$$\rho_B = \tilde{\rho}_{th}(n_1) \otimes \tilde{\rho}_{th}(n_2), \quad (4)$$

where

$$\tilde{\rho}_{th}(n_i) = \frac{1}{2}[1 - (1 + 2n_i)^{-1}\tilde{\sigma}_{zi}] \quad (5)$$

and $1$ is the identity operator.

We then correlate the flying qubits with a unitary transformation $U$ so that their state becomes

$$\rho'_B = U \rho_B U^\dagger. \quad (6)$$

After this initial preparation, each flying qubit collides with a system's qubit (see Fig. 1). The collision, lasting for a time $\tau$, is described by the unitary operator:

$$U_{\text{collision}} = e^{-iH_{\text{tot}}\tau}, \quad (7)$$

where $H_{tot} = H_S + H_B + H_{SB}$ is the total Hamiltonian.

If we call $\rho_S(t)$ the state of the system at time $t$, then its state at time $t + \tau$ after the collision becomes

$$\rho_S(t + \tau) = \text{Tr}_B[U_{\text{collision}} \rho_S(t) \otimes \rho'_B U_{\text{collision}}^\dagger]. \quad (8)$$

In the following, we will always consider the steady state that the system approaches after many collisions which is defined by the relation:

$$\rho_S^{\text{steady}}(t + \tau) = \rho_S^{\text{steady}}(t), \quad (9)$$

and for simplicity, we will drop its time dependence and write simply $\rho_S^{\text{steady}}$.

## III. THERMODYNAMICS

We now discuss the different thermodynamics contributions arising in our setup after the system has reached its steady state. As the state of the system does not change anymore, the internal energy variation is zero:

$$\Delta E = \text{Tr}[H_S(\rho'_{SB} - \rho_{SB})] = 0, \quad (10)$$

where $\rho_{SB} = \rho_S^{\text{steady}} \otimes \rho'_B$ and $\rho'_{SB} = U_{\text{collision}} \rho_{SB} U_{\text{collision}}^\dagger$.

To define heat and work we distinguish two scenarios, called *partial* and *complete* scenarios, which we explain in more details in the following.

### A. Partial scenario

In the partial scenario, we assume that we are provided with the flying qubits in the state $\rho'_B$ and that we are not paying for the work of the unitary $U$. In this scenario, we are provided with a nonequilibrium reservoir and thus the Clausius formulation of the second law of thermodynamics may not apply. This is still an interesting scenario to study
for two reasons: first, it continues the investigation of the functioning of thermodynamic machines in the presence of nonthermal environments [14–20,22–24]; second, it gives us an opportunity to study the open quantum system dynamics in the presence of correlated reservoirs as developed in Sec. IV.

Under these assumptions the work is potentially produced or injected during the system-environment collision [59]. This is given by

$$W_{\text{partial}} = \text{Tr}[(H_S + H_B)(\rho'_{SB} - \rho_{SB})]. \quad (11)$$

Similarly the heat exchanged by the system with the flying qubits is equal to the energy balance of the latter ones:

$$Q_{\text{partial}}^{(i)} = -\text{Tr}[H_{Bi}(\rho'_{SB} - \rho_{SB})], \quad i = 1, 2. \quad (12)$$

At steady state, the first law for the two-qubit system reads

$$\Delta E = Q_{\text{partial}} + W_{\text{partial}} = 0, \quad (13)$$

where $Q_{\text{partial}} = Q_{\text{partial}}^{(1)} + Q_{\text{partial}}^{(2)}$. Notice that throughout this paper we employ the convention that positive work or heat corresponds to energy injected into the system contributing to the increase of the system energy.

An entropic formulation of the second law can be derived in terms of the non-negativity of the entropy production (see, for example, Refs. [61–63] and Chapter 28 by R. Uzdin in Ref. [1]):

$$\Sigma_{\text{partial}} = I(\rho'_{SB}) + D(\rho''_{B}||\rho'_{B}) \ge 0. \quad (14)$$

In the expression of $\Sigma$, we have used the mutual information between two quantum objects $O_1$ and $O_2$ defined as

$$I_{O_1 O_2} = S(\rho_{O_1}) + S(\rho_{O_2}) - S(\rho_{O_1 O_2}), \quad (15)$$

where $S(\rho) = -\text{Tr}\rho \ln \rho$ is the von Neumann entropy and $\rho_i$, $i = O_1, O_2, O_1 O_2$ are the density matrices of the corresponding objects. We have also used the relative entropy: $D(\rho||\sigma) = \text{Tr}\rho \ln \rho - \text{Tr}\rho \ln \sigma$ and we have defined the state of the environment after the collision: $\rho''_{B} = \text{Tr}_{S}\rho'_{SB}$. The non-negativity of $\Sigma_{\text{partial}}$ is ensured by the non-negativity of the mutual information and the relative entropy.

## B. Complete scenario

In the complete scenario, we account for the extra work needed to implement the correlating unitary $U$:

$$W_U = \text{Tr}[H_B(\rho'_B - \rho_B)], \quad (16)$$

so that the total work in the complete scenario is the sum of the two contributions:

$$W_{\text{complete}} = W_{\text{partial}} + W_U. \quad (17)$$

Equally, the heat exchanged is the energy balance of the environment during the whole process: from their preparation into the product thermal state $\rho_B$ to their final state after the two unitaries $U$ and $U_{\text{collision}}$:

$$Q_{\text{complete}}^{(i)} = -\text{Tr}[H_{Bi}(\rho'_{SB} - \rho_S^{\text{steady}} \otimes \rho_B)], \quad i = 1, 2, \quad (18)$$

and $Q_{\text{complete}} = Q_{\text{complete}}^{(1)} + Q_{\text{complete}}^{(2)}$ which in general differs from $Q_{\text{partial}}$.

A modified first law holds also in this scenario:

$$\Delta E = Q_{\text{complete}} + W_{\text{complete}} = 0. \quad (19)$$

It is possible to define the entropy production also in the complete scenario:

$$\Sigma_{\text{complete}} = I(\rho'_{SB}) + D(\rho''_{B}||\rho_{B}) \ge 0. \quad (20)$$

Moreover in this case, since $\rho_B$ is the tensor product of equilibrium states for each reservoir, see Eq. (4), we can connect the change in entropy in the system $\Delta S$ with the entropy production and the heat flow:

$$\Delta S = \Sigma_{\text{complete}} + \sum_{i=1}^{2} \beta_i Q_{\text{complete}}^{(i)} \quad (21)$$

from which we can write the Clausius inequality:

$$\Delta S - \sum_{i=1}^{2} \beta_i Q_{\text{complete}}^{(i)} = \Sigma_{\text{complete}} \ge 0, \quad (22)$$

where at steady state $\Delta S = 0$.

## IV. RESERVOIRS IN A PARTIALLY SWAPPED LOCALLY THERMAL STATE

Here, we consider the special case in which the unitary $U$ correlating the flying qubits is a partial swap operation:

$$S_{\phi} = \exp \left\{ -i \frac{\phi}{2} (\tilde{\sigma}_{x1} \tilde{\sigma}_{y2} - \tilde{\sigma}_{y1} \tilde{\sigma}_{x2}) \right\}, \quad (23)$$

which is a total swap for $\phi = \pi/2$. This leads to the following expression for the bath density matrix after the action of the partial swap $S_{\phi}$:

$$\rho'_{B} = \frac{1}{(1 + 2n_1)(1 + 2n_2)} \quad (24)$$

$$\times \begin{pmatrix} n_1 n_2 & 0 & 0 & 0 \\ 0 & \frac{1}{2}(n_1 + n_2 + 2n_1 n_2 + (n_1 - n_2) \cos 2\phi) & (n_2 - n_1) \sin \phi \cos \phi & 0 \\ 0 & (n_2 - n_1) \sin \phi \cos \phi & \frac{1}{2}(n_1 + n_2 + 2n_1 n_2 - (n_1 - n_2) \cos 2\phi) & 0 \\ 0 & 0 & 0 & (1 + n_1)(1 + n_2) \end{pmatrix}, \quad (25)$$

where the basis of eigenstates of $\sigma_{z1}$ and $\sigma_{z2}$ has been used to write the matrix representation of $\rho'_B$.

The action of the partial swap is to partially exchange populations between the auxiliary qubits. Indeed, the two flying qubits are still locally in a thermal state although at a
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Effective populations of the two flying qubits as a function of $\phi$ for $n_1 = 1$ and $n_2 = 3$. (b) Mutual information between the two flying qubits as a function of $n_1$ and $n_2$ for $\phi = 0.2\pi$. The solid line corresponds to zero entanglement and divides separable and entangled states.</figcaption></figure>

different temperature compared to their state before the partial swap:

$$ \tilde{\rho}'_i = \mathrm{Tr}_i \rho'_B = \tilde{\rho}_{\mathrm{th}}(N_i) \quad (26) $$

where $i = 1, 2$ and $\mathrm{Tr}_i$ represents the partial trace with respect to the qubit other than $i$. The effective population after the application of $S_\phi$ is

$$ N_i = \frac{1 n_1 + n_2 + 4n_1n_2 + (-1)^i(n_2 - n_1) \cos 2\phi}{1 + n_1 + n_2 - (-1)^i(n_2 - n_1) \cos 2\phi}. \quad (27) $$

The effective populations $N_i$ of the two flying qubits are shown in Fig. 2 where it is evident that for $\phi = \pi/2$ the populations are completely swapped.

The partial swap creates classical and quantum correlations, measured for instance by the mutual information $I_{A1A2}$ of $\rho'_B$ among the two auxiliary qubits. The mutual information is always nonzero as shown in Fig. 2 except for the special values of $\phi = m\pi/2$, $m \in \mathbb{Z}$, regardless of the values $n_1$ and $n_2$. This happens because, as it is evident from Eq. (24), $\rho'_B$ is not a product state. Notice how, in Fig. 2, unbalancing the two populations $n_1$ and $n_2$ leads to larger correlations.

The work for implementing the partial swap is given by

$$ W_U = - \frac{2(B_1 - B_2)(n_1 - n_2) \sin^2 \phi}{(1 + 2n_1)(1 + 2n_2)}. \quad (28) $$

Similarly to what is found in Refs. [59,64], the work vanishes for $B_1 = B_2$ or $n_1 = n_2$ and is maximum for the total swap $\phi = \pi/2$.

We set the parameters in such a way that for $\phi = 0$ the two-qubit system acts as an engine ($W < 0$, $Q^{(1)} < 0$, $Q^{(2)} > 0$ for both scenarios). It is possible to find the analytical expression of the steady state, the work and heats exchanged per cycle. While their expressions are quite long and hard to read, it is possible to show that all their ratios are proportional to ratios of the local magnetic fields, $B_1$ and $B_2$. As a consequence, and similarly to other models [30,59,64], the efficiency $\eta = |W|/Q^{(2)}$, if the setup works as an engine, or the coefficient of performance (COP) $\eta_{COP} = Q^{(1)}/W$, if the setup works as a refrigerator, correspond to the Otto values. For example, in the former case, the efficiency is given by ($B_2 < B_1$)

$$ \eta = 1 - \frac{B_2}{B_1} \quad (29) $$

and in the latter case, the COP is given by ($B_1 < B_2$)

$$ \eta_{\mathrm{COP}} = \frac{B_1}{B_2 - B_1}. \quad (30) $$

In the partial scenario, $\phi$ can be used to control the functioning of the thermodynamic machine making it switch from an engine to a refrigerator. In the interval $0 \le \phi \le \pi/2$, the machine behaves as an engine for $0 \le \phi < \pi/4$ and as a refrigerator for $\pi/4 < \phi \le \pi/2$. For values of $\phi$ outside the interval $[0, \pi/2]$, the situation just described repeats periodically. The switching points are found for $N_1 = N_2$ which occurs for $\cos 2\phi = 0$, i.e., $\phi = (2m + 1)\pi/4$, $m \in \mathbb{Z}$. These switching points correspond to effective Carnot points where work and heat exchanged are zero (see also Ref. [59] where the condition $n_1 = n_2$ corresponds to the Carnot point).

Numerical results for the partial scenario are plotted in Fig. 3 which shows the periodic change of the machine

<figure><img src="image_5.png" /><figcaption>FIG. 3. Thermodynamic quantities: work and heat exchanged per collision for the steady state of the system in contact with the flying qubits prepared with the partial swap Eq. (23). (Top) Partial scenario. (Bottom) Complete scenario. Parameters: $\Delta = 1$, $\gamma = 1$, $B_1 = 0.1$, $B_2 = 0.3$, $n_1 = 0.1$, $n_2 = 2$, and $\tau = 0.1$.</figcaption></figure>
operating as an engine and as a refrigerator. Interestingly, the largest absolute values of the thermodynamic quantities are obtained for the values $\phi = m\pi/2$, $m \in \mathbb{Z}$ at which the partial swap operation $S_{\phi}$ corresponds to the identity (even $m$) and the total swap (odd $m$).

Thus the partial swap operation can be employed as a valve to control the direction of the heat flowing between the system's qubits without altering the equilibrium reservoirs that prepare the flying qubits. These can be particularly useful in physical implementations in which one does not have full control of the environment.

Let us now move to the complete scenario. In this case, since we are taking into account all work contributions including those coming from the action of the swap operation, both first and second laws, the latter expressed in terms of heat exchanged, are fulfilled. The corresponding results plotted in

Fig. 3 show that the work and heats exchanged do not change sign and the machine always behaves as an engine. The action of the partial swap is to amplify the work production and the heats exchanged, while keeping the same performance, reaching a maximum for $\phi = m\pi/2, m = \pm 1, 3, 5, \dots$ at which the swap is total.

This shows that given the two flying qubits, initially in equilibrium states with fixed temperatures, the maximum work that can be extracted is achieved by completely swapping their states before making them to collide with the system's qubit. As we will see in Sec. V, this is the maximum value that can be obtained for any unitary operation between the flying qubits.

We end up this section by considering the continuous limit of $\tau \to 0$. Up to first order in $\phi$, one obtains the following Lindblad master equation:

$$ \dot{\rho}_S = -i[H_S, \rho_S] + \gamma \sum_{i=1,2} (1+n_i) \mathcal{L}_{\sigma_{-i}}(\rho_S) + n_i \mathcal{L}_{\sigma_{+i}}(\rho_S) + \frac{\gamma \phi(n_2-n_1)[\mathcal{M}(\sigma_{+1}, \sigma_{-2}, \rho_S) + \mathcal{M}(\sigma_{-2}, \sigma_{+1}, \rho_S) + \text{H.c.}]}{\sqrt{(1+2n_1)(1+2n_2)}}, \quad (31) $$

where $\mathcal{L}_a(\rho) = 2a\rho a^\dagger - a^\dagger a\rho - \rho a^\dagger a$ is the usual Lindblad operator and $\mathcal{M}(a, b, \rho) = 2a\rho b - ba\rho - \rho ba$ is a modified Lindblad-like operator. In Eq. (31), we have also defined the jump operators $\sigma_{\pm i} = \frac{1}{2}(\sigma_{xi} \pm i\sigma_{yi})$.

The presence of quantum correlations in the initial state of the bath is the reason for the appearance in the master equation of the collective term $\mathcal{M}$ proportional to $\phi$ which corresponds to environment-induced processes of emission of an excitation from one qubit and absorption from the other.

In this continuous limit, it is possible to define work power and heat currents both in the partial and complete scenario. However, in the complete scenario, one should assume that the work necessary for the swap operation given in Eq. (28) scales to zero as $\tau \to 0$. This is indeed the case up to first order in $\phi$ consistently with our master equation.

## V. RESERVOIRS PREPARED BY RANDOM UNITARIES

Here we generalize the approach developed in the previous section by considering thermal flying qubits which are subject, before the collision with the system qubits, to a unitary transformation $U_R$. The unitary transformation is always the same for all collisions for a given setup and allows the system to reach a steady state, which we then analyze. We repeat the same procedure for an ensemble of $6 \times 10^6$ random unitaries $U_R$ chosen with uniform distribution according to the Haar measure [65].

For a generic $U_R$, the creation of finite quantum coherences and correlations between the environment qubit may lead to difficulties in the derivation of a consistent master equation in the continuous limit $\tau \to 0$, see Ref. [50]. In this section, we will therefore restrict our analysis to a small but finite $\tau$. Under this assumption, the repeated interactions remain a discrete map for the system which, after many applications reaches a steady state, which we analyze.

We consider a setup which, for $U_R = 1$, corresponding to the situation in which no unitary is applied to the flying qubits before their collision with the system, operates as an engine. Analogous results can be obtained for the refrigerator regime.

In the partial case, as in Sec. IV, work and heat fluxes are proportional to each other with the prefactor being a ratio of the local magnetic fields, $B_1$ and $B_2$, related to the Otto efficiency. For this reason we only show the analysis of the work probability distribution shown in Fig. 4. Its probability density function (PDF) is approximately Gaussian with a very small average statistically compatible with zero.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Histogram of the partial work $W_{\text{partial}}$ with the vertical axis measuring the corresponding probability density function (PDF) for $6 \times 10^6$ random unitaries. The solid line is the best fit normal distribution with average $9.6 \times 10^{-7}$ and standard deviation $1.1 \times 10^{-3}$. The vertical dashed lines signpost the minimum and maximum possible values of $W_{\text{partial}}$ obtained with noncorrelating unitaries (see text). Parameters as in Fig. 3 except for the local magnetic fields which are $B_1 = 0.1$ and $B_2 = 0.15$.</figcaption></figure>
The maximum and minimum values of the partial work $W_{\text{partial}}$ are obtained with unitaries that do not create correlations, classical or quantum, among the auxiliary qubits. As we will see in more detail in the complete scenario, this type of unitary operations do not create any mutual information between the flying qubits but simply rearrange the populations of the four basis states. In the partial scenario, we find that the minimum negative value of the work (largest value of produced work) is obtained with the unitary that inverts the populations of the qubit prepared in the hottest temperature and leaves unchanged the qubit prepared in the coldest temperature: $U_R = \sigma_{x2}$ (corresponding to operation III defined later). In contrast, the maximum positive value of the work (largest value of the work injected), is obtained with the unitary that inverts the populations of the qubit prepared in the coldest temperature and leaves unchanged the qubit prepared in the hottest temperature: $U_R = \sigma_{x1}$ (corresponding to operation VI defined later).

Let us now pass to the complete scenario in which the work produced or extracted is not necessarily proportional to the two heat fluxes. The results show that, by applying different unitary operations, the different modes (accelerator, heater or heat engine) can be achieved with different efficiencies. We plot in Fig. 5 the joint histogram of the complete work and the heat input from the hot environment. The distribution is confined by a nonregular octagon whose eight vertexes correspond to unitary operations that do not create quantum or classical correlations between the environment qubits. Let $\{p_1, p_2, p_3, p_4\}$ be the vector of the populations of the density matrix $\rho_B$ of the flying qubits after their preparation equilibrium states but before undergoing the unitary $U$. These eight unitary operations only affect the populations of the flying qubits according to where $S_{\pi/2}$ is the transformation

<table><thead><tr><th>Label</th><th>Populations</th><th>Unitary</th></tr></thead><tbody><tr><td>I</td><td>$\{p_1, p_2, p_3, p_4\}$</td><td>$1_1 1_2$</td></tr><tr><td>II</td><td>$\{p_1, p_3, p_2, p_4\}$</td><td>$S_{\pi/2}$</td></tr><tr><td>III</td><td>$\{p_2, p_1, p_4, p_3\}$</td><td>$1_1 \sigma_{x2}$</td></tr><tr><td>IV</td><td>$\{p_2, p_4, p_1, p_3\}$</td><td>$S_{\pi/2} 1_1 \sigma_{x2}$</td></tr><tr><td>V</td><td>$\{p_3, p_1, p_4, p_2\}$</td><td>$S_{\pi/2} \sigma_{x1} 1_2$</td></tr><tr><td>VI</td><td>$\{p_3, p_4, p_1, p_2\}$</td><td>$\sigma_{x1} 1_2$</td></tr><tr><td>VII</td><td>$\{p_4, p_2, p_3, p_1\}$</td><td>$S_{\pi/2} \sigma_{x1} \sigma_{x2}$</td></tr><tr><td>VIII</td><td>$\{p_4, p_3, p_2, p_1\}$</td><td>$\sigma_{x1} \sigma_{x2}$</td></tr></tbody></table>

that completely swap the state of the two qubits, see Eq. (23) for $\phi = \pi/2$.

Notice the symmetry between the unitary operations and the position of the corresponding vertex of the octagon. In Fig. 5, two vertexes are opposite on the octagon if their transformations can be obtained one from the other by applying the total swap operation. Moreover, one of the sides of the octagon, delimited by the vertexes I and II, corresponds to the Otto efficiency given in Eq. (29), so that points above the line corresponds to engines operating at a lower efficiency or, for $W_{\text{complete}} \ge 0$ operating as an accelerator, if $Q_{\text{complete}}^{(2)} > 0$ or a heater if $Q_{\text{complete}}^{(2)} < 0$. In the complete scenario, there are no unitaries that lead to refrigeration, similarly to what we discussed in Sec. IV for the partial swap.

<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Joint histogram of the complete work $W_{\text{complete}}$ versus heat $Q_{\text{complete}}^{(2)}$. The dashed gray line represents the Otto efficiency. The blue dots, connected by thin solid lines, are obtained for noncorrelating unitaries marked by the corresponding number. These special points form an extremal octagon within which all other machines obtained for random unitaries must be located. The color code identifies the number of values in each bin of the histogram. Parameters: $\Delta = 1$, $\gamma = 1$, $B_1 = 0.1$, $B_2 = 0.15$, $n_1 = 0.1$, $n_2 = 2$, and $\tau = 0.1$. (b) Extremal octagons obtained for the same parameters of the top panel but different values of $B_2 = 0.15, 0.3, 0.6$, and $0.9$, corresponding to vertexes denoted by circles, squares, triangles, and diamonds, respectively.</figcaption></figure>

The optimal point II yielding the larger amount of work produced (minimum negative value) corresponds to the complete swap operation. This can be related to an engine proposed by Campisi, Pekola and Fazio also based on a complete swap transformation but applied to the qubits of the working fluid [64].

In Fig. 5, we also show how the shape and size of the extremal octagon change when $B_2$ is varied while all other parameters are kept fixed. Notice that, as we increase $B_2$, the region corresponding to the engine increases, leading to a higher probability of achieving this operation mode. This can be understood as follows. First, $Q_{\text{complete}}^{(2)} \propto B_2$, thus if we rescaled the horizontal axis of the lower diagram in Fig. 5, all the vertexes would have the same horizontal coordinates. Second, for $W_{\text{complete}} = W_{\text{partial}} + W_U$ the situation is more involved. While the partial work always fulfils $W_{\text{partial}} \propto B_2 - B_1$, the work $W_U$ needed to implement the noncorrelating unitaries is a linear function of $B_2$ and $B_1$ which depends on the actual transformation. Therefore there exists no rescaling of the vertical axis that would bring the vertexes of different
<figure><img src="image_8.png" /><figcaption>FIG. 6. Histograms of the mutual information and quantum discord versus the total work $W_{complete}$. From left to right: (a) $I_{S1S2}$, (b) $I_{A1A2}$, (c) $I_{AS}$, and (d) $D_{S1S2}$ (see definition in the main text). The blue dots are obtained for noncorrelating unitaries specified by their label I-VIII.</figcaption></figure>

octagons to the same vertical coordinate. Furthermore, the vertex IV, corresponding to an accelerator for $B_2 = 0.15$, turns into an engine for the larger values of $B_2$ we analysed.

To get more insight into the relationship between the functioning of the whole system as an engine and the quantum features of the working medium quantum steady state, we looked at two measures of correlation between two parts of the engine. The first is the mutual information defined in Eq. (15). The mutual information quantifies both classical and quantum correlations. If one subtracts the maximum amount of classical correlations that can be obtained by local measurements, the quantum discord, a genuine measure of quantum correlations, is obtained [66,67]. This can be defined as:

$$D_{O1O2} = I_{O1O2} - J_{O1O2}, \qquad (32)$$

where the classical information $J_{O1O2}$ is the maximum information that can be extracted on $O_2$ if we perform local measurements on $O_1$:

$$J_{O1O2} = S(\rho_{O2}) - \min_{\{\Pi_i\}} \sum_{i=1}^{N} q_i S(\tilde{\rho}_i). \qquad (33)$$

In Eq. (33), we have defined the probabilities $q_i = \text{Tr}[\Pi_i \rho_{O1O2} \Pi_i]$ of the outcome $i$ and the post-measurement states $\tilde{\rho}_i = \text{Tr}_{O1}[\Pi_i \rho_{O1O2} \Pi_i]$ of the object $O_2$. The minimization is done over all possible sets of measurements $\{\Pi_i\}$ on $O_1$, not necessarily orthogonal projectors.

The results are shown in Fig. 6. We start with the distribution of the mutual information $I_{S1S2}$ between the system qubits in the steady state. This shows that, although the optimal point II corresponds to uncorrelated flying qubits, the system qubits are nevertheless correlated as a result of their direct interaction and of reaching a nonequilibrium steady state due to the multiple collisions with the environment. A similar distribution is also obtained for the quantum discord $D_{S1S2}$, which shows that the system qubits are also genuinely quantum correlated.

We have also considered the mutual information $I_{A1A2}$ between the flying qubits after the unitary but before the collision with the system. This shows that the preparation unitaries applied before the correlations do indeed create a lot of correlations, quantum or classical, but these do not necessarily lead to large values of work produced or injected. The extremal operation II in fact corresponds to zero mutual information between the flying qubits, which are thus in a product state. Similar conclusions, although quantitatively different, are reached when analyzing concurrence and discord between the two flying qubits.

Finally, we have analysed the mutual information $I_{AS}$ between the flying qubits and the system after the collision. As before, the state of the system is steady. Thus, although it does not change during the collision, it sustains correlations to be created between system and each pair of environmental qubits. These correlations are necessary for allowing exchange of heat between the two reservoirs through the system. The distribution of $I_{AS}$ plotted in Fig. 6, shows that to achieve large amounts of work, and consequently heat exchange, it is sufficient a small value of $I_{AS}$. In particular the noncorrelating operations IVIII, that leave the flying qubits in a product state, correspond to the smallest values of $I_{AS}$. Finally, in all these figures of merit, the noncorrelating operations IVIII do not necessarily correspond to extremal points since these functions, mutual information and discord, are not linear functions of the state in contrast to the average work and heat.

## VI. CONCLUSIONS

In this paper, we have proposed a general framework to model quantum thermal machines in contact with correlated reservoirs using repeated interactions. Different conclusions are found depending on whether we assume the partial or the complete scenario, the latter one being always consistent with the laws of thermodynamics.
We have shown how, in the partial scenario, the amount of partial swapping among the flying qubits, allows one to control the operating mode of the thermal machine switching it from an engine to a refrigerator.

In the case of random unitaries, we found a complex geometrical structure in the distribution of heat and work bounded by a nonregular octagon whose vertexes correspond to noncorrelating unitaries. We analysed the discord and mutual information between the flying qubits, the system qubits and between system and environment. We found that correlations in the system steady state and between system and environment are necessary to achieve the optimal performance.

Our work leads the way to future studies of open quantum systems with correlated environments consistently with thermodynamics.

All the numerical data presented in this work can be found in Ref. [68].

ACKNOWLEDGMENTS

We thank Marco Cattaneo, Francesco Ciccarello, Adam Hewgill, and Gabriel Landi for useful discussions. G.D.C. thanks the CNRS and the group Theory of Light-Matter and Quantum Phenomena of the Laboratoire Charles Coulomb for hospitality during his stay in Montpellier. We acknowledge support from the UK EPSRC EP/S02994X/1.

[1] Thermodynamics in the Quantum Regime, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer, Cham Switzerland, 2018).

[2] J. Millen and A. Xuereb, New J. Phys. 18, 011002 (2016).

[3] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, J. Phys. A: Math. Theor. 49, 143001 (2016).

[4] S. Vinjanampathy and J. Anders, Contemp. Phys. 57, 545 (2016).

[5] M. T. Mitchison, Contemp. Phys. 60, 164 (2019).

[6] S. Deffner and S. Campbell, Quantum Thermodynamics: An Introduction to the Thermodynamics of Quantum Information (Morgan & Claypool Publishers, 2019).

[7] J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, Science 352, 325 (2016).

[8] A. Ronzani, B. Karimi, J. Senior, Y.-C. Chang, J. T. Peltonen, C. Chen, and J. P. Pekola, Nat. Phys. 14, 991 (2018).

[9] G. Maslennikov, S. Ding, R. Hablützel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Nat. Commun. 10, 202 (2019).

[10] D. von Lindenfels, O. Gräb, C. T. Schmiegelow, V. Kaushal, J. Schulz, M. T. Mitchison, J. Goold, F. Schmidt-Kaler, and U. G. Poschinger, Phys. Rev. Lett. 123, 080602 (2019).

[11] J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl, K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walmsley, R. Uzdin, and E. Poem, Phys. Rev. Lett. 122, 110601 (2019).

[12] J. P. S. Peterson, T. B. Batalhão, M. Herrera, A. M. Souza, R. S. Sarthour, I. S. Oliveira, and R. M. Serra, Phys. Rev. Lett. 123, 240601 (2019).

[13] M. Gluza, J. Sabino, N. H. Y. Ng, G. Vitagliano, M. Pezzutto, Y. Omar, I. Mazets, M. Huber, J. Schmiedmayer, and J. Eisert, arXiv:2006.01177.

[14] B. Leggio, B. Bellomo, and M. Antezza, Phys. Rev. A 91, 012117 (2015).

[15] B. Leggio and M. Antezza, Phys. Rev. E 93, 022122 (2016).

[16] B. Reid, S. Pigeon, M. Antezza, and G. De Chiara, EPL (Europhysics Letters) 120, 60006 (2017).

[17] R. J. de Assis, T. M. de Mendonça, C. J. Villas-Boas, A. M. de Souza, R. S. Sarthour, I. S. Oliveira, and N. G. de Almeida, Phys. Rev. Lett. 122, 240602 (2019).

[18] C. Cherubim, F. Brito, and S. Deffner, Entropy 21, 545 (2019).

[19] M. Pezzutto, M. Paternostro, and Y. Omar, Quantum Sci. Technol. 4, 025002 (2019).

[20] F. Carollo, F. M. Gambetta, K. Brandner, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. Lett. 124, 170602 (2020).

[21] R. Román-Ancheyta, B. Çakmak, and Ö. E. Müstecaplıoğlu, Quantum Sci. Technol. 5, 015003 (2019).

[22] J. Roßnagel, O. Abah, F. Schmidt-Kaler, K. Singer, and E. Lutz, Phys. Rev. Lett. 112, 030602 (2014).

[23] G. Manzano, F. Galve, R. Zambrini, and J. M. R. Parrondo, Phys. Rev. E 93, 052120 (2016).

[24] B. K. Agarwalla, J.-H. Jiang, and D. Segal, Phys. Rev. B 96, 104304 (2017).

[25] V. Singh and Ö. E. Müstecaplıoğlu, arXiv:2006.08311.

[26] W. Niedenzu, V. Mukherjee, A. Ghosh, A. G. Kofman, and G. Kurizki, Nat. Commun. 9, 165 (2018).

[27] P. Doyeux, B. Leggio, R. Messina, and M. Antezza, Phys. Rev. E 93, 022134 (2016).

[28] D. Türkpençe, F. Altintas, M. Paternostro, and Ö. E. Müstecaplıoğlu, EPL (Europhysics Letters) 117, 50002 (2017).

[29] B. Karimi and J. P. Pekola, Phys. Rev. B 96, 115408 (2017).

[30] A. Hewgill, A. Ferraro, and G. De Chiara, Phys. Rev. A 98, 042102 (2018).

[31] C. L. Latune, I. Sinayskiy, and F. Petruccione, Quantum Sci. Technol. 4, 025005 (2019).

[32] G. Manzano, G.-L. Giorgi, R. Fazio, and R. Zambrini, New J. Phys. 21, 123026 (2019).

[33] O. Pusuluk and Ö. E. Müstecaplıoğlu, arXiv:2006.03186.

[34] V. Scarani, M. Ziman, P. Štelmachovič, N. Gisin, and V. Bužek, Phys. Rev. Lett. 88, 097905 (2002).

[35] M. Ziman, P. Štelmachovič, V. Bužek, M. Hillery, V. Scarani, and N. Gisin, Phys. Rev. A 65, 042105 (2002).

[36] D. Karevski and T. Platini, Phys. Rev. Lett. 102, 207207 (2009).

[37] V. Giovannetti and G. M. Palma, Phys. Rev. Lett. 108, 040401 (2012).

[38] F. Ciccarello, G. M. Palma, and V. Giovannetti, Phys. Rev. A 87, 040103(R) (2013).

[39] S. Lorenzo, A. Farace, F. Ciccarello, G. M. Palma, and V. Giovannetti, Phys. Rev. A 91, 022121 (2015).

[40] G. T. Landi, E. Novais, M. J. de Oliveira, and D. Karevski, Phys. Rev. E 90, 042142 (2014).

[41] S. Lorenzo, R. McCloskey, F. Ciccarello, M. Paternostro, and G. M. Palma, Phys. Rev. Lett. 115, 120403 (2015).

[42] F. Barra, Sci. Rep. 5, 14873 (2015).
[43] P. Strasberg, G. Schaller, T. Brandes, and M. Esposito, Phys. Rev. X **7**, 021003 (2017).

[44] F. Ciccarello, *Quantum Measurements and Quantum Metrology* **4**, 53 (2017).

[45] M. Pezzutto, M. Paternostro, and Y. Omar, *New J. Phys.* **18**, 123018 (2016).

[46] S. Cusumano, V. Cavina, M. Keck, A. De Pasquale, and V. Giovannetti, Phys. Rev. A **98**, 032119 (2018).

[47] S. Campbell, F. Ciccarello, G. M. Palma, and B. Vacchini, Phys. Rev. A **98**, 012142 (2018).

[48] J. A. Gross, C. M. Caves, G. J. Milburn, and J. Combes, *Quantum Sci. Technol.* **3**, 024005 (2018).

[49] M. H. Mohammady, H. Choi, M. E. Trusheim, A. Bayat, D. Englund, and Y. Omar, Phys. Rev. A **97**, 042124 (2018).

[50] F. L. S. Rodrigues, G. De Chiara, M. Paternostro, and G. T. Landi, Phys. Rev. Lett. **123**, 140601 (2019).

[51] A. Manatuly, W. Niedenzu, R. Román-Ancheyta, B. Çakmak, Ö. E. Müstecaplıoğlu, and G. Kurizki, Phys. Rev. E **99**, 042145 (2019).

[52] P. Strasberg, Phys. Rev. Lett. **123**, 180604 (2019).

[53] S. Seah, S. Nimmrichter, and V. Scarani, Phys. Rev. E **99**, 042103 (2019).

[54] B. Çakmak, S. Campbell, B. Vacchini, Ö. E. Müstecaplıoğlu, and M. Paternostro, Phys. Rev. A **99**, 012319 (2019).

[55] G. Guarnieri, D. Morrone, B. Çakmak, F. Plastina, and S. Campbell, Phys. Lett. A **384**, 126576 (2020).

[56] G. García-Pérez, M. A. Rossi, and S. Maniscalco, npj Quantum Inf. **6**, 1 (2020).

[57] X.-M. Li, Y.-X. Chen, Y.-J. Xia, Q. Zhang, and Z.-X. Man, Chin. Phys. B **29**, 060302 (2020).

[58] D. Cilluffo, A. Carollo, S. Lorenzo, J. A. Gross, G. M. Palma, and F. Ciccarello, arXiv:2006.08631.

[59] G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, and M. Antezza, New J. Phys. **20**, 113024 (2018).

[60] S. Daryanoosh, B. Q. Baragiola, T. Guff, and A. Gilchrist, Phys. Rev. A **98**, 062104 (2018).

[61] M. Esposito, K. Lindenberg, and C. V. den Broeck, New J. Phys. **12**, 013013 (2010).

[62] D. Reeb and M. M. Wolf, New J. Phys. **16**, 103011 (2014).

[63] K. Ptaszyński and M. Esposito, Phys. Rev. Lett. **123**, 200603 (2019).

[64] M. Campisi, J. Pekola, and R. Fazio, New J. Phys. **17**, 035012 (2015).

[65] M. Ozols (2009).

[66] H. Ollivier and W. H. Zurek, Phys. Rev. Lett. **88**, 017901 (2001).

[67] L. Henderson and V. Vedral, J. Phys. A: Math. Gen. **34**, 6899 (2001).

[68] G. De Chiara and M. Antezza (2020), doi:10.17034/46387195-ff10-4b6c-8725-a1fed5bad759.