Strong system-bath coupling induces negative differential thermal conductance and heat amplification in nonequilibrium two-qubits systems

TL;DR

This paper investigates quantum heat transfer in two-qubits systems, revealing that strong system-bath coupling can cause negative differential thermal conductance and significant heat amplification, advancing understanding of nonequilibrium quantum thermal devices.

Contribution

It demonstrates that strong qubit-bath coupling induces negative differential conductance and heat amplification in nonequilibrium two-qubits systems, unifying weak and strong coupling regimes.

Findings

Negative differential thermal conductance occurs with strong qubit-bath coupling.

Giant heat amplification factor observed in three-terminal setups.

Partial strong spin-boson interaction suffices for heat conductance effects.

Abstract

Quantum heat transfer is analyzed in nonequilibrium two-qubits systems by applying the nonequilibrium polaron-transformed Redfield equation combined with full counting statistics. Steady state heat currents with weak and strong qubit-bath couplings are clearly unified. Within the two-terminal setup, the negative differential thermal conductance is unraveled with strong qubit-bath coupling and finite qubit splitting energy. The partially strong spin-boson interaction is sufficient to show the negative differential thermal conductance. Based on the three-terminal setup, that two-qubits are asymmetrically coupled to three thermal baths, a giant heat amplification factor is observed with strong qubit-bath coupling. Moreover, the strong interaction of either the left or right spin-boson coupling is able to exhibit the apparent heat amplification effect.

Equations

$$\hat{H}=U{\hat{\sigma }}_z^L{\hat{\sigma }}_z^R+\sum _{v=L,R}{\hat{H}}_v,$$

$${\hat{H}}_v={\hat{H}}_b^v+\frac{{\epsilon }_v}{2}{\hat{\sigma }}_z^v+\frac{{\Delta }_v}{2}{\hat{\sigma }}_x^v+\sum _k(g_{kv}{\hat{b}}_{kv}^{†}+g_{kv}^{\ast }{\hat{b}}_{kv}){\hat{\sigma }}_z^v,$$

$${\hat{H}}_s^{\prime }=U{\hat{\sigma }}_z^L{\hat{\sigma }}_z^R+\sum _{v=L,R}(\frac{{\epsilon }_v}{2}{\hat{\sigma }}_z^v+\frac{{\eta }_v{\Delta }_v}{2}{\hat{\sigma }}_x^v),$$

$${\hat{V}}_v^{\prime }=\frac{{\Delta }_v}{2}[{\hat{\sigma }}_x^v(\cos \widehat{B_v}-{\eta }_v)+{\hat{\sigma }}_y^v\sin {\hat{B}}_v].$$

$$\frac{d{\hat{\rho }}_s(t)}{dt}$$

$${\Gamma }_{a,\pm }^v(\omega )={\int }_0^{\infty }d\tau C_a^v(\pm \tau )e^{i\omega \tau }.$$

$$C_x^v(\tau )$$

$$C_y^v(\tau )$$

$$Q_v(\tau )=4\sum _k|\frac{g_{kv}}{{\omega }_k}{|}^2[\cos {\omega }_k\tau (2n_k+1)-i\sin {\omega }_k\tau ],$$

$$\frac{d{\rho }_{nn}}{dt}$$

$$\mathcal{Z}({\chi }_v,t)$$

$$G({\chi }_v)=\underset{t\to \infty }{\lim }\frac{1}{t}\text{ln}\mathcal{Z}({\chi }_v,t).$$

$$\displaystyle~{}I^{(n)}=\frac{{\partial}^{n}G(\chi_{v})}{{\partial}(i\chi_{v})^{n}}\Big{|}_{\chi_{v}=0}.$$

$$\displaystyle I=\frac{{\partial}G(\chi_{v})}{{\partial}(i\chi_{v})}\Big{|}_{\chi_{v}=0}.$$

$${\hat{H}}_v(\chi )$$

$${\hat{H}}_{\chi }^{\prime }={\hat{H}}_s^{\prime }+\sum _{v=L,R}({\hat{H}}_b^v+{\hat{V}}_v(\chi )),$$

$${\hat{V}}_v(\chi )=\frac{{\Delta }_v}{2}[{\hat{\sigma }}_x^v(\cos \widehat{B_v^{\chi }}-{\eta }_v)+{\hat{\sigma }}_y^v\sin {\hat{B}}_v^{\chi }],$$

$$\frac{d{\hat{\rho }}_{\chi }(t)}{dt}$$

$${\Gamma }_{a,\pm }^v(\chi ,\omega )={\int }_0^{\infty }d\tau C_a^v(\chi ,\pm \tau )e^{i\omega \tau }.$$

$$I_{\mathrm{NIBA}}$$

$$\mathcal{A}$$

$${\kappa }_v^{ij}=(\frac{{\eta }_v{\Delta }_v}{2}{)}^2{\int }_{-\infty }^{\infty }d\tau e^{Q_v(\tau )}{e}^{i{E}_{ij}\tau }$$

$$I_{\mathrm{NIBA}}\approx \frac{(4U{)}^2\Delta T}{\mathcal{A}{T}_0^2}{\kappa }_L^{24}{\kappa }_R^{43}{\kappa }_L^{31}{\kappa }_R^{12}.$$

$${\beta }_{L_u}=|\partial I_{L_u}/\partial I_R|,u=h,c.$$

$${\beta }_{L_c}=|{\beta }_{L_h}+(-1{)}^{\theta }|,$$

$${\hat{H}}_s=\sum _v(\frac{{\epsilon }_v}{2}{\hat{\sigma }}_z^v+\frac{{\Delta }_v}{2}{\hat{\sigma }}_x^v)+U{\hat{\sigma }}_z^L{\hat{\sigma }}_z^R,$$

$${\hat{V}}_v=\sum _k(g_{kv}{\hat{b}}_{kv}^{†}+g_{kv}^{\ast }{\hat{b}}_{kv}){\hat{\sigma }}_z^v.$$

$${\hat{V}}_v^{\chi }=\sum _k(g_{kv}{e}^{i\frac{\chi {\omega }_k}{2}{\delta }_{v,R}}{\hat{b}}_{kv}^{†}+g_{kv}^{\ast }{e}^{-i\frac{\chi {\omega }_k}{2}{\delta }_{v,R}}{\hat{b}}_{kv}){\hat{\sigma }}_z^v.$$

$$\frac{d{\hat{\rho }}_{\chi }(t)}{dt}$$

$$\frac{d{\rho }_{n{n}^{\prime }}^{\chi }}{dt}$$

