Effect of Inter-System Coupling on Heat Transport in a Microscopic Collision Model

Abstract

In this paper we consider a bipartite system composed of two subsystems each coupled to its own thermal environment. Based on a collision model, we mainly study whether the approximation (i.e., the inter-system coupling is ignored when modeling the system–environment interaction) is valid or not. We also address the problem of heat transport unitedly for both excitation-conserving system–environment interactions and non-excitation-conserving system–environment interactions. For the former interaction, as the inter-system interaction strength increases, at first this approximation gets worse as expected, but then counter-intuitively gets better even for a stronger inter-system coupling. For the latter interaction with asymmetry, this approximation gets progressively worse. In this case we realize a perfect thermal rectification, and we cannot find an apparent rectification effect for the former interaction. Finally and more importantly, our results show that whether this approximation is valid or not is closely related to the quantum correlations between the subsystems, i.e., the weaker the quantum correlations, the more justified the approximation and vice versa.

1. Introduction

In most practical situations, a quantum system inevitably interacts with its environment, which induces decoherence and dissipation [1]. In this case, its dynamics are usually described by the Gorini–Kossakowski–Lindblad–Sudarshan (GKLS) quantum master equation with a series of approximations [2]. When deriving a master equation for an open quantum system, one may always obtain a global master equation, which considers the full system Hamiltonian, i.e., including the direct coupling between the subsystems [3,4,5,6]. However, such derivation is very complicated when the system is composed of two or more interacting subsystems. Hence a local master equation that ignores the direct interactions between the subsystems is often used as a substitute [3,4,5,6,7,8,9,10,11]. By means of these two kinds of master equations, a lot of efforts has recently been devoted to the heat transport for thermodynamic systems [7,8,9,10,11,12,13,14,15]. However, by comparing the dynamics resulting from the corresponding master equations with exact numerical simulations, both approaches may lead to seeming thermodynamic inconsistencies or just suit to some parameter regimes, as pointed out in References [15,16,17,18,19,20,21,22]. In Reference [3], the local description for two coupled quantum nodes may predict heat currents from a cold to a hot thermal reservoir, or the existence of currents even in the absence of a temperature gradient. The origin of these effects, as discussed in References [5,23], lies in the fact that there is an external work cost related to the breaking of global detailed balance. By including this work cost, this inconsistencies can be resolved. For weak inter-system coupling, it was shown that the global approach fails in non-equilibrium situations, whereas the local approach agrees with the exact solution [19]. In Reference [21], due to a failure of the secular approximation, it was reported that the global master equation erroneously gave a vanishing heat current through a spin-$\frac{1}{2}$ Heisenberg chain in the presence of a finite-temperature gradient. For two coupled qubits interacting with common and separate baths, it was shown that the global approach with partial secular rather than full secular approximation always provides the most accurate choice for the master equation [24]. In Reference [25], it was shown that the completely positive version of the Redfield equation and an appropriate time-dependent convex mixture of the local and global solutions gives rise to the most accurate semigroup approximations of the exact system dynamics. Recently, it was proved that local master equations are consistent with thermodynamics without resorting to a microscopic model [26]. Moreover, the two approaches above, from the viewpoint of the system, often rely on approximate Markovian master equations derived under the assumptions of weak system–environment coupling, which would become challenging under strong coupling.

Furthermore, the manipulation of heat transport in non-equilibrium steady-state has been identified as one of the crucial studies of quantum thermodynamics, which gives us an improved understanding of classical thermodynamics in quantum domain [27,28,29,30,31,32,33,34,35]. For example, heat transport between two bosonic reservoirs was predicted for a coupled two-state system, and a formula for thermal conductance was derived based on a rate equation formalism [33]. Bandyopadhyay et al. proposed a numerical scheme for exactly simulating the heat transport in a quantum harmonic chain with self-consistent reservoirs [34]. By means of an effective harmonic Hamiltonian, a quantum thermal transport through anharmonic systems was studied within the framework of the nonequilibrium Green’s function method [35]. Besides, quantum devices such as heat rectifier, thermal memory, and thermal ratchet, have also become goals of controlling thermal transport in quantum thermodynamics [36,37,38,39,40,41,42,43,44,45]. It was found that, by using the quantum master equation, thermal rectification in anisotropic Heisenberg spin chains could change sign when the external homogeneous magnetic field was varied [43]. An optimal rectification in the ultrastrong-coupling regime of two coupled two-level systems was shown [44]. In Reference [45], Jose et al. studied two interacting spin-like systems characterized by different excitation frequencies, which can be used as a quantum thermal diode.

Recently, the collision model, also called repeated interactions, has drawn attention for its potential advantage in simulating open quantum system [46,47,48,49,50,51,52,53,54,55,56,57]. It was assumed that the environment consists of a large collection ancillae and the system of interest interacts, or collides, with an ancilla at each time step. In the framework of collision model, a continuous-time description in terms of a Lindblad master equation can be derived in the short-time limit provided some assumptions are made about the system–ancilla interaction [40,49]. For instance, in Reference [58] the system’s dynamics embodied by the stroboscopic map can be approximated by a Lindblad master equation in a short-time limit. In a similar way, in Reference [5] a local master equations with a Lindblad form for two coupled harmonic oscillator was derived by using the method of repeated interactions. Moreover, also from the viewpoint of the system, the corresponding reduced dynamics can be obtained in many cases without any approximations [59,60,61,62,63,64,65]. This is because collision model allows for the possibility to decompose a complicated open dynamics in terms of discrete elementary processes. It is particularly suited for addressing the thermodynamics of engineered reservoirs [23,48,49,62,66,67]. For example, under an energy-preserving system–environment interaction within the framework of collision model, it has been found that the non-monotonic time behavior of the heat exchange between the system and environment could serve as an indicator of non-Markovian behaviors [62]. Besides, based on a microscopic collision model, heat current between a coupled system can flow from the cold nonthermal reservoir to the hot one due to the contribution of coherence [66]. A link between information and thermodynamics, for a multipartite open quantum system with a finite temperature reservoir, was displayed in term of a collision model [67].

