Defining and generating current in open quantum systems

Defining current in open quantum systems can be problematic: No general description exists for the current of operators not conserved by the system-environment interaction. We fill this gap by deriving a general formula for probability current on an arbitrary graph, universally applicable to any system-environment interaction. We furthermore provide a representation of the average current, whereby the operator is first measured weakly, then strongly. When the dynamics is of Lindblad form, we derive an explicit formula for the current. We exemplify our theory by analysing a simple Smoluchowski-Feynman-type ratchet, operating deep in the quantum regime. Consisting of only two interacting particles, each moving on a three-site ring, the ratchet displays several novel quantum effects, such as tunnelling-induced current inversion, which we relate to the onset of quantum contextuality, and steady-state entanglement generation in the presence of arbitrarily hot environment. The role of spatial symmetry in current generation is also studied.


I. INTRODUCTION
Transforming thermal, "disordered" motion into an "ordered form" is the purpose of heat engines and thermal motors [1]. In order to operate a heat engine, one only needs to create non-equilibrium conditions, e.g., in the form of two thermal baths at differing temperatures. Putting a thermal motor in motion, however, requires some kind of symmetry breaking in addition to thermal disequilibrium [2][3][4][5]. In the classical regime, thermal motors are rather well-understood [2][3][4][5][6][7][8][9][10][11][12]. For instance, autonomous motors propelled by temperature gradient and broken spatial symmetry have been investigated in the context of Brownian motors [9,10,12], with a design similar in spirit to the original Feynman ratchet [2].
On the other hand, despite their fundamental and technological importance, only now quantum thermal motors are receiving due attention in the literature [13][14][15][16][17][18]. Motivated by [12], in this work, we introduce a minimalistic model of an autonomous quantum rotor, consisting of two particles with 3-dimensional Hilbert spaces. Although we will work in the approximation of weak coupling to the baths, our analysis starts from a microscopic Hamiltonian of the global -rotor plus baths -system, in order to properly account for the symmetries (and their breaking). The model under consideration is a twoparticle 3-component Potts model (q = 3) [19], which is a "higher-spin" (spin-1, in our case) generalization of the Ising model, complemented with quantum tunnelling. Alternatively, one can think of two atoms on an optical lattice [20], each of which is confined to 3 lattice sites.
Here we restrict to the 3-position model as it is the minimal case exhibiting the symmetry breaking necessary for ordered motion to occur. Importantly, our motor does not require three-body interactions [13,15,17] in order to operate.
Despite the considerable attention quantum walks and quantum transport [21][22][23][24][25][26][27][28][29] in open quantum systems have received, the problem at hand, i.e., defining local particle current in a strongly interacting system undergoing a global dissipative evolution, has never been studied in the literature before. We fill this gap by deriving a surprisingly simple, yet universal, formula for particle current on a general quantum network, applicable well beyond the scope of the present paper.
Using standard techniques from the theory of open quantum systems, we fully characterize the nonequilibrium steady states of the rotor. This allows us to make a rigorous connection between the symmetries of the rotor's Hamiltonian and the main transport properties of the system: particle current and heat flux. Moreover, we find that an interesting current inversion phenomenon takes place: despite the fact that tunnelling in our model is rotationally symmetric and thus does not favour any specific direction of particle propagation (see Fig. 1 and Eq. (3)), the particle current can change its direction as the intensity of tunnelling is varied.
In addition to generating steady current, our machine can perform other useful tasks. Enabled by symmetry and fueled by local coherence, our machine is capable of converting uncorrelated states of the rotor into entangled steady states, for arbitrarily high temperatures of the baths. Moreover, in the presence of sufficiently strong tunnelling, powered by temperature difference, the machine can "charge" the initially "empty" rotor with extractable work, where by empty we mean that no work The rotor consists of two partitions a and b, each consisting of three locations a particle (represented by a green circle) can occupy. The particles can tunnel between sites at rate τ . The interaction between the particles is illustrated by the connecting green line. The blue and red semicircles illustrate the fact that each partition is coupled to its own thermal bath.
can be extracted from the state.

II. ROTOR
Our model rotor consists of two subsystems, a and b, each living in a 3-dimensional Hilbert space. These can be thought of as two distinguishable spin-1 particles. The self-Hamiltonians of these particles are taken to be zero in order to ensure translation invariance and periodicitythe state space of the particle is but three identical levels of the same energy. We then add a classical interaction potential of Potts form [19] between the particles, which corresponds to two spins pointing to one of the 3 equispaced directions specified by the angles 2πj/3. Such a potential mimics a dipole-dipole interaction. It can be immediately seen that the potential is periodic in each coordinate in the sense that U (j a + 3, j b ) = U (j a , j b + 3) = U (j a , j b ). However, unless φ = 0, the potential breaks the exchange symmetry in that swapping the particles bears an energetic cost: U ja,j b = U j b ,ja . The breaking of this type of symmetry is responsible for the occurrence of transport, as detailed in Sec. V.
The quantum analogue of this discrete stochastic classical model is described by the Hamiltonian which, in the basis of vectors |j a , j b , has matrix elements where |j is the state of the particle being at the position j. Let us now adopt the view where the quantum states are locations in space (see Fig. 1 for an illustration) to which the particle is confined. One can think of, e.g., an atom in an optical lattice [20]. Since the confining potential cannot be infinite, tunnelling between the sites will be present. We add tunnelling to the Hamiltonian in the form of hopping terms: where I is the identity operator in the corresponding Hilbert space and Here, the index is cyclic in the sense that |k + 3 ≡ |k .
In the first order in τ , such a Hamiltonian describes tunnelling subject to the following kinematic constraints: (i) there can be no tunnelling between a and b and (ii) the particles jump only one by one and never simultaneously. Note that simultaneous jumps are present in the higher orders in τ .