$$\frac{d{\rho }_{nn}^{\chi }}{dt}$$

Full text

I Introduction

Understanding mechanism and efficient manipulation of heat energy at nanoscale is of fundamental significance in the scientific community slepri2003pr ; gchen2005book ; ydubi2011rmp ; nbli2012rmp ; lee2013nature . One long-standing goal is to establish basic units of functional thermal devices for the implementation of heat transport and information processing bli2006prl ; lwang2007prl ; lwang2008prl ; kstan2012prl ; jren2015aip . Many theoretical models have been proposed for thermal operations, ranging from the thermal diode bli2004prl ; lfzhang2009prb ; afornieri2014apl ; rsanchez2015njp , thermal transistor bli2006prl ; tojanen2008prl ; truokola2011prb ; pba2014prl ; jhjiang2015prb ; kjoulain2016prl ; afornieri2016prb ; rsanchez2017prb ; gcraven2017prl ; yczhang2017apl ; yczhang2018epl , thermal memory lwang2007prl ; vkubytskyi2014prl ; cguarcello2018pra and even thermal computer lwang2008prl . And tremendous experiments were conducted to realize these novel thermal operations cwchang2006science ; fgiazotto2012nature ; syigen2014nl ; kjoulain2015apl ; bdutta2017prl ; dzhao2017nc . In particular for the thermal transistor, negative differential thermal conductance (NDTC) and heat amplification are considered as two key ingredients nbli2012rmp ; dhhe2009prb ; dhhe2010pre ; dhhe2014pre .

The NDTC was originally exploited in nonlinear phononic lattices coupled to two terminals by B. Li and his colleges bli2006prl , where the heat current shows suppression by the increase of the temperature bias between thermal baths. Based on the NDTC, the heat amplification was also exploited within the three-terminal phononic setup. Consequently, these concepts have been extensively investigated in various quantum thermal transistors. Particularly for the quantum dots, the capacitively Coulomb interaction allows the electronic fluctuation to transfer heat, and the asymmetric Coulomb blockade is crucial for the NDTC truokola2011prb ; yczhang2017apl . For the Josephson junctions, the coherent phase modulates heat transport between superconductors, and the analogous singularity-matching-peak effect originates the NDTC and heat amplification afornieri2016prb ; jren2013prb11 ; afornieri2017nnano . While for the near-field thermal transistor, a tiny change around the critical temperature of the insulator-metal transition material dramatically reduces the flux into the drain pba2014prl .

Recently, the quantum thermal transistor was theoretically unravelled in the nonequilibrium spin-boson models agleggett1987rmp ; dsegal2005prl ; uweiss2008book . Specifically, in the qubit-qutrit system, proper energy of the qutrit is crucial for the realization of thermal transistor effect bqguo2018pre . In resonant three-qubits systems, the NDTC and giant heat amplification factor were observed both in weak and strong coupling regimes kjoulain2016prl ; cwang2018pra , which are attributed to the three-terminal cooperative processes. Moreover, the strong spin-boson interaction is unravelled to be crucial to exhibit these far-from equilibrium features cwang2018pra . While for the nonequilibrium two-qubits system, in which the steady state heat transport has been preliminarily investigated lawu2011pra ; akato2015jcp , the NDTC and heat amplification is still lack of exploration, particularly in strong system-bath coupling regime. Based on the two-qubits system we raise questions: (i) can the NDTC be found with strong qubit-bath coupling within the two-bath setup? (ii) what is the role of the strong coupling on the NDTC? Moreover, considering the analysis of heat flow in the two-qubit system asymmetrically coupled to three-thermal baths zxman2016pre , can we observe heat amplification effect with strong qubit-bath coupling?

In the present paper, we investigate steady state quantum heat transfer in nonequilibrium two-qubits systems in Fig. 1, by applying the nonequilibrium polaron-transformed-Redfield equation (NE-PTRE) cklee2015jcp ; cwang2015sp ; dzxu2016fp combined with full counting statistics (FCS) llevitov1992jetp ; llevitov1996jmp ; mesposito2009rmp , detailed at Sec. II. The NE-PTRE can be reduced to the Redfield and noninteracting blip approximation (NIBA) schemes with weak and strong system-bath couplings, respectively. At Sec. III A, the NDTC is analyzed based on the two-terminal setup, and the mechanism of the NDTC is proposed in strong qubit-bath coupling limit. Moreover, the partially strong qubit-bath interaction is sufficient to exhibit the NDTC. At Sec. III B, the heat amplification is analyzed within the three-terminal setup, and the strong interaction of either the left or right qubit with the corresponding thermal bath can exhibit the dramatic heat amplification effect. Finally, we give a summary at Sec. IV.

II Model and method

In this section, we first describe the nonequilibrium two-qubits model. Then, by applying the polaron transformation we derive the nonequilibrium polaron-transformed Redfield equation to analyze the dynamics of qubits density matrix. Next, we obtain the expression of steady state heat flux by combining the NE-PTRE and full counting statistics. Under the thermodynamic bias of thermal baths, the steady state heat transfer can be exploited.

II.1 Nonequilibrium two-qubits model

The nonequilibrium two-qubits system consists of two coupled two-level-system (TLS), each interacting with an individual thermal bath in Fig. 1(a). The two-qubits model is described as

[TABLE]

where $U$ is the inter-qubit coupling strength, the Pauli operators are expressed as $\hat{\sigma}^{v}_{z}=|\uparrow_{v}{\rangle}{\langle}\uparrow_{v}|-|\downarrow_{v}{\rangle}{\langle}\downarrow_{v}|$ and $\hat{\sigma}^{v}_{x}=|\uparrow_{v}{\rangle}{\langle}\downarrow_{v}|+|\downarrow_{v}{\rangle}{\langle}\uparrow_{v}|$, with $|\uparrow_{v}{\rangle}~{}(|\downarrow_{v}{\rangle})$ the excited (ground) state of the $v$th qubit. $\hat{H}_{v}$ is the seminal spin-boson model

[TABLE]