In this paper, we consider a two coupled qubit system that interacts with its local environment consisting of a large collection of identical ancillae. Inspired by the previous work on the local and global master equations, we consider an decoupling approximation within the framework of collision model, specifically, environment acts on subsystem without considering the inter-system coupling, i.e., ignores the direct coupling between the subsystems when modeling the system–environment interaction. It is noted that although this approximation in this collision model is not strictly a substitute for the local master equation we described above, we expect that the results derived from this simple and solvable model can provide a reference. We mainly study whether this decoupling approximation is valid or not, i.e., whether or when this direct coupling can be ignored. We systematically examine the heat transport through the system in both weak and strong system–ancilla coupling. Within the framework of collision model, we give a more general definition of heat current under non-energy preserving system–ancilla interactions (there is external work performing on the whole system–ancilla compound). In the case of excitation-conserving system–ancilla interactions, we find that the results predicted with the decoupling approximation is valid even at strong inter-system interactions. In particular, we realize a perfect rectification in the asymmetric systems and give a discussion about the mechanism of this phenomenon. Moreover, we find that whether or not this decoupling approximation is valid is closely related to the quantum correlations between the subsystems.

2. Model and Methods

We consider a bipartite system S consisting of two identical two-level subsystems $S_i$$(i=1,2)$, and each subsystem is coupled to its local thermal subenvironment $E_i$ set at temperature $T_i$. Here the subenvironment $E_i$ is a sequence of non-interacting ancillae $(E_1^i,E_2^i,...,E_n^i)$ all in the same initial state ${\eta }_n^i$. Subsystem $S_i$ interacts with the connected subenvironment $E_i$ via a series of short subsystem–ancilla interactions. The joint state of system and environment is initially factorized:

$${\rho }^{SE}\left(0\right)={\rho }^S\left(0\right)\left({\otimes }_{j=1}^n{\eta }_j^1\right)\left({\otimes }_{j=1}^n{\eta }_j^2\right)$$

(1)

where ${\rho }^S\left(0\right)$ and ${\otimes }_{j=1}^n{\eta }_j^i$ are the initial states of the open system and subenvironment $E_i$, respectively.

In Figure 1 we show a schematic sketch of the collision model considered. It is composed of a series of repeated rounds. In any one round of the collision process, first the whole system S undergoes a free evolution lasting a time interval $\tau$, subsequently $S_1$ and $S_2$ locally interact with only one ancilla of its environment, respectively. Then by tracing out the environment’s degrees of freedom and repeating the above process in the next round, we can obtain the reduced state of the system in the full time evolution. Therefore, each round is composed of three steps: a free evolution of S and two subsystem–ancilla interactions. It is noted that the subenvironment is assumed to be large enough, so that the subsystem never collides twice with the same ancilla. As a consequence, at each collision round n, the subsystem $S_i$ collides with a ”fresh” $E_n^i$.

We assume throughout this paper that S and each ancilla $E_n^i$ of the subenvironment $E_i$ are qubits with logical states $\left\{\right|0\rangle ,|1\rangle \}$. The corresponding free Hamiltonians for the subsystem $S_i$$(i=1,2)$ and the ancilla $E_n^i$ are ${\widehat{H}}_{S_i}=\frac{1}{2}{\omega }_i{\widehat{\sigma }}_z$ and ${\widehat{H}}_{E_n^i}=\frac{1}{2}{\omega }_0{\widehat{\sigma }}_z$, respectively. Here ${\widehat{\sigma }}_z$ is the usual Pauli operator (we set $\hslash =1$). The free dynamics of the whole system S are described by the unitary evolution operator

$${\widehat{U}}_{S_1,S_2}=exp[-i({\widehat{H}}_0+{\widehat{H}}_{int}^{S_1,S_2})\tau ],$$

(2)

with ${\widehat{H}}_0={\widehat{H}}_{S_1}+{\widehat{H}}_{S_2}$ and the inter-system interaction Hamiltonian ${\widehat{H}}_{int}^{S_1,S_2}$. We use unitary operator ${\widehat{V}}_{S_i,E_n^i}$ to model the collisions between the subsystem $S_i$ and ancilla $E_n^i$ of subenvironment $E_i$ (Its exact definition will be given later in different cases); it is assumed that all the subsystem–ancilla collisions have the same duration $\tau$. So the dynamical map ${\widehat{\mathsf{\Lambda }}}_{S_1,S_2}$ that governs the free evolution of the system S and map ${\widehat{\Psi }}_{S_1,E_n^i}$ that governs the system–subenvironment interaction at the nth round can be written as

$${\widehat{\mathsf{\Lambda }}}_{S_1,S_2}\left(\rho \right)={\widehat{U}}_{S_1,S_2}\rho {\widehat{U}}_{S_1,S_2}^{†},$$

(3)

$${\widehat{\Psi }}_{S_i,E_n^i}\left(\rho \right)={\widehat{V}}_{S_i,E_n^i}\rho {\widehat{V}}_{S_i,E_n^i}^{†},$$

(4)