III. INTERACTION WITH THE ENVIRONMENT
Let us now take the next step and add dissipation to the picture. More specifically, we are interested in a setup where the spin a is coupled to a bath at temperature T a and spin b is connected to a bath at temperature T b . In the limit of weak interaction with the envronment, the standard way of describing the evolution of a quantum system is by means of Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation (ME). There are two standard approaches towards prescribing a GKSL ME to an open quantum system's dynamics. In one approach, one derives an equation for the state of the system from the exact dynamics of the system-plus-environment composite under the Born-Markov and secular approximations [30]. In the case of a multipartite system, the equation takes into account the "global" properties of the system's Hamiltonian even if the environment acts locally on the subsytems. The other approach, often referred to as "local" GKSL, ascribes local deissipators to the subsystems without taking into account the global properties of the system's Hamiltonian. This is an incomplete description and generally the thermal state is not a steadystate solution of the local ME when the temperatures of the baths are equal, which brings about thermodynamically inconsistent behaviour [31][32][33][34][35] such as spontaneous heat flow against the temperature gradient [31] or nonzero heat flow between two thermal baths at the same temperature [34,35]. This is contrasted with the global approach, where the dissipation operates on the global Hamiltonian and accounts for the interaction in a thermodynamically consistent manner. It is important to note, however, that, in the weak coupling regime of some models [36,37], the local ME leads to a steady state closer to that obtained via an exact solution of the system as a whole, including the baths. Moreover, as opposed to the global ME, it recovers the correct classical dynamics (see below).

A. Local master equation
The classical analogue of the system [12] consists of two three-state particles interacting via the potential in Eq. (1). In order for local detailed balance to be welldefined, the particles are not allowed to jump simultaneously. Denoting the probability of particle a being in state j a and particle b in state j b via p ja,j b , we can write the general ME in the form where W ja,j b |j a ,j b are the transition rates satisfying local detailed balance conditions: This picture directly translates into a local GKSL ME in the quantum analogue without tunnelling. Indeed, let us write the following quantum ME: where S[L][ρ] ≡ LρL † − 1 2 {L † L, ρ} (with {·, ·} denoting the anticommutator) and L ja,j a ,j b = |j a , j b j a , j b |, L ja,j b ,j b = |j a , j b j a , j b |. This is the standard quantum ME used in the theory of quantum random walks [27,28]. When the Hamiltonian is classical (i.e., H = H cl ), the steady-state solution of Eq. (7) is a diagonal density matrix with the diagonal being the steady-state solution of the classical ME (given by Eq. (5)). Moreover, if the initial density matrix is also diagonal, the probability vector composed of the its diagonal elements will evolve according to Eq. (5).
As mentioned above, the local GKSL ME may not be thermodynamically consistent. In our case, it becomes problematic when τ = 0. Indeed, when β a = β b = β, the steady-state solution of Eq. (7), ρ loc st is not a thermal state at inverse temperature β. In fact, ρ loc st − 1 Z e −βH ∝ τ (where Z = Tr e −βH ), so the deviation cannot be considered to be small.

B. Global master equation
In order to derive a consistent global GKSL equation, we will start with a microscopic Hamiltonian and use the Born, Markov, rotating-wave, and secular approximations [30]. This is a standard procedure [30], therefore we will only briefly summarize the relevant formulas.
The total system Hamiltonian is standardly chosen to be of the form where A α and B α are operators living, respectively, in the Hilbert spaces of the rotor and the baths, and H a and H b are the Hamiltonians of the baths. As a generic example, one can think of bosonic baths, composed of a large number of harmonic oscillators [30]. The resulting master equation will then be Here, the "jump operators" are given by where E m is an eigenvalue of H and Π Em is the eigenprojector corresponding to it. Note that the ω's are all the gaps of the system's bare Hamiltonian H. Note also that, as a further approximation, we omit the Lamb shift correction to the Hamiltonian in Eq. (9) [30]. Being composed of operators Λ α (ω) and commuting with H [30], the Lamb shift term has the same symmetry properties as H, which makes this omission "safe" (especially since we are working in the weak system-bath coupling regime, and the Lamb shift is a second-order effect [30]). The bath correlation functions, γ α (ω), ("transition rates", in the language of the preceding subsection) are given by γ α (ω) = ∞ −∞ dse iωs B † α (s)B α (0) , and, for a generic bosonic bath living in a single spatial dimension, amount to [30] γ α (ω) = g|ω| 1 − e −βα|ω| × e βαω when ω ≤ 0 1 when ω > 0 .
These functions satisfy the detailed balance condition: γα(ω) = e −βαω , which, combined with Eq. (10), guarantees that, when β a = β b , the thermal state is a steady solution for the global ME given by Eq. (9) [30]. This means that, unlike the local ME, the global ME is thermodynamically consistent. We will use γ α (ω) also when dealing with the classical and local GKSL MEs, which means that the transition rate from, say, j a j b to j a j b , will be given by We emphasize that Eq. (9) is qualitatively different from the local GKSL ME (Eq. (7)) in that the jump operators now operate between the eigenspaces of H, which are not diagonal in the "location" basis. Moreover, since the spectrum of H contains degeneracies and some of the gaps may be repeated, the jump operators (Eq. (10)) will generally not be of rank-1 (like operators L are) and will have non-zero non-diagonal elements. Moreover, the global ME does not coincide with the local ME even when the inter-particle interaction, which is controlled by K, goes to zero. Neither it does in the "classical" limit of τ → 0. Indeed, the jump operators of the global ME, due to their dependence on the spectrum of H (which is degenerate for all the values of the parameters), have a coherence structure which is different from that of the jump operators of the local ME. The discrepancy between the two descriptions is in fact exacerbated when K or τ go to zero, which is due to increase in the degeneracy of H. Given that there is no interaction between the particles when K 1, it is obvious that it is the description provided by the global ME that fails. This is related to the fact that, as K decreases, so do also some of the gaps in the spectrum of H, which, for small enough K's, leads to the breakdown of the secular approximation vital for deriving the global ME [30,36].