where $\hat{H}^{v}_{b}=\sum_{k,v}\omega_{k}\hat{b}^{{\dagger}}_{kv}\hat{b}_{kv}$ stands for the $v$th noninteracting bosonic bath, with $\hat{b}^{{\dagger}}_{kv}~{}(\hat{b}_{kv})$ creating (annihilating) one phonon with energy $\omega_{k}$, $\varepsilon_{v}$ is the qubit splitting energy, $\Delta_{v}$ is the tunneling strength and $g_{kv}$ is the coupling strength between the $v$th qubit and the corresponding bath. The $v$th thermal bath is typically characterized by the spectral function $J_{v}(\omega)=4\pi\sum_{k}|g_{k,v}|^{2}\delta(\omega-\omega_{k})$. Here, we specify the spectral function as super-Ohmic form $J_{v}(\omega)=\pi\alpha_{v}{\omega^{3}}e^{-\omega/\omega_{c}}/{\omega^{2}_{c}}$, with $\alpha_{v}$ the qubit-bath coupling strength and $\omega_{c}$ the cutoff frequency of the bath, which has been extensively included in the quantum energy transfer dzxu2016fp ; anazir2009 ; sjang2013njp . It is well-known that the spin-boson model was initially proposed to analyze the dissipative dynamics of quantum systems agleggett1987rmp , and later extended to the nonequilibrium regime to analyze the quantum thermal transfer dsegal2005prl ; dsegal2006prb ; dsegal2008prl ; dsegal2011prb ; ksaito2013prl ; dsegal2014pre ; dsegal2015prl ; akato2015jcp ; jliu2017pre ; zqiang2017prb .

From Eq. (2), it is known that qubits linearly interact with thermal bath. Hence, to capture the multi-phonons effect in thermal transfer processes, we apply a canonical transformation (”polaron transformation”) to obtain the modified Hamiltonian $\hat{H}^{\prime}=\hat{U}^{{\dagger}}\hat{H}\hat{U}$ rsilbey1984jcp ; rharris1985jcp ; hzheng2004epjb , where the unitary operator is $\hat{U}=\exp{(i\sum_{v}\hat{B}_{v}\hat{\sigma}^{v}_{z}/2)}$ , with the collective phononic momentum $\hat{B}_{v}=2i\sum_{k}(g_{kv}/\omega_{k}\hat{b}^{{\dagger}}_{kv}-g^{*}_{kv}/\omega_{k}\hat{b}_{kv})$. Then, the modified system Hamiltonian is given by $\hat{H}^{\prime}=\hat{H}^{\prime}_{s}+\sum_{v}(\hat{H}^{v}_{b}+\hat{V}^{\prime}_{v})$, where the transformed qubits Hamiltonian is

[TABLE]

where the renormalization factor $\eta_{v}={\langle}\cos\hat{B}_{v}{\rangle}$ is specified as $\eta_{v}=\exp[-2\sum_{k}|g_{kv}/\omega_{k}|^{2}(2n_{k}+1)]$. Moreover, the modified qubit-bath interaction is given by

[TABLE]

By expanding the interaction terms at Eq. (4) as $\hat{\sigma}^{v}_{x}(\cos\hat{B_{v}}-\eta_{v})=\hat{\sigma}^{v}_{x}[(1-\eta_{v})+\sum_{n=1}(\hat{B_{v}})^{2n}/(2n)!]$ and $\hat{\sigma}^{v}_{y}\sin\hat{B}_{v}=\hat{\sigma}^{v}_{y}[\sum_{n=0}(\hat{B_{v}})^{2n+1}/(2n+1)!]$, it is easily found that multiple phonon effect is now explicitly included in heat transfer processes.

II.2 Nonequilibrium polaron-transformed Redfield equation

We apply the nonequilibrium polaron-transformed Redfield equation, one type of quantum master equation, to study the dynamics of qubits density matrix. Traditionally, the qubit-bath interaction at Eq. (2) is directly perturbed to obtain the Redfield equation, which is proper for weak qubit-bath coupling dsegal2005prl ; dsegal2006prb . While for NE-PTRE, it has been intensively applied in quantum heat transfer and heat engines cklee2015jcp ; cwang2015sp ; dzxu2016fp . By perturbing the modified qubit-bath interaction at Eq. (4), we are able to analyze dynamical behaviors in the broad coupling region. Specifically, based on the Born approximation, the total density matrix in the polaron framework is decomposed as $\hat{\rho}_{tot}(t){\approx}\hat{\rho}_{s}(t){\otimes}\hat{\rho_{b}}$, where $\hat{\rho}_{s}(t)$ is the reduced density operator of qubits and $\hat{\rho_{b}}=\exp{(-\sum_{v}\hat{H}^{v}_{b}/k_{B}T_{v})}/\textrm{Tr}_{b}\{e^{-\sum_{v}\hat{H}^{v}_{b}/k_{B}T_{v}}\}$ is the equilibrium distribution operator of thermal baths. Moreover, by considering the Markovian approximation, the NE-PTRE is obtained as

[TABLE]

where the projecting operators are $\hat{\sigma}^{v}_{a}(-\tau)=\sum_{\omega}\hat{P}^{v}_{a}(\omega)e^{i\omega\tau}$, transition rates are

[TABLE]

and the correlation functions are given by

[TABLE]

with the propagating phase of phonons

[TABLE]

and the Bose-Einstein distribution function $n_{kv}=1/[\exp{(\omega_{k}/k_{B}T_{v})}-1]$. The transition rate $\Gamma^{v}_{y,\pm}(\omega)$ describes that odd phonons are involved in the heat exchange processes between the $v$th qubit and the corresponding bath with transition energy $\omega$ at Eq. (6). For the lowest order of correlation function $C^{v}_{y}(\pm\tau)$ as $(\frac{\eta_{v}\Delta_{v}}{2})^{2}Q_{v}(\pm\tau)$, the rates are reduced to $\textrm{Re}[\Gamma^{v}_{y,-}(\omega)]=(\frac{\eta_{v}\Delta_{v}}{2\omega})^{2}J_{v}(\omega)n_{v}(\omega)$ and $\textrm{Re}[\Gamma^{v}_{y,+}(\omega)]=(\frac{\eta_{v}\Delta_{v}}{2\omega})^{2}J_{v}(\omega)(1+n_{v}(\omega))$, with the spectral function of the $v$th bath as $J_{v}(\omega)=4\pi\sum_{k}|g_{k,v}|^{2}\delta(\omega-\omega_{k})$ and $n_{v}(\omega)=1/[e^{\omega/k_{B}T_{v}}-1]$, which exhibits sequential excitation and relaxation of the $v$th qubit by absorbing and emitting one phonon, accordingly. While the transition rate $\Gamma^{v}_{x,\pm}(\omega)$ shows even phonons involved transfer processes. The lowest order of correlation function $C^{v}_{x}(\pm\tau)$ is given by $(\frac{\eta_{v}\Delta_{v}}{2})^{2}Q^{2}_{v}(\pm\tau)$. Hence, the reduced rates are obtained as $\textrm{Re}[\Gamma^{v}_{x,\pm}(\omega)]=(\frac{\eta_{v}\Delta_{v}}{2})^{2}\int^{\infty}_{0}\frac{d{\omega_{1}}}{4\pi}Q_{v}(\pm\omega_{1})Q_{v}(\pm\omega\mp\omega_{1})$, with $Q_{v}(\omega_{1})=\int{d\tau}e^{i\omega_{1}\tau}Q_{v}(\tau)$. For $\textrm{Re}[\Gamma^{v}_{x,\pm}(\omega)]$, as the $v$th qubit is excited(relaxed) with energy gap $\omega$, it simultaneously absorbs(emits) two phonons with energy $\omega_{1}$ and $\omega-\omega_{1}$ from(to) the $v$th bath.