respectively. Following the repeated interaction approach mentioned above, the joint state of system and ancillae is brought from the nth round to the $(n+1)$th round through the process

$${\rho }_{n-1}^S\otimes {\eta }_n^1\otimes {\eta }_n^2\to {\rho }_n^{SE}=\widehat{\mathcal{U}}[{\rho }_{n-1}^S\otimes {\eta }_n^1\otimes {\eta }_n^2]{\widehat{\mathcal{U}}}^{†},$$

(5)

where $\widehat{\mathcal{U}}={\widehat{V}}_{S_2,E_n^2}{\widehat{V}}_{S_1,E_n^1}{\widehat{U}}_{S_1,S_2}$. Hence after round n, the reduced density matrix of the system ${\rho }_n^S$ is

$${\rho }_n^S={\mathrm{Tr}}_{E_n^1,E_n^2}\left[{\rho }_n^{SE}\right],$$

(6)

where ${\mathrm{Tr}}_{E_n^1,E_n^2}[·]$ denotes the partial trace over the two ancillae $E_n^1$ and $E_n^2$. Similarly, the reduced state ${\tilde{\eta }}_n^{1\left(2\right)}$ of the nth ancilla of subenvironment $E_{1\left(2\right)}$ is

$${\tilde{\eta }}_n^{1\left(2\right)}={\mathrm{Tr}}_{S{E}_n^{2\left(1\right)}}\left[{\rho }_n^{SE}\right].$$

(7)

Throughout we assume each ancilla of subenvironment $E_i$ to be initially in a thermal state with inverse temperature ${\beta }_i=1/T_i$ (we set $k=1$), namely,

$${\eta }_n^i=\frac{1}{Z}exp(-{\beta }_i{\widehat{H}}_{E_n^i}),$$

(8)

where $Z=\mathrm{Tr}[exp(-{\beta }_i{\widehat{H}}_{E_n^i})]$ is the partition function.

In the nth round of the dynamics, the whole system S undergoes a free evolution; $S_1$ interacts with $E_n^1$; next, $S_2$ interacts with $E_n^2$; then $S_1$ and $S_2$ shift by one site to the $(n+1)$th round in which also after the free evolution of system $S_1$-$E_{n+1}^1$ and $S_2$-$E_{n+1}^2$ interactions subsequently take place. This process will repeat to the $(n+2)$th round. The free evolution is governed by Equations (2), and now we begin to give the subsystem–ancilla evolution operator. When ancilla $E_n^i$ collides with subsystem $S_i$ not ignoring the direct interaction between the two subsystems, the corresponding unitary time evolution operator can be written as

$${\widehat{V}}_{S_i,E_n^i}=exp[-i\left({\widehat{H}}_{S_1}+{\widehat{H}}_{int}^{S_1,S_2}+{\widehat{H}}_{S_2}+{\widehat{H}}_{E_n^i}+{\widehat{H}}_{int}^{S_i,E_n^i}\right)\tau ],$$

(9)

where ${\widehat{H}}_{int}^{S_i,E_n^i}$ is the interaction Hamiltonian between the subsystem $S_i$ and ancilla $E_n^i$. However, it is convenient to ignore the direct interaction between the two subsystems, namely, a decoupling approximation neglects the coupling between the subsystems when modeling the system–environment interaction. In this case, the corresponding unitary time evolution operator is written as

$${\widehat{V}}_{S_i,E_n^i}^{app}=exp[-i\left({\widehat{H}}_{S_i}+{\widehat{H}}_{E_n^i}+{\widehat{H}}_{int}^{S_i,E_n^i}\right)\tau ].$$

(10)

It is obvious that the difference between Equations (9) and (10) resides in the inter-system interaction in each round: Equation (9) arises naturally when modeling the subenvironment–subsystem coupling from a microscopic model considering the full system Hamiltonian, i.e., not ignore the direct interaction between its subsystems, while Equation (10) ignores this interaction, i.e., the decoupling approximation.

3. Symmetric System

In this section we consider a symmetric system, from which to investigate whether the inter-system interaction can be ignored or not. Among the possible choices for the interaction between the subsystem $S_i$ and their connected ancilla $E_n^i$, we choose the excitation-conserving interaction Hamiltonian as

$${\widehat{H}}_{int}^{S_i,E_n^i}=\gamma ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y)$$

(11)

with subsystem–subenvironment interaction strength $\gamma$. The interaction between $S_1$ and $S_2$ takes the same form

$${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y)$$

(12)

with coupling strength $\delta$. We consider that the energy gaps of subsystems are the same, i.e., ${\omega }_1={\omega }_2={\omega }_0$.

3.1. Heat Current

Now we investigate the heat currents to evaluate the performance of the decoupling approximation. We also consider the heat currents in the case of not ignoring the inter-system interactions (i.e., without the decoupling approximation), which serves as a benchmark. In the case of decoupling approximation, for the system characterized by Hamiltonians of Equations (11) and (12), its corresponding unitary system–environment operator preserves the energy, i.e.,

$$[{\widehat{V}}_{S_i,E_n^i}^{ign},{\widehat{H}}_{S_i}+{\widehat{H}}_{E_i}]=0.$$

(13)

This implies that all the energy leaving the ancilla enters the system. Hence, during the $(n+1)$th round, the heat exchange $▵Q$ between system $S_{1\left(2\right)}$ and ancilla $E_n^{1\left(2\right)}$ is given by

$$▵Q_{E_n^{1\left(2\right)}}=\mathrm{Tr}\left[{\widehat{H}}_{E_n^{1\left(2\right)}}({\tilde{\eta }}_n^{1\left(2\right)}-{\eta }_n^{1\left(2\right)})\right].$$