IV. CURRENT
Let us now turn to the problem of defining particle current on each individual side of the rotor. In the classical case (Eqs. (5)-(6)), the local particle current, say, on the side a, between neighbouring sites j a and j a (note that, since there are only three states and the system is periodic, all sites are neighbours) will be given by In the steady state, it is straightforward to show that J In the quantum regime, when working with local GKSL ME (Eq. (7)), the current is given by Eq. (12) whenever the Hamiltonian is diagonal (for instance, when it is equal to H cl ). This is due to the fact that the quantum dynamics is exactly the same as the classical stochastic dynamics. However, whenever the global GKSL ME is used or the state is not diagonal (e.g., when it is the steady solution of the dynamics when the system Hamiltonian is non-diagonal), defining the current is less straightforward. Indeed, although both the probabilities to find the particles at given sites and transition rates between these sites are always well-defined [21,26,27,29], defining the current in a state with non-zero coherences via Eq. (12) would be to ignore the coherent part of that state. Instead, in such situations, standard quantum transport theory defines the current through the continuity equation. For example, if one measures the transport in terms of particle velocity, then one considers the continuity equation for the position operator in the Heisenberg picture: dx/dt = −divJ [23-25, 38, 39]. Following that logic, let us introduce the position operator of the particle a at location j a : where I b is the identity operator acting on the Hilbert space of b. Now, in order to maintain full generality and simplify notation, let us write the ME in the general GKSL form Then, in the open-system Heisenberg picture, the position operator will evolve according tȯ where the dot over the object denotes the time derivative.
First of all, we notice immediately that the Hamiltonian part can be represented as where (17) is the tunnelling-induced current. This object is essentially the discretized version of the standard current associated with Schrödinger equation: 2mi Ψ † ∇Ψ − ∇Ψ † Ψ (where Ψ is the state in the "position" representation), and the associated current operator is 1 2m {p, |x x|} [40] (cf. Eq. (26)). A similar operator (17) is also used to describe the tunnelling spin current in spin-1 2 lattice systems [25,38,41].
Let us turn to the dissipative part of the RHS of Eq. (15), Keeping in mind the trivial identities Λ λ = Λ λ I and Λ † λ = IΛ † λ I, we use the identity resolution which takes into account the periodicity condition j + 3 ≡ j and holds for any j a , to represent D * [x ja ] as a sum of terms of the form x j a Λ † λ x j a Λ λ x j a . Now, as can be checked by direct inspection, all these terms can be rearranged in such a way that we can write where Here, the symbol (j a , j a , j a ) is introduced for ease of notation and stands for x j a Λ † λ x j a Λ λ x j a . Using the identity resolution (18) again, we can further simplify the above expression into Moreover, reading from Eqs. (15) and (17), it is a simple exercise to show that the total current operator, ja→ja+1 , can be written as Note that this expression is in the Heisenberg picture, and one can return to the Schrödinger picture by substituting the operator derivatives from Eq. (15) and considering all operators constant. Specifically, the average current at the moment of time t is given by ja→ja+1 are defined in, respectively, Eqs. (17) and (25), with all the operators in the Schrödinger picture.
Note also that the current defined by Eq. (26) trivially satisfies the continuity equation dx ja /dt = −divJ = J ja−1→ja − J ja→ja+1 for any dynamics that respects x ja−1 + x ja + x ja+1 = I for all times. This is a very general condition encompassing all conceivable physically meaningful dynamical maps. Indeed, in the Schrödinger picture, it is equivalent to the dynamical map being trace-preserving. In particular, this means that Eq. (26) describes the current also when the dynamics is non-Markovian. In such cases, however, the dependence of the state at a given time on (generally, all) the preceding states [30] will generically render using Eq. (26) impractical.
The remarkably simple Eq. (26) for the local current of a general compound system undergoing general dissipative evolution under the influence of its own (possibly) strongly interacting Hamiltonian and the environment, is the central result of this section. To the best of our knowledge, this formula has not been reported in the literature before. Moreover, quantum weak measurements [42,43] offer a neat interpretation for Eq. (26). Indeed, let us rewrite J ja→j a at the moment of time t as so that the average current reads where the symbol (j a ↔ j a ) denotes the repeating of the same term as on its left, but with the j a and j a indices interchanged. Here, p j a (t + ) = Tr(ρx j a (t + )) is the probability of detecting the particle at site j a at the moment of time t + , and is the real part of the weak value of x ja (t) on the preselected state ρ and postselected on x j a (t + ) [42][43][44]. More specifically, it can be interpreted as the conditional average of x ja at the moment of time t conditioned on the measurement of x j a at a later moment t + [44]. Since x ja is the projector on the location j a , the conditional average is the same as the conditional probability.
In other words, P x ja (t)|x j a (t+ ) is the probability of finding the particle a on the site j a at the moment of time t, as a result of a weak (minimal-disturbance) measurement of x ja (t) on the system in the state ρ, conditioned on a measurement of the particle at the site j a at a later moment of time t + . In turn, this prompts the interpretation of P x ja (t)|x j a (t+ ) p j a (t+ ) as the joint probability of finding the particle at location x ja at the moment of time t and at location x j a at the moment of time t + . Therefore, the formula for the current (28) fits the classical intuition of "flow forward minus flow backward", encapsulated in Eq. (12). Importantly, the location at t is measured weakly so that the state is disturbed minimally before the second measurement. In this way, the system's state is not lost while we observe it jump from j a to j a . Note that the distribution P x ja (t)|x j a (t+ ) p j a (t + ) = Re Tr x ja (t)x j a (t + )ρ is also known as Terletsky-Margenau-Hill quasiprobability distribution [45,46] (see also [47] for a further discussion on the connection between weak value and Terletsky-Margenau-Hill distribution). We also note that, besides providing an appropriate theoretical language for interpreting Eq. (26), weak values can also be directly accessed experimentally (see the relevant discussion in Ref. [43]). Despite the widespread applications of weak measurements in various areas of quantum physics (see, e.g., [43,[48][49][50][51]), to the best of our knowledge, the above connection with quantum transport has not been previously made in the literature.
In order to appreciate the importance of the first measurement being weak, let us see what the current looks like when both position measurements are non-weak (i.e., are standard quantum von Neumann measurements [30,52]). Let us turn to the Schrödinger picture for convenience, and measure the particle a at site j a at the moment of time t. The probability of finding the particle at site j a that will be given by p ja (t) = Tr(x ja ρ(t)), and, after the measurement, the state of the system will collapse . Within the period of time , the state will evolve into ρ ja + L[ρ ja ]+O( 2 ). Therefore, the probability of finding the particle at site x j a , at the moment of time t + , will be Tr x j a ρ ja + L[ρ ja ] + O( 2 ) . Taking the reverse flow into account, we thus find the two-strong-measurement (TSM) current to be given by Taking into account that x ja x j a = 0 (since j a = j a ), we immediately see that J coincides with the RHS of Eq. (20). Hence, making the first measurement strong, eliminates all coherent contributions to the current, given by Eqs. (17), (21)- (24). Note that the latter include also the tunnelling current.
Interestingly, when the ME is given by the local GKSL equation (see Eq. (7)), it is straightforward to show that the total current (26) is given by the sum of J (tun) ja→j a (given by Eq. (17)) and J (TSM) ja→j a (given by Eq. (30)), i.e., the environment-induced current can be measured by performing two strong measurements of the position. Substituting the local jump operators L ja,j a ,j b and L ja,j b ,j b into Eq. (20), we find that the thermal current, J (th) ja→j a , is given by ja→j a obtained from the local GKSL ME has the current in Eq. (12) as its classical limit. Recalling that the transition to the classical regime is indeed provided by the local GKSL equation, we thus summarize that the current in Eq. (26) reverts to the classical current (given by Eq. (12)) in the appropriate classical limit.
Lastly, let us remark that Eq. (26) is applicable to the most general situation of a particle undergoing an arbitrary trace-preserving evolution on an arbitrary N -vertex graph. Indeed, in such a case, the continuity equation takes the form where x j denotes the projector onto the vertex j. (Note that the divergences in Eqs. (16) and (19) are special cases of Eq. (31).) Now, keeping in mind that N j=1 x j = I, it is straightforward to show that defining the vertexvertex currents J k→j via Eq. (26) turns Eq. (31) into an identity. Furthermore, the interpretation in terms of weak measurements, which we developed above, applies in this most general situation without modifications. And, as above, when the dynamics is Markovian and can be described by a GKSL equation, the current operators J k→j can be directly read off from the master equation.