In previous works of nonequilibrium spin-boson model, the NE-PTRE was found to unify steady state behaviors in the weak and strong qubit-bath coupling limit, respectively cwang2015sp ; dzxu2016fp ; cwang2017pra . Here, we note that in the two-qubits system the NE-PTRE can also find such correspondence. Specifically, in the weak coupling regime ($|g_{k,v}/\omega_{k}|^{2}{\approx}0$), the renormalization factor becomes $\eta_{v}{\approx}1$. Hence, the qubits Hamiltonian is simplified as $\hat{H}^{\prime}_{s}=U\hat{\sigma}^{1}_{z}\hat{\sigma}^{2}_{z}+\sum_{v}(\frac{\varepsilon_{v}}{2}\hat{\sigma}^{v}_{z}+\frac{\Delta_{v}}{2}\hat{\sigma}^{v}_{x})$. And the lowest order of transformed qubit-bath interaction $\hat{V}^{\prime}_{v}{\approx}\hat{\sigma}^{v}_{y}\hat{B}_{v}$ dominates the dynamics, which results in the correlation functions $C^{v}_{x}(\tau){\approx}0$ and $C^{v}_{y}(\tau){\approx}(\frac{\Delta_{v}}{2})^{2}Q_{v}(\tau)$. Under the basis $\{|\phi_{n}{\rangle}\}$ with $\hat{H}^{\prime}_{s}|\phi_{n}{\rangle}=E_{n}|\phi_{n}{\rangle}$, the dynamical equation at Eq. (II.2) after long-time evolution can be reduced to

[TABLE]

with the density matrix element $\rho_{nn}={\langle}\phi_{n}|\hat{\rho}_{s}|\phi_{n}{\rangle}$, the energy gap $E_{nm}=E_{n}-E_{m}$, $n_{v}(E_{nm})=1/[\exp{(E_{nm}/k_{B}T_{v})}-1]$ and $\sigma^{nm}_{y,v}={\langle}\phi_{n}|\hat{\sigma}^{v}_{y}|\phi_{m}{\rangle}$. Moreover, from the commutating relationship ${\langle}\phi_{n}|[\hat{\sigma}^{v}_{z},\hat{H}^{\prime}_{s}]|\phi_{m}{\rangle}=i\Delta_{v}\sigma^{nm}_{y,v}=E_{nm}\sigma^{nm}_{z,v}$, we conclude that Eq. (II.2) is equivalent with the Redfield equation at Eq. (28) (see appendix A for the detail).

While in the strong qubit-bath interaction regime ($|g_{k,v}/\omega_{k}|^{2}{\gg}1$), the renormalization factor becomes $\eta_{v}{\approx}0$. Hence the qubits Hamiltonian is simplified as $\hat{H}^{\prime}_{s}=U\hat{\sigma}^{1}_{z}\hat{\sigma}^{2}_{z}+\sum_{v}\frac{\varepsilon_{v}}{2}\hat{\sigma}^{v}_{z}$, and the qubit-bath interaction is give by $\hat{V}^{\prime}_{v}=\frac{\Delta_{v}}{2}[\hat{\sigma}^{v}_{x}\cos{B}_{v}+\hat{\sigma}^{v}_{y}\sin{\hat{B}_{v}}]$. The correlation functions are reduced to $C^{v}_{x}(\tau)=C^{v}_{y}(\tau)=\Delta^{2}_{v}\eta^{2}_{v}e^{Q_{v}(\tau)}/8$. It should be noted that though $\eta_{v}{\approx}0$, $\eta^{2}_{v}e^{Q_{v}(\tau)}$ will be kept as a finite value. Consequently, NE-PTRE at Eq. (II.2) is reduced to the nonequilibrium NIBA limit at Eq. (B) (see appendix B for the detail).

II.3 Steady state heat current

We combine NE-PTRE with FCS to count the heat flow into thermal baths mesposito2009rmp ; cwang2015sp . FCS is considered as a powerful method to investigate the full information of current fluctuations, which was initially proposed to analyze the charge fluctuations llevitov1992jetp ; llevitov1996jmp . For the nonequilibrium spin-boson system, FCS has been extensively applied to study quantum heat transport at steady state or driven by the geometric phase jren2010prl ; jren2013prb ; cuchiyama2014pre ; cwang2017pra ; jliu2018jcp . Generally, by employing a two-time measurement protocol during the time interval $[0,t]$, the generating function to characterize heat transfer is expressed as hmfriedman2018njp $\mathcal{Z}(\chi_{v},t)=\textrm{Tr}\{e^{i\chi_{v}\hat{H}^{v}_{b}(0)}e^{-i\chi_{v}\hat{H}^{v}_{b}(t)}\hat{\rho}_{tot}(0)\}$, where $\chi_{v}$ is the counting field parameter to count energy into the $v$th bath, $\hat{\rho}_{tot}(0)$ is the initial density operator of the whole system, and the operator $\hat{A}(t)$ in the Heisenberg representation is $\hat{A}(t)=e^{i\hat{H}t}\hat{A}e^{-i\hat{H}t}$. Moreover, by assuming $[\hat{\rho}_{tot}(0),\hat{H}^{v}_{b}]=0$, the generating function can be re-expressed as

[TABLE]

where the generalized propagators are $\hat{M}_{-\chi_{v}}(t)=e^{-i\hat{H}_{-\chi_{v}}t}$ and $\hat{M}^{{\dagger}}_{\chi_{v}}(t)=e^{i\hat{H}_{\chi_{v}}t}$, with $\hat{H}_{\chi_{v}}=e^{\frac{i\chi_{v}}{2}\hat{H}^{v}_{b}}\hat{H}e^{-\frac{i\chi_{v}}{2}\hat{H}^{v}_{b}}$. Then, the cumulant generating function at steady state can be obtained as