(14)

Due to $▵Q_{E_n^1}=-▵Q_{E_n^2}$ for steady state, the stationary heat current flowing from subenvironment $E_1$ to the system $S_1$ can be defined as

$$J_h=-▵Q_{E_n^1}.$$

(15)

This definition was often used when Equation (13) is satisfied [42,43,44,45,58,66]. However, in the case of not ignoring the inter-system interactions, the corresponding unitary system–environment operator no longer preserves the energy because of

$$[{\widehat{V}}_{S_i,E_n^i}^{con},{\widehat{H}}_{S_1}+{\widehat{H}}_{S_2}+{\widehat{H}}_{int}^{S_1,S_2}+{\widehat{H}}_{E_i}]\ne 0.$$

(16)

That is, not all the energy leaving the ancilla enters the system. So an external work is required to turn the system–ancilla coupling on and off [67,68], which we refer to as the switching work (labeled $W_{SW}$). In each round of collisions illustrated in Figure 1, the corresponding work $W_{SW}$ on the system can be written as

$$W_{SW}=W_1+W_2.$$

(17)

The first term $W_1$ is from the sudden on and off switching of $S_1-E_n^1$ interaction, and it reads

$$W_1=\mathrm{Tr}\left[{\widehat{H}}_{int}^{S_1,E_n^1}({\rho }_{S_1{S}_2}^{\prime }\otimes {\eta }_n^1-{\rho }_{S_1{S}_2{E}_n^1})\right],$$

(18)

where ${\rho }_{S_1{S}_2}^{\prime }$ is the state of system before $S_1-E_n^1$ collision, and ${\rho }_{S_1{S}_2{E}_n^1}$ is the global state of system S and ancilla $E_n^1$ after $S_2-E_n^1$ collision. Similarly, the second term associated with $S_2-E_n^2$ interaction takes the form

$$W_2=\mathrm{Tr}\left[{\widehat{H}}_{int}^{S_2,E_n^2}({\rho }_{S_1{S}_2}^{\prime \prime }\otimes {\eta }_n^2-{\rho }_{S_1{S}_2{E}_n^2})\right],$$

(19)

where ${\rho }_{S_1{S}_2}^{\prime \prime }$ is the reduced state of system after $S_1-E_n^1$ collision, and ${\rho }_{S_1{S}_2{E}_n^2}$ is the global state of system S and ancilla $E_n^2$ after $S_2-E_n^2$ collision. As a result, Equation (15) is no longer valid to calculate the heat current in the case of not ignoring the inter-system interactions.

Based on this collision model, we focus on deriving the expression of heat current that can be applicable in the cases of ignoring and not ignoring the inter-system interactions, as well as weak and strong system–environment couplings. Physically, heat current is defined as the energy going through the system. We consider the energy change of the system in a complete round (i.e., free evolution of S, $S_1-E_n^1$ and $S_2-E_n^2$ interactions). The energy change $▵E_S^1$ of system S in each $S_1-E_n^1$ interaction can be written as

$$▵E_S^1=\mathrm{Tr}\left[({\widehat{H}}_0+{\widehat{H}}_{int}^{S_1{S}_2})({\rho }_{S_1{S}_2}^{\prime \prime }-{\rho }_{S_1{S}_2}^{\prime })\right].$$

(20)

Similar expressions can be obtained for $S_2-E_n^2$ interaction, and the corresponding energy change $▵E_S^2$ of system S is given by

$$▵E_S^2=\mathrm{Tr}\left[({\widehat{H}}_0+{\widehat{H}}_{int}^{S_1{S}_2})({\rho }_{S_1{S}_2}^{\prime \prime \prime }-{\rho }_{S_1{S}_2}^{\prime \prime })\right],$$

(21)

where ${\rho }_{S_1{S}_2}^{\prime \prime \prime }$ is the reduced state of system S after $S_2-E_n^2$ collision. Since the inter-system dynamics is unitary, the energy of system is preserved in each free evolution. For steady state, the state of the system cannot change in a complete round, leading to

$$▵E_S^1=-▵E_S^2.$$

(22)

From the above discussion, it is clear that $▵E_S^1$ is the energy going through the system. Thus, the stationary heat current $J_h$ flowing from $E_1$ to $E_2$ can be rewritten as

$$J_h=▵E_S^1.$$

(23)

In the case of ignoring the inter-system interactions, if the corresponding system satisfies the commutation condition Equation (13), i.e., all energy changes in the system can be attributed to energy flowing to or from the ancilla, it is thus clear that Equations (23) reduces to Equation (15).

Here we should mention that the exchange of energy between two quantum systems, characterized by particular commutation relation between the local Hamiltonians and the interaction operator, can always be split into work and heat [69,70]. This formalism was extended to situations where one, or both, subsystems are coupled to a thermal environment. In the case of not ignoring the inter-system interaction, according to Reference [70] the stationary heat current and work in each $S_1-E_n^1$ collision are calculated as

$$J_h^{\prime }=-i{\int }_0^{\tau }\mathrm{Tr}[\mathbb{I}\otimes ({\widehat{H}}_0+{\widehat{H}}_{int}^{S_1,S_2}),{\widehat{H}}_{int}^{S_1,E_n^1}\otimes \mathbb{I}]C_{12}dt,$$

(24)

$$W=-i{\int }_0^{\tau }\mathrm{Tr}[{\widehat{H}}_S^{eff},{\widehat{H}}_S]{\tilde{\eta }}_n^{1}dt,$$

(25)