V. STEADY STATES AND SYMMETRIES
Transport is controlled by symmetry [4,23,[53][54][55]. This is universally observed in the literature, and the reference is made to the symmetries of the generator of the evolution. In this section, we will describe the symmetries of the Hamiltonian and show how they relate to the transport properties of the machine.

A. Classical ME
We are interested in the steady-state rotation of our system, and we will start by analysing the classical setting. There, the evolution is given by Eq. (5), which, upon collecting p ja,j b 's and the transition rates into, respectively, the vector p and matrix W, can be rewritten as dp dt = Wp. The stochastic matrix W, in turn, depends on the potential U through Eq. (6) (Eq. (11) will be our specific choice).
Since the phase space is discrete, all transformations are given by permutation matrices. And a configuration function, e.g., the steady state p, is symmetric under a transformation P if P p = p. Importantly, Eqs. (5)- (6) ensure that U and p have the same symmetries. Note, however, that, although W and U might have common symmetry transformations, due to the nonlinear relation between W and U , not all symmetries of U will in general be respected by W. The global "rigid" rotation symmetry of the potential (namely, U i+k,j+k = U i,j ) provides an example of the first kind: W is also symmetric under the global rotation, which is expressed in the fact that W ja+1 j b +1,j a +1 j b +1 = W jaj b ,j a j b . Interestingly, W maintains a single steady state despite this symmetry. As discussed in the next subsection, this is related to the fact that, although the local GKSL equation giving rise to the classical evolution (see Sec. III A) is, as a whole, symmetric under global rotations, its individual jump operators are not. Following the terminology of Ref. [53], one may say that the classical dynamics is hence only "weakly" symmetric under global rotations. Now, we observe that there is no current in the system whenever φ = kπ/3 (k is an arbitaray integer), and the current is non-zero whenever φ = kπ/3 (see Fig. 2). At the same time, we notice that, whenever φ = kπ/3, U j,i = U i,j+k . In other words, when φ = kπ/3, the potential is symmetric under the combination of particle exchange and a unilateral rotation. The disappearance of the current in this case can be explained by a simple physical argument. Indeed, exchanging the particles is equivalent to swapping the temperatures of the baths. On the other hand, taking, for simplicty, φ = 0, we see that the system and its state remain intact (recall that the steady state inherits all U 's symmetries). Put otherwise, inverting the direction of the temperature gradient does not alter the currents. Now, drawing our intuition from phenomenological non-equilibrium thermodynamics [56], where, for small temperature differences, the current is a linear function of the temperature gradient, we thus conclude that the current must be zero.
Interestingly, the generator of the evolution, W, is not symmetric under the above "generalized" exchange symmetry. Nevertheless, the increased degeneracy in U 's spectrum (U has 2 distinct eigenvalues when φ = kπ/3, and 3 when φ = kπ/3) brings about degeneracies in W, leaving it with only 5 distinct eigenvalues (contrasting W's having no degeneracies when φ = kπ/3), which in turn means that W is more symmetric when φ = kπ/3. B. Quantum ME

Local ME
The quantum and classical cases share many similarities when one considers the local GKSL ME. As discussed above, it is only weakly symmetric under global rotation. Indeed, the latter is given by the unitary operator R = ja,j b |j a + 1, j b + 1 j a , j b | (and its integer powers), and it is easy to see that R does not commute with the jump operators L µ . On the other hand, it can be seen by direct inspection that, denoting the right-hand side of Eq. (7) [53], and our local ME, given by Eq. (7), indeed has a unique steady state for any value of τ . This can be seen as an artifact of the global (system plus baths) Hamiltonian not being symmetric under global rotation.
As we can see in Fig. 2, where the average of J ja→ja+1 in the steady state (the same for all j a 's) is plotted against φ, the current is again zero whenever φ = kπ/3. The generalized exchange symmetry of the potential, discussed in the previous section, is also a symmetry of H. Indeed, the unitary operator corresponding to the generalized exchange is given by Ξ k = SWAP · I ⊗ R k b , where SWAP is the swap operator (for ∀ |ψ and |ξ , SWAP|ψξ = |ξψ ) and R b = 3 k=1 |k + 1 k| rotates the particle b by one step. It can be easily seen that [H, Ξ k ] = 0. Similarly to the classical case, although the generator of the evolution is not symmetric under Ξ k (Ξ k L loc [ρ]Ξ † k = L loc [Ξ k ρΞ † k ]), the steady state is. Moreover, taking into account that, in the steady state, unilateral rotation does not affect the average local current (J j b →j b = J j b +1→j b +1 ), we conclude that the transformation Ξ k is equivalent to swapping the temperatures of the baths. Hence, by the same physical argument as in the classical case, we restore the intuition as to why the current is zero for φ = kπ/3. Importantly, we note that the average steady-state current is 2π/3-periodic. Indeed, the potential satisfies U ja,j b (φ + 2π/3) = U ja,j b +1 (φ), which, along with the fact that the tunnelling parts of the Hamiltonian are local-rotation-invariant, means that a 2π/3 phase shift is equivalent to a unilateral rotation. On the other hand, the average local current in the steady state is invariant under local rotation, which proves the 2π/3-periodicity of the average current. This periodicity is evident in Fig. 2.
An interesting purely quantum effect can be read off from Fig. 2: the current changes its direction as the tunnelling rate increases. This is a purely quantum phenomenon, and it occurs only when T a is sufficiently low. Current inversion, albeit through the fundamentally different mechanism of external periodic driving, was also reported in Brownian ratchet systems [3,57].
As a final remark, let us note that the thermodynamic Local GKSL ME: The average current of particle a, Ja , versus the phase φ in the steady state. The blue curve corresponds to zero tunnelling and is equivalent to the classical case. The orange and green curves correspond to, respectively, τ = 0.2 and τ = 0.6. The other parameters are: where the temperature is lower, the quantum effects of the tunnelling become more dominant, and, with the increase of the tunnelling rate, are expressed by changes in the current's direction. No similar inversions in the direction of the current are observed for particle b, which is attached to a hottertemperature bath. As it can be clearly seen, the current is zero whenever φ = kπ/3. inconsistency of the local GKSL manifests itself when T a = T b . In such a case, we find that, although J a + J b = 0, J a = 0. This is, however, only a second-order effect in terms of the tunnelling rate: J a ∝ τ 2 .

Global ME
In order to properly account for the role of the symmetry breaking in the global GKSL ME, in this section, we will study the global, microscopically derived quantum ME (see Eq. (9)) and require the total Hamiltonian in Eq. (8) to be symmetric under global rotation. Given that H already commutes with R, we thus need to choose A a and A b also commuting with R. As shown in Ref. [53], such "strong" symmetry necessarily implies that the global GKSL ME has multiple steady states. More specifically, there are at least as many linearly independent steady states of the ME as there are distinct eigenvalues of R. Since R has 3 distinct eigenvalues, the steady-state subspace of L glob will be of 3 or more dimensions. Importantly, this ambiguity means that the steady state of the evolution bears some memory on the initial state [53,54]. Let us choose where X is given by Eq. (4). Being the arithmetic mean of the Gell-Mann matrices λ 1 , λ 4 , and λ 6 [58], X indeed generalizes the x-component of spin (the standard 2-dimensional Pauli X operator) from SU (2) to SU (3). With such a choice, we have [H tot , R] = 0, and, for k = kπ/3, our numerical analysis reveals that the steady state depends on two free parameters. Given that the strong symmetry with respect to R guarantees 3 trace-1 solutions to L glob [ρ] = 0 [53], we thus conclude that each eigensubspace of R contains a unique steady state and that there are no non-zero traceless solutions to L glob [ρ] = 0. Let us now write the eigenresolution of R as R = 3 k=1 e 2πi 3 k R k , where the cubic roots of 1 are the eigenvalues and R k are the eigneprojectors (R 2 = R * 1 ). By numerically finding an arbitrary solution for the degenerate linear system L glob [ρ] = 0 -call itρ -we determine the three mutually orthogonal steady states as ρ st k = R kρ R k Tr(R kρ ) , which of course do not depend onρ. An important property of the steady states is that, although they are not invariant under any of the three Ξ k 's, their marginals coincide for all choices of the parameters: Now, if the rotor's initial state is ρ 0 , then the steady state it eventually evolves into, will be Indeed, keeping in mind that k R k = I, we can write Since there are no non-zero traceless solutions to L glob [ρ] = 0 and R k ρ 0 R k (k = k ) are traceless (and remain so because the evolution is trace preserving), the second sum in the decomposition will vanish as the system converges to its steady state. Whereas each of the states R k ρ0R k Tr(R k ρ0) will evolve into its corresponding ρ st k , leaving us with Eq. (34). We can now use Eqs. (17) and (25) to determine the currents associated to all three ρ st k 's. First of all, for all the values of the parameters, all the average currents are independent of the site (e.g., J ja→ja+1 [ρ st ] is the same for all j a 's). Furthermore, we find (numerically) several other general properties that also hold for all values of the parameters. First, the average current in ρ st 3 is zero for both particles a and b: Next,  (35), and (36), and using the numerical values of R 1 and R 2 , for an arbitrary initial state ρ 0 , we find The total average current versus the tunnelling rate τ , with the phase φ = π/6. The total current starts negative, and, depending on the initial state, can turn positive as τ increases. For both plots, the rest of the parameters are: where θ 0 = 3 ja,j b =1 j a , j b |ρ 0 |j a + 1, j b + 1 is the sum of some of the non-diagonal elements of the initial density matrix.
An important feature the global GKSL ME shares with the local GKSL ME, is that, whenever φ = kπ/3, J (th) = 0 for any initial state. In Fig. 3a, we plot the total (thermal plus tunnelling) current of particle a against the phase φ, for different values of τ . We see that, unlike the local ME, the tunnelling current is not zero when φ = kπ/3. Moreover, in stark contrast to the local ME, the average current for nearly zero tunnelling rate is noticeably smaller than the average currents for larger tunnelling rates.
When φ = kπ/3, the Hamiltonian of the rotor is symmetric under Ξ k , whereas the coupling operators, A a and A b (and therefore also H tot ), are not. Nevertheless, φ = kπ/3 is marked by an increased degeneracy in the spectrum of H tot , and as a consequence, of L glob , which is expressed in the fact that the steady state now depends on 5 free parameters, meaning that the steady-state space of L glob is 6-dimensional. The steady states are all symmetric under global rotation, but are not invariant under Ξ k . By noticing that, similarly to Sec. V B 1, due to the invariance of the local average current under local rotation, the Ξ k transformation is equivalent to swapping the temperatures of the baths, for any k. Hence, we can invoke the same physical argument as in the cases of classical and local GKSL MEs in order to interpret the nullification of the currents. We emphasize that, as in the above two cases, the appearance of non-zero current in the system is related to breaking of (some of) the symmetries of the generator of the evolution, taking place whenever φ = kπ/3. Finally, we note that the current is a 2π/3-periodic function of φ. As in Sec. V B 1, this is due to the fact that a 2π/3 shift in the phase is equivalent to a local rotation, and the local average current in the steady state is invariant under local rotation.
Furthermore, the current inversion phenomenon, observed in Sec. V B 1, is also present here. In Fig. 3b, we plot the total average current for all three basis steady states, as a function of τ . The phase is chosen to be φ = π/6. Keeping Eqs. (34) and (37) in mind, we see that, unless the initial state is completely in the subspace of R 3 , J < 0 for τ 1, and if Im θ 0 < 0, J will become positive for large enough τ . The inversion of the direction of the current is counterintuitive in that the tunnelling is symmetric with respect to local rotation and therefore does not favour a particular direction for particle flow.
Lastly, we note that, whenever T a = T b , J (th) α = 0 for any initial state. The tunnelling currents, on the other hand, do not necessarily turn to zero when T a = T b . This however does not contradict the second law of thermodynamics, and is simply due to the fact that, although the thermal state -for which the tunnelling currents are indeed zero -is a steady-state solution for the global ME, the basis steady states, ρ st k , are not thermal.

C. Heat
In the absence of external driving, heat always flows from hot to cold. Due to clear separation between the two baths in the total Hamiltonian (see Eq. (8) Since the system does not exchange energy with external media, in the steady state, the energy conservation requiresQ a +Q b = 0. Moreover, at thermal equilibrium, i.e., when T a = T b , the heat flow must be zero. In our For the global ME, we have numerically established the following properties of the heat flow: This means that, when the tunnelling is not zero, the heat flow can be controlled by changing the weights of the eigensubspaces of R in the initial state of the rotor. Such symmetrycontrolled manipulation is well-known in the literature (see, e.g., Ref. [55]). In Fig. 4, we plotQ[ρ st 1 ] (Q[ρ st 3 ] has the same behaviour) against φ, for various values of τ . First, we notice that the heat flow decreases with increasing the tunnelling rate, which is related to the fact that the higher the τ , the higher the relative weight of the non-interacting component of H. We further notice that the zeroes of J (th) , i.e., φ = kπ/3, correspond to global extrema ofQ. Lastly, analogously to the current, the heat flow is also a 2π/3-periodic function. This fact is evident from Fig. 4, and is proven as follows. As viewed from the perspective of the Hamiltonian, a 2π/3 phase shift is equivalent to a unilateral rotation; specifically, Being the eigenprojectors of H, Π Em undergo the same transformation, which, combined with the fact that [X, (10) and (32)). Hence, it also holds that We conclude this section by emphasizing that caution must be maintained when using the global ME for computing the heat flow for small K. As discussed in Sec. III B, the secular approximation, indispensable for the applicability of the global ME, is compromised when K approaches 0. Given the results in Refs. [36,37] for the Caldeira-Leggett model, it is reasonable to expect that the local ME would provide a more reliable description in that regime. This issue can be decisively settled only upon solving the global, system-plus-baths dynamics in the limit of infinitely large baths, doing which is however unfeasible (see, e.g., [30]).

VI. UTILISING THE STEADY STATE
In this section, we will study several useful properties of the steady states of the machine. Specifically, we will focus on the entanglement and the amount of unitarily extractable work (also known as ergotropy) in the rotor.

A. Entanglement
Entanglement is a key resource in basically all aspects of quantum technology [52,59]. Although thermal noise can sometimes be beneficial for creating and maintaining entanglement (see, e.g., [60][61][62]), at high temperatures, the effect of environment will typically be detrimental. Ineed, the maximally mixed state in any finitedimensional Hilbert space is not entangled, and there exists a finite-volume convex set around the maximally mixed state in which all states are not entangled [59]. In particular, this means that, for any Hamiltonian, there exists a finite temperature T [H], above which all thermal states are non-entangled. On the other hand, thermal disequilibrium is known to be beneficial for entanglement generation [61][62][63][64][65][66][67][68][69]. Therefore, a key question in this regard is whether a thermal disequilibrium created by temperatures higher than T [H] can generate entangled non-equilibrium states. Of especial interest would be the steady states, since these can be prepared and maintained reliably and robustly. In this section, we will show that a thermal initial state can be transformed into an entangled steady state only if it is already entangled. Thus, for thermal initial states the above question has a negative answer. However, we will show that, for initially uncorrelated states, local coherence can be traded for entanglement in the steady state, for any temperatures of the baths. In other words, powered by global rotation symmetry, our machine can convert a local quantum resource (coherence) into a global quantum resource (entanglement), even when the surroundings have such high temperatures that neither resource would survive thermalization.
Our rotor is a 3 × 3-dimensional system, and in 9 (or higher) dimensions is geneally very hard to tell whether a state is entangled or not [59]. In this section, we will use a very strong necessary condition, called positive partial transpose (PPT) criterion [59]. The partial transpose of ρ with respect to partition b, ρ Γ b , is defined as the matrix with entries j a , j b |ρ Γ b |j a , j b = j a , j b |ρ|j a , j b . Now, if ρ is not entangled, then ρ Γ b will be a non-negative matrix. Hence, ρ will necessarily be entangled if ρ Γ b has a negative eigenvalue. This condition is also sufficient Global GKSL ME. The steady-state entanglement between particles a and b, as measured by the negativity, versus the phase φ. Only the negativity for ρ st 1 (which is also equal to that of ρ st 0 ) is plotted since N (ρ st 2 ) < N (ρ st 1 ) has the same behaviour. The (inset) shows the steady-state negativity when the initial state is chosen to be the thermal state at the temperature of the hot bath, as a function of τ , with φ = π/3. The rest of the parameters are the same as those in Fig. 3.
when the joint Hilbert space dimension is ≤ 6, however, in higher dimensions, e.g., for our rotor, this condition is not sufficient: ρ can be entangled even when ρ Γ b ≥ 0 [59]. The entanglement detected by the PPT criterion can be measured by the negativity, N [ρ], defined as the sum of the absolute values of all negative eigenvalues of ρ Γ b [70]: where ||O|| = Tr √ O † O is the trace norm [52]. Similarly to the heat flux (see Sec. V C), we find numerically that (i) N (ρ st 1 ) = N (ρ st 2 ) > N (ρ st 3 ) > 0 and (ii) the global extrema of the negativity correspond to the zeroes of the thermal current (φ = kπ/3). In Fig. 5, we plot N (ρ st 1 ) against φ, for three different values of tunnelling: τ = 10 −3 , τ = 0.2, and τ = 0.4. It shows, in particular, that the negativity is a monotonically decreasing function of τ . This seemingly counterintuitive fact can be easily understood by first observing that the presence of entanglement even when τ → 0 is caused by the fact that, in order to obtain the basis steady state, we project on the highly entangled eigensubspaces of the global rotation matrix R (cf. Sec. V); and the decrease of entanglement with the increase of τ is caused by the increased weight of the non-interacting part of the Hamiltonian H − H cl (see Eq. (3)). The standard intuition is recovered when one considers realistic initial states. We show this in the inset of Fig. 5, where we plot the steady-state negativity versus τ . There, we fix K = 2 and φ = π/3, and, for any τ , choose the initial state to be the thermal state at the temperature of the hot bath: ρ 0 ∝ e −β b H . As expected, we see that there is no entanglement for weak tunnelling. Then, having reached a maximum at an intermediate value of τ , the entanglement decays as H becomes more and more dominated by the non-interacting tunnelling term.
The dependence of the entanglement in the steadystate of a global GKSL ME on the temperature difference of the baths has been extensively studied in qubit systems [61][62][63][64][65][66][67]. In highly asymmetric problems, situations when entanglement can be generated only when there is temperature difference have been reported [62]. In our model, we find that, for any values of all the other parameters fixed, if we set the cold temperature T a constant, then the negativity of all three basis steady states is a monotonically decreasing function of T b . Furthermore, the negativity at T b = T a is a monotonically decreasing function of T a . Interestingly, however, we find that, when both keeping T a constant and varying T b or fixing T b = T a and varying T a , the entanglement does not experience sudden death: the negativity decays gradually, reaching zero only asymptotically. Therefore, by choosing the initial state appropriately, we can guarantee the presence of entanglement in the steady state for arbitrarily large values T a and T b . More importantly, we can achieve this with uncorrelated initial states, albeit paying for that with local coherence. Indeed, as discussed around Eq. (37), for arbitrary initial state ρ 0 , On the other hand, when mixing the basis steady states, the entanglement vanishes, especially when the states are only weakly entangled. Therefore, in the limit of high T a and T b , in order for ρ st to be entangled, one of the λ k 's has to be close to 1. Let us pick for definiteness λ 1 (the analysis for λ 2 and λ 3 will be identical). Then requiring λ 1 ≈ 1 will be equivalent to |θ σ ||θ κ | cos(arg θ σ + arg θ κ + 4π/3) ≈ 1.
Finally, let us show that, analogously to the current and heat flux, the negativity is also a 2π/3-periodic function of φ (see Fig. 5). Indeed, as we have established in Sec. V C, shifting the phase by 2π/3 transforms the steady state On the other hand, negativity is invariant under local unitary transformations on the state [59], which means that a 2π/3 phase shift does not alter the negativity.

B. Work content
The steady states of our machine can be useful not only for quantum information, but also for thermodynamics. In particular, the non-equilibrium steady states of the rotor store work that can be extracted from it in a cyclic Hamiltonian process. In other words, the average energy of the rotor, with respect to its Hamiltonian H, can be decreased by unitarily evolving its state ρ st . The unitary Global GKSL ME. The ergotropy of the steady state when the initial state is a thermal state at the temperature of the cold bath, versus the tunnelling coefficient τ . The plot also shows the total current in order to show that the amount of work stored in the state is not related to the current flowing in it. As an aside, the plot for the current once again illustrates the phenomenon of current inversion. The other parameters are the same as those in Fig. 3.
operation that extracts maximal work is the one that diagonalizes ρ st in the basis of H in such a way that, if H's eigenbasis is chosen so that E 1 ≤ E 2 ≤ · · · , then the eigenvalues of ρ st , {r k }, have the reverse ordering: r 1 ≥ r 2 ≥ · · · [72]. The average work thus extracted, E, is called ergotropy [72] and is equal to The states with zero ergotropy are called passive [72,73], with the most prominent example being provided by the thermal states [73].
When the rotor's initial state is thermal at temperature T in , the ergotropy of the steady state will be a decreasing function of T in , turning to zero somewhere in between T a and T b . Now, given the infinite size of the baths, thermal states at the temperatures of the baths can be considered to be for free -one just puts the system in contact with the bath and waits for it to thermalize. Therefore, if T in = T a , the machine essentially takes in a free (and useless from the perspective of work extraction) state and evolves (and maintains for an arbitrarily long time) it into a state "charged" by work. In Fig. 6, we plot E[ρ st [τ Ta ]] as a function of τ . We see that, in addition to a temperature gradient, a sufficiently high tunnelling rate is necessary for the machine to work as a charger. The plot also shows the total steady-state current, in order to indicate the surprising fact that the work content of the steady state is not related to (let alone being caused by) the current it harbours.

VII. CONCLUSIONS
We have studied a model thermal rotor, consisting of two strongly interacting partitions, in the quantum regime. Our main purpose was to study the transport properties of the rotor, specifically, the conditions under which there is non-zero particle current in the steady state of the rotor.
Before proceeding to analyzing the specifics of our model, we first had to define local particle current in the given situation, i.e., when the particle is a part of a larger quantum system immersed in a dissipative environment. We derive a novel general formula for the particle current applicable not only in our setup, but also to arbitrary open quantum networks. We furthermore provide an interpretation for our formula in terms of quantum weak measurements. Given the generality of the formula, this interpretation thus applies to such widely used definitions of current as the tunnelling current in Eq. (17).
Having defined the main quantity of interest, we go on analyzing the machine when its dynamics is described by both local and global quantum GKSL master equations. The former is applicable when the interaction between the particles is weak, while the latter is necessary when the particles are coupled strongly. In both cases, we study the steady state current from the perspective of the symmetries of the rotor Hamiltonian, and find that the particle current due to thermal hopping (J (th) ) can occur only when, in addition to a non-zero thermal gradient, the particle exchange symmetry of the Hamiltonian is broken. We also study the classical limit of the machine, and find that the same symmetry-breaking mechanism applies also there. Lastly, we show that the steady state current can reverse its direction as the tunnelling rate is varied. This is purely quantum effect, hitherto thought to take place only in the presence of strong external driving [3,57]. This is a somewhat counterintuitive effect, especially in the light of the fact that the tunnelling in our model is rotationally symmetric and thus, in itself, does not introduce a preferred direction for particle flow.
By studying the entanglement and ergotropy of the steady states, we furthermore prove that our machine can also be useful in such practical tasks as (i) converting uncorrelated states into entangled states in the presence of hot dissipative environment and (ii) evolving initially thermal states, which as useless for work extraction, into non-passive steady states, from which non-zero work can be extracted. Both tasks are enabled by purely quantum resources. In the first case, the entanglement generation is made possible by feeding the machine with coherent local states. In other words, the evolution to the steady state is a realization of a completely positive trace preserving map, converting local coherence into entanglement such that the entanglement on the output depends on the amount of coherence at the input, akin to the analysis in Ref. [74]. In the second case, "charging" of the initially thermal rotor is possible only if the tunnelling rate is sufficiently high.