[TABLE]

Consequently, the $n$th-order current cumulant is given by

[TABLE]

In particular, the lowest-order cumulant is the heat current

[TABLE]

Therefore, once we obtain the dynamics of the counting-field dependent density matrix of the qubits system, we are able to obtain full information of current fluctuations.

Then, we apply FCS to study the heat transfer in two-qubit system. The counting-field dependent Hamiltonian $\hat{H}(\chi)=e^{\frac{i\chi}{2}\hat{H}_{b,R}}\hat{H}e^{-\frac{i\chi}{2}\hat{H}_{b,R}}$ is given by $\hat{H}_{{\chi}}=U\hat{\sigma}^{L}_{z}\hat{\sigma}^{R}_{z}+\sum_{v=L,R}\hat{H}_{v}(\chi)$, where the sub-Hamiltonian is

[TABLE]

with $\hat{H}^{v}_{b}=\sum_{k}\omega_{k}\hat{b}^{{\dagger}}_{k,v}\hat{b}_{k,v}$. Furthermore, we apply the modified canonical operator $\hat{U}_{\chi}=\exp[\frac{i}{2}(\hat{B}_{L}\hat{\sigma}^{L}_{z}+\hat{B}_{R}(\chi)\hat{\sigma}^{R}_{z})]$ to $\hat{H}_{\chi}$ as $\hat{H}^{\prime}_{\chi}=\hat{U}^{{\dagger}}_{\chi}\hat{H}_{\chi}\hat{U}_{\chi}$, resulting in

[TABLE]

with

[TABLE]

and $\hat{B}^{\chi}_{v}=2i\sum_{k}(\frac{g_{kv}}{\omega_{k}}e^{\frac{i\chi}{2}\omega_{k}\delta_{v,R}}\hat{b}^{{\dagger}}_{kv}-H.c.)$. Following a similar procedure, the total density matrix is decomposed as $\hat{\rho}^{\chi}_{tot}(t)=\hat{\rho}_{\chi}(t){\otimes}\hat{\rho}_{b}$. The counting-field dependent NE-PTRE is obtained as

[TABLE]

where the modified transition rates are

[TABLE]

and the correlation functions are $C^{v}_{a}(\chi,\pm\tau)=C^{v}_{a}(\pm\tau-\chi\delta_{v,R})$. As the counting field parameter vanishes, the dynamical equation at Eq. (II.3) is reduced to the standard version at Eq. (II.2).

Now, based on the definition of generating function at Eq. (II.3) and current fluctuations at Eq. (12), it is ready to analyze the steady state heat transfer by applying the dynamics of qubits density matrix at Eq. (II.3). In the following, we simplify system parameters as $\Delta_{v}=\Delta$, $\varepsilon_{v}=\varepsilon$. Here, the first question we should answer is: does the heat current derived from the counting-field dependent NE-PTRE unifies the results based on the Redfield and nonequilibrium NIBA schemes? In Fig. 2, it clearly demonstrates that both at resonance ($\varepsilon=0$) and off-resonance (e.g., $\varepsilon=1$), heat currents perfectly bridge the weak and strong coupling results. Therefore, it is safe to apply the NE-PTRE in a wide qubit-bath coupling regime to analyze nonequilibrium steady state behaviors.

III Results and discussions

III.1 Negative differential thermal conductance

Negative differential thermal conductance (NDTC), as analogue of the negative differential conductance in electronics, traditionally describes the process that heat current is suppressed by increasing the temperature bias in a two-terminal setup bli2006prl . It has been extensively applied to investigate quantum thermal transport in phononics nbli2012rmp . However, the NDTC is less exploited for nonequilibrium two-qubits systems.

Here, we first analyze the NDTC in a biased two-qubits model by tuning qubit-bath coupling within the NE-PTRE in Fig. 3(a). For the weak coupling case ($e.g.,\alpha=0.05$), It is shown that the normalized heat current $I/I_{max}(\alpha)$ becomes monotonically enhanced by tuning the temperature bias $\Delta{T}$. As the coupling strength is enlarged, the NDTC appears and finally becomes significant ($e.g.,\alpha=5$). Hence, we conclude that the NDTC exists in the strongly-coupled two-spin-boson model. To give a comprehensive picture of the NDTC under the modulation of the temperature bias and the coupling strength, we plot Fig. 3(b). It is found that at the strong coupling regime ($\alpha{\gtrsim}0.5$) the $I/I_{max}(\alpha)$ shows the turnover behavior, corresponding to the NDTC.

In previous works with strong qubit-bath coupling, it is known that the single qubit is unable to exhibit the NDTC cwang2015sp . While for the three-qubits case, the NDTC is clearly shown cwang2018pra . Therefore, considering the result obtained in the present paper with the two-qubits system, we conclude that multi-qubits and strong qubit-bath interaction are necessary conditions to generate the NDTC. It should be noted that the current conclusion is limited to such (qubit-qubit and qubit-bath) interactions of the Hamiltonian at Eq. (1) and Eqs. (1-2) in Ref. cwang2018pra . While for other kinds of quantum spin systems, the NDTC can also be unravelled with weak qubit-bath coupling, e.g., quantum thermal transistor kjoulain2016prl and spin seebeck transistor jren2013prb ; jren2013prb2 ; jren2013prb3 .

III.1.1 The mechanism of NDTC

Next, we try to provide some insight into the underlying mechanism of the NDTC in strong qubit-bath coupling limit. It is known that with strong coupling the nonequilbrium NIBA can at least qualitatively describe heat transfer processes. Under the condition $\varepsilon_{v}=\varepsilon$, the steady state heat current can be approximately expressed as (see Eq. (38) at appendix B)

[TABLE]

where the coefficient $\mathcal{A}$ is

[TABLE]

which is invariant by exchanging the subindex $L$ with $R$. $\kappa^{2i}_{v}=\kappa^{3i}_{v}$, $\kappa^{i2}_{v}=\kappa^{i3}_{v}$ and the transition rate is given by

[TABLE]