where $C_{12}$ represents the correlations between the system S and $E_n^1$ after $S_1-E_n^1$ collision, and ${\widehat{H}}_S^{eff}$ is the diagonal part of matrix ${\mathrm{Tr}}_{E_n^1}\left[{\widehat{H}}_{int}^{S_1,E_n^1}({\tilde{\eta }}_n^1\otimes \mathbb{I})\right]$. From numerical calculations, we check that the heat current $J_h^{\prime }$ given by Equation (24) is equal to that given by Equation (23). Thus our previous definition of Equation (23) is consistent with Reference [70].

In Figure 2 we plot $J_h$ as a function of the inter-system coupling strength $\delta$ for various subsystem–ancilla coupling strength $\gamma$ in the cases with and without the decoupling approximation. We fix the initial state of the system to be $|11\rangle$ and set $T_1=5{\omega }_0$ and $T_2={\omega }_0$ for subenvironments $E_1$ and $E_2$, respectively. It can be seen from Figure 2 that, for fixed $\gamma$, the heat currents in both cases increase rapidly with the increase of $\delta$ and they eventually reach the same steady value, before which the heat current from the decoupling approximation is always smaller than that without the approximation. In general, for weak inter-system coupling, one would expect smaller deviation from the case without the approximation. As expected, in the limit $\delta \to 0$, it gives the same (zero) heat current as that without the approximation. For fixed $\gamma$, when $\delta$ increases we observe that the heat current predicted with the decoupling approximation gradually deviates from that without the approximation, and this deviation quickly increases from zero to its maximum. However, it might seem counter-intuitive that this deviation gradually decreases with $\delta$ and eventually vanishes at larger values of $\delta$. That is, for larger $\delta$, the heat currents from the decoupling approximation can still be consistent with that without the approximation, and this decoupling approximation consequently may not necessarily break down.

Moreover, we find that the difference of the heat currents between those with and without the decoupling approximation strongly depends on $\gamma$, as shown in Figure 2. A bigger difference is obtained for stronger coupling strength $\gamma$. For instance, for $\gamma =0.8{\omega }_0$, the difference that maintains nonzero values is within a larger region of $\delta$, compared with $\gamma =0.2{\omega }_0$ and $0.5{\omega }_0$. That is, when assessing the heat currents the results from the decoupling approximation is inconsistent with that without this approximation. Smaller difference is obtained for weaker coupling strength $\gamma$. For $\gamma =0.2{\omega }_0$, it can be seen from the inset of Figure 2 that the difference between two cases maintaining nonzero values is within a very small region below $\delta \sim 0.04{\omega }_0$. In other words, for smaller values of $\gamma$, the results from the decoupling approximation agrees well with that without this approximation, and consequently it also approximately predicts the correct heat currents within almost all region of $\delta$. Physically, this can be easily understood as following: decreasing the system–ancilla coupling strength would weaken the influence of the inter-system coupling ignored.

3.2. Trace Distance

When investigating the accuracy of the local and global master equations, steady state is often used as a reference for predicting the results of the nonequilibrium dynamics [16,19]. To further assess the validity of decoupling approximation, we also consider the obtained steady states. As a measure of distinguish ability, here we examine the trace distance between the steady state obtained from the case with the decoupling approximation and the steady state obtained without this decoupling approximation:

$${\mathcal{D}}_{\mathcal{T}}=\frac{1}{2}{\Vert {\rho }^S-{\rho }^{Sapp}\Vert }_1,$$

(26)

where ${\Vert ·\Vert }_1$ is the trace norm while ${\rho }^S$ and ${\rho }^{Sapp}$ are the reduced steady states of the system with and without the decoupling approximation, respectively. The trace distance is equal to unity for fully distinguishable states, while it is null for identical states. Figure 3 shows the dependence of ${\mathcal{D}}_{\mathcal{T}}$ on $\delta$ for different strength $\gamma$. For fixed $\gamma$, it can be seen that ${\mathcal{D}}_{\mathcal{T}}$ quickly increases and then gradually decreases with the increase of $\delta$ (e.g., for $\gamma =0.1{\omega }_0$). That is, when $\delta$ increases, the steady state ${\rho }_n^{Sapp}$ predicted with the decoupling approximation gradually deviates from the steady state ${\rho }_n^S$ predicted without this approximation, then this deviation reaches its maximum and eventually decreases. It is obvious that this behavior of the steady state from this approximation is indeed similar to that of heat current. Here we also arrive at a conclusion similar to that of heat current: even for a large $\delta$, decoupling approximation can be well justified when predicting the steady state. The inset of Figure 3 plots the value of ${\delta }^{\prime }$ for the maximum of ${\mathcal{D}}_{\mathcal{T}}$ as a function of $\gamma$ (i.e., when $\delta ={\delta }^{\prime }$, ${\mathcal{D}}_{\mathcal{T}}$ reaches its maximum for fixed $\gamma$). It can be seen that ${\delta }^{\prime }$ nonlinearly increases with the subenvironment–ancilla interaction strength $\gamma$.

Figure 3 also shows the effect of system–ancilla coupling strength $\gamma$ on ${\mathcal{D}}_{\mathcal{T}}$. For $\gamma =0.2{\omega }_0$, it can be seen that ${\mathcal{D}}_{\mathcal{T}}$ maintains a larger value only within a smaller region of $\delta$, and quickly decreases compared to $\gamma =0.5{\omega }_0$ and $\gamma =0.8{\omega }_0$. Again, for weak $\gamma$, steady state predicted with the decoupling approximation agrees well with that without the approximation even at strong interaction strengths $\delta$. In other words, this approximation is well justified to describe steady state in this case.