with the energy gap is $E_{ij}=E_{i}-E_{j}$ with $(U\hat{\sigma}^{L}_{z}\hat{\sigma}^{R}_{z}+\sum_{v}\frac{\varepsilon_{v}}{2}\hat{\sigma}^{v}_{z})|i{\rangle}=E_{i}|i{\rangle}~{}(i=1,2,3,4)$. From Eq. (19), it is clearly shown that the heat current is contributed by two loop currents ${4U}\kappa^{24}_{L}\kappa^{43}_{R}\kappa^{31}_{L}\kappa^{12}_{R}/{\mathcal{A}}$ and ${4U}\kappa^{21}_{R}\kappa^{13}_{L}\kappa^{34}_{R}\kappa^{42}_{L}/{\mathcal{A}}$.

Then, we analyze the heat current in two limiting temperature bias regimes to exploit the NDTC. For the low bias ${\Delta}T{\ll}T_{0}$, the heat current is approximately expressed as

[TABLE]

Moreover, the term $\kappa^{24}_{L}\kappa^{43}_{R}\kappa^{31}_{L}\kappa^{12}_{R}/{\mathcal{A}}$ is nearly constant, shown with blue-dark dashed line with circles in Fig. 4(c). Hence, $I_{\mathrm{NIBA}}$ is proportional to ${\Delta}T$, which is exhibited with black solid line. While for the large temperature bias ${\Delta}T{\approx}2T_{0}$, i.e. $T_{L}{\approx}2T_{0}$ and $T_{R}{\approx}0$, The transition from $|\downarrow\downarrow{\rangle}$ to $|\downarrow\uparrow{\rangle}$ assisted by the right bath is significantly suppressed($\kappa^{43}_{R}{\approx}0$), which is schematically demonstrated in Fig. 4(a) and also exhibited in Fig. 4(b). It is mainly attributed to the fact that phonons are difficult to excite in low temperature limit($T_{R}{\approx}0$). Hence, the loop current ${4U}\kappa^{24}_{L}\kappa^{43}_{R}\kappa^{31}_{L}\kappa^{12}_{R}/{\mathcal{A}}$ becomes negligible, shown in Fig. 4(c). Similarly, the loop current $\frac{4U}{\mathcal{A}}\kappa^{21}_{R}\kappa^{13}_{L}\kappa^{34}_{R}\kappa^{42}_{L}$ is also blocked, due to the monotonic suppression of the transition rate $\kappa^{21}_{R}$(not shown here). Therefore, we conclude that the heat current is strongly suppressed by the large temperature bias. Combining the linear increase of the current in the linear response regime and dramatic suppression at large temperature bias, the heat current should exhibit the turnover feature, which is the signature of the NDTC. It should be noted that the analysis of the origin of the NDTC is approximate based on the nonequilibrium NIBA.

III.1.2 Partially strong qubit-bath interaction exhibits NDTC

Next, we raise the question that whether the strong coupling between the $v$th qubit and the corresponding thermal bath is sufficient to induce the negative differential thermal conductance? By setting the weak left qubit-bath coupling (e.g., $\alpha_{L}=0.05$) in Fig. 5(a), we investigate the influence of the right qubit-bath interaction in the steady state heat current. It is interesting to find that in strong coupling regime of $\alpha_{R}$, the NDTC also appears. Thus, the strong coupling between the right qubit and bath is sufficient to exhibit the NDTC. To give a comprehensive picture of the NDTC, we plot Fig. 5(b). It is found that as $\alpha_{R}{\gtrsim}1.2$, the NDTC can be observed by tuning the temperature bias. Hence, partially strong qubit-bath interaction is sufficient to exhibit the NDTC.

III.2 Heat amplification

Heat amplification, the key ingredient of quantum thermal transistor, shows a power to dramatically amplify heat flow via a tiny change of one terminal current within a three-terminal setup bli2006prl . It has been exploited in nonequilibrium three-qubits systems with both weak and strong qubit-bath coupling cases kjoulain2016prl ; cwang2018pra . In Fig. 1(b), the amplification factor can be expressed as

[TABLE]

Considering the energy conservation $\sum_{v=L_{h},L_{c},R}I_{v}=0$, the factor relationship with the $L_{h}$ and $L_{c}$ terminals can be established as

[TABLE]

where the phase $\theta=0$ as ${\partial}I_{L_{u}}/{\partial}I_{R}>0$ and $\theta=1$ as ${\partial}I_{L_{u}}/{\partial}I_{R}<0$. Traditionally, the thermal transistor is able to work as $\beta_{L_{u}}>1$. Here, we focus on the amplification effect in a two-qubits system, asymmetrically coupled to three thermal baths (see the derivation of the dynamical equation at appendix C).

We first investigate the influence of the qubit-bath interaction in the heat amplification factor $\beta_{L_{h}}$ by tuning $T_{R}$ in Fig. 6(a). In the weak coupling regime(e.g., $\alpha=0.05$), $\beta_{L_{h}}$ is below the unit, which implies the breakdown of heat amplification. While in strong coupling regime (e.g., $\alpha=5$), it is interesting to observe a giant amplification factor in comparatively low temperature regime around $T_{R}{\approx}0.5$, and is dramatically suppressed as qubits is driven away from this critical temperature regime. Moreover, in Fig. 6(b) the heat current into the $R$th bath shows nonmonotonic behavior with strong coupling by modulating $T_{R}$. From the definition of $\beta_{L_{h}}$ at Eq. (22), it is understandable the amplification factor becomes divergent at the turnover point of $I_{R}$ (${\partial}I_{R}/{\partial}T_{R}=0$). Therefore, we conclude that strong qubit-bath coupling is able to exhibit the significant heat amplification effect.

Next, we try to find out separately that whether the strong coupling of left two thermal baths or the right bath exhibits the generation of the heat amplification. In Fig. 7(a) with weak left qubit-bath couplings, it is interesting to find that the heat amplification occurs with strong coupling between the right qubit and corresponding bath. While in Fig. 7(b) with strong left qubit-bath couplings, a giant amplification factor is shown even at weak right qubit-bath coupling case ($e.g.,\alpha_{R}=0.1$). Hence, we conclude that strong coupling either from the left qubit-baths or right counterpart results in the heat amplification effect. It should be noted these results are different from the previous work of nonequilibrium three-qubits system cwang2018pra , where the strong coupling from the thermal baths with moderate temperature (corresponding to $T_{R}$ in Fig. 1(b)) is necessary to exhibit the amplification effect.

IV Conclusion

To summarize, we study the quantum thermal transfer in the nonequilibrium two-qubits system, by applying the nonequilibrium polaron-transformed Redfield equation combined with full counting statistics. The heat current is unified in a broad coupling regime, which reduces to the Redfield and NIBA schemes in weak and strong interaction limits. The negative differential thermal conductance is investigated within the two-terminal setup in Fig. 1(a), and the NDTC occurs with strong qubit-bath coupling and large temperature bias, which clearly answers the first question raised in the introduction. The underlying mechanism is exploited based on the expression of the heat current at Eq. (19). Moreover, by setting the weak left-qubit coupling, it is interesting to find that the NDTC sustains in strong right qubit-bath coupling regime. This clearly demonstrates that the partially strong qubit-bath interaction is sufficient to show the appearance of the NDTC. Next, we analyze the heat amplification effect within the three-terminal setup in Fig. 1(b), where two qubits are asymmetrically coupled to three thermal baths. A giant heat amplification factor is clearly observed in strong qubit-bath coupling regime with comparatively low temperature of the $R$th bath, which is the answer for the second question in the introduction. It is tightly related with the turnover behavior of the heat current into the $R$th bath. Moreover, it is discovered that either the strong interaction of the left qubit-bath interactions or the right counterpart can exhibit such significant amplification effect. We hope these findings can provide some insight into designing the functional thermal transistor in coupled-qubits systems.

V Acknowledgement

H.L and C.W. are supported by the National Natural Science Foundation of China under Grant No. 11704093. L.Q.W and J.R. acknowledge the support by the National Natural Science Foundation of China (No. 11775159), Natural Science Foundation of Shanghai (No. 18ZR1442800), and the National Youth 1000 Talents Program in China.

Appendix A Steady state heat current within Redfield

The two-coupled spin-boson system is expressed as $\hat{H}=\hat{H}_{s}+\sum_{v}(\hat{H}^{v}_{b}+\hat{V}_{v})$, where the qubits Hamiltonian is given by

[TABLE]

the $v$th thermal bath is $\hat{H}^{v}_{b}=\sum_{k}\omega_{k}\hat{b}^{{\dagger}}_{k,v}\hat{b}_{k,v}$, and the system-bath interaction is

[TABLE]

To count the energy flow into the right bath including full counting statistics, the total Hamiltonian is changed to $\hat{H}(\chi)=e^{i\chi\hat{H}^{R}_{b}/2}\hat{H}e^{-i\chi\hat{H}^{R}_{b}/2}=\hat{H}_{s}+\sum_{v}(\hat{H}^{v}_{b}+\hat{V}^{\chi}_{v})$, where

[TABLE]

Considering the weak qubit-bath interaction, we perturb the system-bath interaction at Eq. (25) up to the second order. Then, based on the Born-Markov approximation the Redfield equation is given by

[TABLE]

where the commutating relation is $[\hat{A}_{\chi},\hat{B}_{\chi}]_{\chi}=\hat{A}_{\chi}\hat{B}_{\chi}-\hat{B}_{\chi}\hat{A}_{-\chi}$. In the eigen-basis $\{|\phi_{n}{\rangle}\}$ with $\hat{H}_{s}|\phi_{n}{\rangle}=E_{n}|\phi_{n}{\rangle}$, the dynamical equation of the system density matrix element is given by

[TABLE]

Moreover, from Fig. 8 it is known that off-diagonal elements becomes negligible at steady state. Hence, the Redfield equation after the long time evolution is reduced to

[TABLE]

Finally, the steady state heat flux is given by

[TABLE]

Then, we briefly discuss the absence of the NDTC within the Redfield scheme in Fig. 9. We first define the transition current from the population in the eigenspace $P_{n}$ to $P_{m}$ as $J_{n{\rightarrow}m}=\frac{E_{mn}}{2}J_{R}(E_{mn})(1+n_{R}(E_{mn}))\sigma^{R}_{nm}\sigma^{R}_{mn}P_{m}$, and the net transition current $J_{nm}=J_{n{\rightarrow}m}+J_{m{\rightarrow}n}$. With the weak qubit-qubit coupling(e.g., $U=0.1$ in Fig. 9), it is found that $J_{14}$ and $J_{23}$ are negligible. Thus, the heat current can be simplified as

[TABLE]

By increasing the temperature bias, the contributing net currents at Eq. (30) show monotonic increase, resulting in the enhancement of the $J_{\mathrm{Redfield}}$. Hence, $J_{\mathrm{Redfield}}$ will not exhibit the NDTC.

Appendix B Steady state heat current within NIBA

Following the similar treatment in Sec. II (B), after the polaron transformation the system Hamiltonian is given by

[TABLE]

and the system-bath interaction is

[TABLE]

with the collective phonon momentum $\hat{B}_{v}=2i\sum_{k}(g_{k,v}/\omega_{k}\hat{b}^{{\dagger}}_{k,v}-H.c.)$. By applying the full counting statistics, the counting field dependent qubit-bath interaction is changed to $V^{\chi}_{v}=\hat{\sigma}^{v}_{x}\cos\hat{B}^{\chi}_{v}+\hat{\sigma}^{v}_{y}\sin\hat{B}^{\chi}_{v}$, with $\hat{B}^{\chi}_{v}=2i\sum_{k}(\frac{g_{k,v}}{\omega_{k}}e^{\frac{i\chi}{2}\omega_{k}\delta_{v,R}}\hat{b}^{{\dagger}}_{k,v}-H.c.)$. Consequently, the quantum master equation is given by

[TABLE]

where the transition correlation function is

[TABLE]

Specifically, the population dynamics in local basis $\{|\phi_{1}{\rangle}=|\uparrow\uparrow{\rangle},|\phi_{2}{\rangle}=|\uparrow\downarrow{\rangle},|\phi_{3}{\rangle}=|\downarrow\uparrow{\rangle},|\phi_{4}{\rangle}=|\downarrow\downarrow{\rangle},\}$ are

[TABLE]

where $P^{\chi}_{i}={\langle}\phi_{i}|\hat{\rho}_{\chi}|\phi_{i}{\rangle}$ the counting field dependent transition rates are $\kappa^{ij}_{v}({\chi})=e^{iE_{ij}\chi}\kappa^{ij}_{v}$, and the standard transition rates are

[TABLE]

Consequently, the heat current is expressed as

[TABLE]

Under the condition $\varepsilon_{v}=\varepsilon$, the steady state populations in absence of the counting field parameter are given by

[TABLE]

where the coefficient $\mathcal{A}$ is

[TABLE]

And the current is simplified as

[TABLE]