3.3. Quantum Discord

Why can the decoupling approximation give a good estimate of heat current and steady state even for a large $\delta$ ? In References [71,72], it was found that the compositeness (two particles behave like a single particle) is closely related to the quantum correlations between the constituent particles. It is thus very interesting to consider the quantum correlations between the bipartite system $S_1$ and $S_2$, which can be captured by the discord [73]

$$\mathcal{D}\left({\rho }_S\right)={min}_{{\mathsf{\Pi }}_A}\{I\left({\rho }_S\right)-I\left({\mathsf{\Pi }}_A{\rho }_S\right)\}.$$

(27)

Here ${\mathsf{\Pi }}_A$ is a set of rank-one POVM projectors on system $S_1$, and $I\left({\rho }_S\right)=S\left({\rho }_S\right)-S\left({\rho }_{S_1}\right)-S\left({\rho }_{S_2}\right)$ is the quantum mutual information associated with the von Neumann entropy. Without the decoupling approximation, in Figure 4 we plot $\mathcal{D}$ as a function of $\delta$. As inter-system coupling strength $\delta$ increases it can be seen that, for fixed $\gamma$, $\mathcal{D}$ first increases from zero to its maximum, then it gradually decreases, i.e., the quantum correlations between two subsystems first increases and then decreases. It is obvious that such behavior is associated with the heat current or the steady state, i.e., the greater the quantum correlation $\mathcal{D}$ between two subsystems, the greater the deviation of the heat current and steady state predicted with the decoupling approximation, and vice versa. In general, one would think that whether this approximation is valid or not should depend on inter-system coupling strength. In fact, inter-system interactions are really necessary and they serve as physical means to create quantum correlations, while it is noted that such correlations can still be present in spatially separated subsystem that no longer interact after it is created. When the inter-system coupling produces a stronger correlations, these two subsystems behave more like a composite and can not be treated separately, so this decoupling approximation gets worse. The fact that decoupling approximation can give an good description even at stronger inter-system coupling is because the stronger inter-system interaction does not produce a strong enough quantum correlations, as shown in Figure 4. To summarize, this approximation depends strongly on quantum correlations between the two qubits, i.e., the higher the quantum correlations, the worse the approximation, and vice versa.

Moreover, for this model, we also consider the other situations by replacing Equations (11) and (12) with any of the following three combinations of excitation-conserving interactions:

• ${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y+{\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z)$ and ${\widehat{H}}_{int}^{S_i,E_n^i}=\gamma ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y)$.
• ${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y)$ and ${\widehat{H}}_{int}^{S_i,E_n^i}=\gamma ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y+{\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z)$.
• ${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y+{\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z)$ and ${\widehat{H}}_{int}^{S_i,E_n^i}=\gamma ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y+{\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z)$.

For these three combinations, we also investigate the heat current, trace distance, and quantum discord with and without ignoring the inter-system interactions. The corresponding results are similar to those from Equations (11) and (12). It is noted that, with the decoupling approximation, the corresponding unitary operators of all three combinations above, like Equations (11) and (12), also satisfy the condition $[{\widehat{V}}_{S_i,E_n^i}^{ign},{\widehat{H}}_{S_i}+{\widehat{H}}_{E_i}]=0$, i.e., they correspond to energy-preserving system–environment interactions and no external work exists. Without the approximation, in each round (i.e., free evolution of S, $S_1-E_n^1$ and $S_2-E_n^2$ interactions) although there is an external work due to $[{\widehat{V}}_{S_i,E_n^i}^{con},{\widehat{H}}_{S_1}+{\widehat{H}}_{S_2}+{\widehat{H}}_{int}^{S_1,S_2}+{\widehat{H}}_{E_i}]\ne 0$, after some calculations we find that this external work in each round is much smaller than the heat flowing through the system in each round.

4. Asymmetric System

In this section, we turn our attention to the systems featured by distinct asymmetries, such as the qubits with different frequencies or asymmetric couplings to their reservoirs. Indeed, inspired by the manipulation and control of the thermal transport in micro-scale, such asymmetric systems were widely investigated in the thermal diode and thermal transistor [37,38,39,40,41,42,43,44,45]. In the following, we will study the heat transport and, at the same time, study whether the inter-system interaction can be ignored or not when modeling the system–environment interaction.

4.1. Off-Resonant Interacting Qubits

First we consider two qubits are characterized by different energy gaps ${\omega }_1$ and ${\omega }_2$, i.e., ${\omega }_1\ne {\omega }_2$. The interaction Hamiltonian between two subsystems $S_1$ and $S_2$ is given as

$${\widehat{H}}_{int}^{S_1,S_2}=\delta {\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z.$$

(28)

Such system may arise, for example, when a nonuniform magnetic field is applied to a pair of interacting spins in the z directions. The coupling to each ancilla of environment is described by the following non-excitation-conserving Hamiltonian:

$${\widehat{H}}_{int}^{S_i,E_n^i}=\gamma {\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x.$$

(29)

For this interaction, in both cases with and without the decoupling approximation, not all the energy leaving the environment enters the system due to $[{\widehat{V}}_{S_i,E_n^i}^{ign},{\widehat{H}}_{S_i}+{\widehat{H}}_{E_i}]\ne 0$ and $[{\widehat{V}}_{S_i,E_n^i}^{con},{\widehat{H}}_{S_1}+{\widehat{H}}_{S_2}+{\widehat{H}}_{int}^{S_1,S_2}+{\widehat{H}}_{E_i}]\ne 0$.

From Equation (23) we calculate the stationary heat currents for system characterized by Hamiltonians of Equations (28) and (29), also addressing with and without the decoupling approximation separately. From our study we find that the heat current from the decoupling approximation is always zero even in the presence of a finite temperature gradient $T_1>T_2$ or $T_2>T_1$, regardless of the coupling strength $\delta$ and $\gamma$. This result deviates from that without the approximation, in which heat currents increases as $\delta$ increases, so the decoupling approximation gets progressively worse with the increase of $\delta$.

Moreover, for this model we compare the heat currents for various temperature difference $▵T$ between two environments, which shows an asymmetric heat transport of the system (i.e., heat rectification). To quantify this rectification efficiency, we also give a rectification factor R defined as follows [43,44,45]:

$$R=\frac{|J_h\left(▵T\right)+J_h(-▵T)|}{max\left(\right|J_h\left(▵T\right)|,|J_h(-▵T)\left|\right)},$$

(30)

where $J_h\left(▵T\right)$ is the forward heat current for $T_1>T_2$, and $J_h(-▵T)$ is the backward heat current when the temperature gradient is reversed. It can be seen from Figure 5 that in the case $T_1=10{\omega }_0$ and $T_2=8{\omega }_0$, its corresponding forward heat current is greater than the backward heat currents ($T_2=10{\omega }_0$ and $T_1=8{\omega }_0$). Hence an asymmetric conduction ($|J_h\left(▵T\right)|\ne |J_h(-▵T)|$) emerges, which varies with the increase of $\delta$. If the temperature gradient raises, such as $T_1=10{\omega }_0$ and $T_2=4{\omega }_0$, this asymmetric conduction increases, i.e., the higher the gradient, the stronger the heat rectification. Especially, when $T_1=10{\omega }_0$ and $T_2=0.1{\omega }_0$, its corresponding backward heat current (cf. green dashed line corresponds to $T_1=0.1{\omega }_0$ and $T_2=10{\omega }_0$ in Figure 5) is almost zero for any $\delta$, i.e., an optimal rectification is realized with $R>0.996$. Insets in Figure 5 plot factor R as a function of $\delta$ for various temperature differences $▵T$. It can be seen that R increases with the increase of $▵T$.

The above rectification phenomena can be explained as follows. There are two ingredients affecting the rectification effect. One is nonzero heat currents in the absence of temperature gradient and the other is the temperature difference $▵T$ between two environments. (i) It can be seen from Figure 6 that the case without the decoupling approximation can give $J_h>0$ even when the two environments are at the same temperature (cf. black dashed line corresponds to $T_1=T_2=10{\omega }_0$ in Figure 6). It seems to violate the second law of thermodynamics but is justified. This is because there is an external work which is closely related to the internal subsystem–ancilla interaction. By taking this into account, the thermodynamic consistency can be guaranteed. We also find that the heat transport through the two-qubit system is dramatically affected by this work cost and can even be inverted, leading to $J_h>0$ for $T_1<T_2$, i.e., heat current from a cold to a hot bath (cf. dashed lines in Figure 5). Similarly, in Reference [5] heat transport through a chain of harmonic oscillator is investigated by using a local master equation based on repeated interactions, and all thermodynamic inconsistencies (such as a heat current from the cold to the hot bath) can be resolved correctly with the consideration of external work turning the subsystem–ancilla coupling on and off. (ii) Because there is a non-zero heat current in the absence of temperature gradient (labeled $J_0$), so when $T_1>T_2$, it can be seen from Figure 6 that the corresponding heat current ($J_h>J_0)$ increases with the increase of temperature gradient (cf. blue line in Figure 6). However, when $T_1<T_2$, as the temperature gradient increases, the corresponding heat current decreases from positive to zero (cf. red line in Figure 6). If the temperature difference increases further (if it is large enough), the corresponding heat current can change from positive to negative (cf. green line in Figure 6).

Figure 7 plots the trace distance ${\mathcal{D}}_{\mathcal{T}}$ between the density matrices with and without the decoupling approximation against $\delta$ for steady state. As expected, for $\delta =0$, we observe that ${\mathcal{D}}_{\mathcal{T}}=0$, which means that this approximation gives the same steady state as that without the approximation. For $\delta \ne 0$, ${\mathcal{D}}_{\mathcal{T}}$ increases with the increase of $\delta$. Therefore, the steady state from the decoupling approximation would deviates from that without the approximation, and this deviation increases as $\delta$ increases.

Now we consider quantum discord of the bipartite system. Without the decoupling approximation, in Figure 8 we plot $\mathcal{D}$, i.e., quantum correlations of two qubits, as a function of $\delta$ at steady state. It is clear that $\mathcal{D}$ increases with the increase of $\delta$. This behavior can also be closely related to that of Figure 5 (or Figure 7), i.e., as $\delta$ increases, the heat currents (or the steady state) predicted with the decoupling approximation gradually deviates from that without the approximation. Unlike Section 3.3 where quantum correlations first increase and then decrease with the increase of inter-system coupling, it can be seen from Figure 8 that increasing inter-system coupling enforces stronger correlations, which means these two qubits behave more like a composite and must be treated as a whole, so the decoupling approximation gets worse with the increase of the inter-system interaction. This further confirms our conclusion in Section 3.3, i.e., the higher the quantum correlations between two qubits, the greater the deviation from not ignoring the inter-system interaction.

In addition, we also consider the other situations by replacing Equation (28) with ${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y)$ or ${\widehat{H}}_{int}^{S_1,S_2}=\delta ({\widehat{\sigma }}_x\otimes {\widehat{\sigma }}_x+{\widehat{\sigma }}_y\otimes {\widehat{\sigma }}_y+{\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_z)$. In spite of a quantitative difference, the corresponding results are similar to those obtained above.

In Section 3 we only consider a resonant two-qubit system (it is noted that the system–environment interactions we chose in Section 3 are different from Equation (29) in this section). For a more comprehensive study of heat transfer, now we consider two off-resonant qubits under the same interactions as in Section 3. From our numerical calculations we find that the corresponding results are similar to those obtained in Section 3, and are different from the results in this section above. That is, as $\delta$ increases, heat current predicted with the decoupling approximation gradually deviates from that without the approximation, then this deviation decreases at larger values of $\delta$. Moreover, for the system–environment interactions in Section 3, we find that although there is also an asymmetric heat conduction in the case of without the approximation, no obvious thermal rectification emerges. This shows that thermal rectification are not only related to asymmetry of the system, but also related to the system–environment interactions.

4.2. Anisotropic Interacting Qubits

Connecting two qubits $S_1$ and $S_2$ by anisotropic interactions is another possible way to introduce asymmetry that may also exhibit rectification. Here we consider the case where the inter-system coupling is taken as an anisotropic exchange interaction [42]

$${\widehat{H}}_{int}^{S_1,S_2}=\delta {\widehat{\sigma }}_z\otimes {\widehat{\sigma }}_x,$$

(31)

and the subsystem–ancilla coupling is chosen as Equation (29). From Equation (23), in Figure 9 we plot the heat current as a function of $\delta$ for these two cases. It can be seen that the heat current increases with the increase of $\delta$. When $T_1>T_2$, the heat current obtained from the decoupling approximation (cf. red dashed line in Figure 9) gradually deviates from the corresponding one obtained without the approximation (cf. red solid line in Figure 9) with the increase of $\delta$. We also investigate the corresponding trace distance ${\mathcal{D}}_{\mathcal{T}}$, and despite quantitative difference, its behaviors are similar to those of Figure 7, more specifically, as $\delta$ increases, the steady state predicted with the decoupling approximation gradually deviates from that without the approximation. So this approximation gets progressively worse as $\delta$ increases.

Moreover, in Figure 9, when $T_1<T_2$ we also consider the heat current without the decoupling approximation. It can be seen that in the case $T_1=10{\omega }_0$ and $T_2={\omega }_0$ (cf. red solid line in Figure 9), its corresponding forward heat current is greater than the backward heat currents ($T_2=10{\omega }_0$ and $T_1={\omega }_0$) (cf. blue solid line in Figure 9). So there is also a heat rectification, and the corresponding rectification factor is plotted in the inset of Figure 9, which decreases with the increase of $\delta$. It can be seen that an optimal rectification is almost realized with a high rectification factor $R>0.985$. Besides, under the same excitation-conserving system–environment interactions as those of Section 3, we also compare the forward and backward heat currents for this asymmetric system. The corresponding result shows that, despite an asymmetric heat conduction, there is no obvious thermal rectification. From numerical calculations, we find that as long as the system–environment interactions are chosen to be like those in Section 3, no matter how we introduce asymmetry into the system, there is no apparent thermal rectification effect. Therefore, we conclude that the rectification effect strongly depends on the form of the system–environment interaction. It might be the reason why most work [42,43,44,45] chose system–bath interactions like Equation (29) to achieve the thermal rectification for various kinds of asymmetric systems.

Next we move to the quantum correlation $\mathcal{D}$ between two subsystems. Without the decoupling approximation, in Figure 10 we plot $\mathcal{D}$ as a function of $\delta$ at steady state. It is clear that its behavior can be closely related to the failure of the approximation. We again confirm the same conclusion: the weaker the correlation between the two subsystems, the more justified the decoupling approximation and vice versa.

5. Conclusions

In this paper we have considered a system of two coupled qubits that each interact with their own heat bath. Based on the collision model, we mainly study whether the decoupling approximation (i.e., the inter-system coupling is ignored when modeling the system–environment interaction) is valid or not for describing the nonequilibrium dynamics. We have first investigated the stationary heat current and steady state for symmetric systems characterized by excitation-conserving system–ancilla interactions. In this situation, we have found that the decoupling approximation is still valid to describe either heat current or steady state even at strong inter-system couplings. Then, we have turned to the asymmetric systems in the case of non-excitation-conserving system–ancilla interactions. In this case we have found that the heat current or (steady state) predicted with the decoupling approximation gradually deviates from that without this approximation as $\delta$ increases, namely, this approximation gets progressively worse with the increase of $\delta$. In particular, without the decoupling approximation we have realized a perfect thermal rectification effect in this asymmetric systems. Through the analysis of different forms of interaction between system and environment, we have found that the thermal rectification effect is closely related to the form of system–ancilla interaction.

We have also considered the quantum discord between the two qubits without the approximation. For symmetric systems, we have found that quantum correlation between two subsystems first increases and then decreases as the inter-system coupling strength increases. The deviation of heat current (or steady state) from the decoupling approximation also increases and then decreases in this case. For asymmetric systems, the quantum correlation between two subsystems increases with the increase of inter-system coupling strength. The deviation of heat current (or steady state) from the decoupling approximation in this case increases with the increase of inter-system coupling strength. Although one would think that whether the decoupling approximation is valid or not should depend on inter-system coupling strength, we have found that it is closely related to the quantum correlations between the two subsystems. Specifically, the higher the quantum correlations means these two qubits behave more like a composite and must be treated as a whole, so the decoupling approximation gets worse, and vice versa.