From Fig. 3, it is known that the NDTC occurs with finite qubit splitting energy with strong qubit-bath coupling. Here, we analyze the influence of the splitting energy on the steady state heat current in Fig. 10(a). At resonance ($\varepsilon=0$), the Normalized heat current shows monotonic enhancement by increasing temperature bias (${\Delta}T$) without the NDTC. From the Eq. (38), the heat current is contributed by two loop currents $\frac{4U}{\mathcal{A}}\kappa^{24}_{L}\kappa^{43}_{R}\kappa^{31}_{L}\kappa^{12}_{R}$ and $\frac{4U}{\mathcal{A}}\kappa^{21}_{R}\kappa^{13}_{L}\kappa^{34}_{R}\kappa^{42}_{L}$. In Fig. 10, it is found that the second component $\frac{4U}{\mathcal{A}}\kappa^{21}_{R}\kappa^{13}_{L}\kappa^{34}_{R}\kappa^{42}_{L}$ is strongly suppressed by increasing ${\Delta}T$, which enhances the difference between two loop currents. Hence, $I_{\mathrm{NIBA}}$ exhibits monotonic behavior. Next by turning on $\varepsilon$ (e.g., $\varepsilon=0.5,0.8$), it is found the NDTC finally occurs. Therefore, we conclude that qubit splitting energy is necessary to show the NDTC in the nonequilibrium two-qubit system.

Then, we study the effect of quibt-qubit interaction on the heat current in Fig. 10(b). It is found that by increasing the repulsion strength $U$, the behavior of NDTC is suppressed accordingly, and completely vanished with strong repulsion strength(e.g., $U=0.8$). With strong qubit-qubit repulsion, the current $\frac{4U}{\mathcal{A}}\kappa^{24}_{L}\kappa^{43}_{R}\kappa^{31}_{L}\kappa^{12}_{R}$ is apparently enhanced by increasing ${\Delta}T$, whereas $\frac{4U}{\mathcal{A}}\kappa^{21}_{R}\kappa^{13}_{L}\kappa^{34}_{R}\kappa^{42}_{L}$ shows decrease in Fig. 10(d). It finally contributes to the monotonic increase of $I_{\mathrm{NIBA}}$.

Appendix C Heat current in two-qubits asymmetrically coupled to three thermal baths

The two qubits system asymmetrically coupled to three thermal baths, which is shown in Fig. 1(b), is modeled as

[TABLE]

The system Hamiltonian is given by

[TABLE]

The $v$th thermal bath is described as $\hat{H}^{v}_{b}=\sum_{k}\omega_{k}\hat{b}^{{\dagger}}_{k,v}\hat{b}_{k,v}$. The $L_{u}$th qubit-bath interaction is given by

[TABLE]

and the $R$th qubit-bath interaction is

[TABLE]

Then, we apply the full counting statistics to count the heat flow into the $v$th thermal bath as

[TABLE]

where the counting parameter set is $\{\chi\}=(\chi_{L_{h}},\chi_{L_{c}},\chi_{R})$, $\hat{V}_{L_{u}}(\chi_{L_{u}})=\hat{\sigma}^{L}_{z}\sum_{k}(g_{k,L_{u}}e^{\frac{i\omega_{k}\chi_{L_{u}}}{2}}\hat{b}^{{\dagger}}_{k,L_{u}}+H.c.)$ and $\hat{V}_{R}(\chi_{R})=\hat{\sigma}^{R}_{z}\sum_{k}(g_{k,R}e^{\frac{i\omega_{k}\chi_{R}}{2}}\hat{b}^{{\dagger}}_{k,R}+H.c.)$.

Similar to the treatment at Eq. (15), we apply the canonical transformation to obtain the modified Hamiltonian as

[TABLE]

where the unitary operator is given by $\hat{U}_{\{\chi\}}=\exp\{i[(\sum_{u}\hat{B}_{L_{u}}(\chi_{L_{u}}))\hat{\sigma}^{L}_{z}+\hat{B}_{R}(\chi_{R})\hat{\sigma}^{R}_{z}]$, with the collective phononic momentum $\hat{B}_{v}(\chi_{v})=2i\sum_{k}(g_{kv}e^{\frac{i\omega_{k}\chi_{v}}{2}}/\omega_{k}\hat{b}^{{\dagger}}_{kv}-H.c.)$. The modified Hamiltonian is specified as $\hat{H}^{\prime}(\{\chi\})=\hat{H}^{\prime}_{s}+\hat{V}^{\prime}(\{\chi\})+\sum_{v}\hat{H}^{v}_{b}$, where the transformed qubits Hamiltonian is

[TABLE]

with the renormalization factors

[TABLE]

The system-bath interaction is given by

[TABLE]

Then, based on the Born-Markov approximation and perturbing the system-bath interaction up to the second order, we obtain the quantum master equation as

[TABLE]

where the transition rates are

[TABLE]

and the counting-field dependent correlation functions are

[TABLE]

where the correlation phase is

[TABLE]

Finally, the cumulant generating function after a long time evolution is given by $\mathcal{Z}(\{\chi\})=\lim_{t{\rightarrow}\infty}\frac{1}{t}\ln\textrm{Tr}\{\hat{\rho}_{\{\chi\}}(t)\}$, and the steady state heat current into the $v$th bath is

[TABLE]

Comparing the dynamical equation of qubits density matrix at Eq. (C) with two-terminal case at Eq. (II.3), it is interesting to find that two equations are almost identical except expressions of the transition rates related with left bath(s) $\Gamma^{L}_{a,\pm}(\chi,\omega)$ at Eq. (18) and $\Gamma^{L}_{a,\pm}(\{\chi\},\omega)$ at Eq. (48). Hence, considering strong qubit-bath interaction, the heat flux into the $R$th bath can be expressed as

[TABLE]

where $\kappa^{ij}_{R}$ is identical with the counterpart at Eq. (35).

Then, we briefly discuss the nonmonotonic behavior of $I_{R}$ with strong qubit-bath coupling in Fig. 6(b). From Eq. (C), it is known that $I_{R}$ is composed by two components $E_{12}(\kappa^{12}_{R}P_{1}-\kappa^{21}_{R}P_{2})$ and $E_{34}(\kappa^{43}_{R}P_{4}-\kappa^{34}_{R}P_{3})$. Both show nonmonotonic features by increasing temperature $T_{R}$, shown in Fig. 11. In particular, the gap between two current components is first enhanced in the regime $T_{R}{\in}(0.2,0.5)$, and then decreases in the regime $T_{R}{\in}(0.5,1.5)$. This directly results in the turnover behavior of $I_{R}$ by tuning the temperature $T_{R}$.