Unifying quantum heat transfer and superradiant signature in a nonequilibrium collective-qubit system: a polaron-transformed Redfield approach

TL;DR

This paper develops a unified theoretical framework to analyze quantum heat transfer and superradiant effects in a collective-qubit system, revealing how system size influences superradiance and bridging weak and strong coupling regimes.

Contribution

It introduces a polaron-transformed Redfield approach with an auxiliary counting field to study heat transfer and fluctuations across different coupling strengths.

Findings

Heat current and fluctuations unify weak and strong coupling results.

Superradiant heat transfer depends on system size and diminishes with more qubits.

Large qubit numbers suppress superradiant signatures in nonequilibrium conditions.

Abstract

We investigate full counting statistics of quantum heat transfer in a collective-qubit system, constructed by multi-qubits interacting with two thermal baths. The nonequilibrium polaron-transformed Redfield approach embedded with an auxiliary counting field is applied to obtain the steady state heat current and fluctuations, which enables us to study the impact of the qubit-bath interaction in a wide regime. The heat current, current noise and skewness are all found to clearly unify the limiting results in the weak and strong couplings, respectively. Moreover, the superradiant heat transfer is clarified as a system-size-dependent effect, and large number of qubits dramatically suppresses the nonequilibrium superradiant signature.

Equations

$${\hat{H}}_s=\epsilon {\hat{J}}_z+\Delta {\hat{J}}_x,$$

$${\hat{V}}_u=2{\hat{J}}_z\sum _k(g_{k,u}{\hat{b}}_{k,u}^{†}+g_{k,u}^{\ast }{\hat{b}}_{k,u}),$$

$${\hat{H}}_s^{\prime }=\epsilon {\hat{J}}_z+\eta \Delta {\hat{J}}_x-\xi {\hat{J}}_z^2,$$

$$\eta =\exp [-2\sum _{k,u=L,R}|g_{k,u}/{\omega }_k{|}^2(2n_{k,u}+1)],$$

$${\hat{V}}_u^{\prime }=\Delta [(\cos \hat{B}-\eta )\widehat{J_x}+\sin \hat{B}{\hat{J}}_y].$$

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

$$C_x(\tau )$$

$$C_y(\tau )$$

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

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

$$\frac{d{P}_n}{dt}=-({\kappa }_{n-1}^{-}+{\kappa }_n^{+})P_n+{\kappa }_n^{-}{P}_{n+1}+{\kappa }_{n-1}^{+}{P}_{n-1},$$

$$Z(\{{\chi }_u\},t)=\text{Tr}[e^{i{\sum }_u{\chi }_u{\hat{H}}_u(0)}{e}^{-i{\sum }_u{\chi }_u{\hat{H}}_u(t)}{\hat{\rho }}_{tot}(0)],$$

$$Z(\{{\chi }_u\},t)=\text{Tr}[{\hat{U}}_{-\{{\chi }_u\}}{\hat{\rho }}_{tot}(0){\hat{U}}_{\{{\chi }_u\}}^{†}]=\text{Tr}[{\hat{\rho }}_{tot}^{\{{\chi }_u\}}(t)].$$

$$\Theta (\{{\chi }_u\})=\underset{t\to \infty }{\lim }\frac{1}{t}\ln Z(\{{\chi }_u\},t).$$

$$\displaystyle~{}J^{(n)}_{u}={{\frac{{\partial}^{n}\Theta(\{\chi_{u}\})}{{\partial}(i\chi_{u})^{n}}}}\bigg{|}_{\{\chi_{u}=0\}}.$$

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

$${\hat{V}}_u(\chi )=2{\hat{J}}_z\sum _k(g_{k,u}{e}^{i\frac{{\omega }_k\chi }{2}{\delta }_{u,R}}{\hat{b}}_{k,u}^{†}+H.c.),$$

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

$${\hat{V}}_u^{\prime }(\chi )=\Delta [(\cos {\hat{B}}_{\chi }-\eta ){\hat{J}}_x+\sin {\hat{B}}_{\chi }{\hat{J}}_y].$$

$$\frac{d{\hat{\rho }}_{\chi }}{dt}$$

$$C_x(\chi ,\pm \tau )$$

$$C_y(\chi ,\pm \tau )$$

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

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

$$\frac{d}{dt}|{\rho }_{\chi }(t)⟩=\hat{\mathcal{L}}(\chi )|{\rho }_{\chi }(t)⟩,$$

$$\Theta (\chi ,t)=\frac{1}{t}\ln [⟨I|e^{\hat{\mathcal{L}}(\chi )t}|\rho (0)⟩],$$

$$J_{\text{weak}}=4[n_L(\Delta )-n_R(\Delta )]\frac{G_L(\Delta )G_R(\Delta )}{G_L(\Delta )+G_R(\Delta )}\frac{(2-x)(1+x^3)}{(1+x+x^2)},$$

$$\frac{d{P}_n^{\chi }(t)}{dt}$$

$${\kappa }_n^{\pm }(\chi )$$

$${\Theta }_{\text{NIBA}}(\chi )$$

$$J_{\text{strong}}=\frac{{\Delta }^2}{2\pi }{\int }_{-\infty }^{\infty }d\omega \omega \left[\frac{{\kappa }_0^{-}}{2{\kappa }_0^{+}+{\kappa }_0^{-}}{C}_R(\omega )C_L(-\omega +\xi )-\frac{{\kappa }_0^{+}}{2{\kappa }_0^{+}+{\kappa }_0^{-}}{C}_R(-\omega )C_L(\omega -\xi )\right],$$

Full text

I Introduction

Understanding mechanism of nonequilibrium quantum transport in low dimensional systems is a long-standing challenge, which has been extensively investigated from solid-state physics apjauho2008 , quantum thermodynamics rkosloff2013entropy , molecular electronics anitzan2014 to quantum biology mmohseni2014 . In particular, quantum thermal transport where the particle and heat flows are modulated by the temperature bias, has triggered the emergence of the phononics jswang2008epjb ; nbli2012rmp ; jren2015aip . Phonons have been successfully utilized to fabricate various functional devices, e.g., thermal diode, memory and transistor bli2004prl ; bli2006apl ; lwang2007prl ; lwang2008prl . In analogy with phononics, quantum heat transfer mediated by the spin(e.g., anharmonic molecule or qubit) has also been intensively analyzed in a parallel way, which leads to the realization of quantum spin-thermal transistor kjoulain2016prl ; cwang2018pra , spin heat engine dsegal2008prl ; bsothmann2012epl ; jren2013prb and nonequilibrium spin network gtcraven2016pnas ; gtcraven2017prl .

As a generic model to describe quantum heat transfer at the nanoscale, the nonequilibrium spin-boson model(NESB) is composed by a two-level-system(TLS) coupled to two thermal baths dsegal2005prl , which was originally proposed to study the quantum dissipation ajleggett1987rmp ; pao1989prl ; uweiss2008 . Many approaches have been explored to investigate the underlying mechanism of heat transfer in the NESB dsegal2006prb ; mgalperin2007prb ; kvelizhanin2008cpl ; ksaito2013prl ; yyao2015prb ; etaylor2015prl ; akato2015jcp ; akato2016jcp ; dzxu2016fp ; cwang2017pra ; jliu2017pre ; lferialdi2017pra . Analytically, the Redfield scheme is properly introduced to investigate the sequential process in the weak spin-bath interaction regime, where two thermal baths show additive contribution to the heat transfer dsegal2008prl ; jren2010prl . However, the Redfield approach breaks down in strong spin-bath interaction regime, where the multiphonon excitations should be involved to characterize the nonequilibrium heat-exchange. Then, the nonequilibrium noninteracting blip approximation(NIBA) can be applied to study baths induced nonadditive and cooperative transfer processes dsegal2006prb ; lnicolin2011jcp ; lnicolin2011prb ; tchen2013prb ; dsegal2014pre . However, both methods are found to have their own limitations, i.e., the Redfield approach is unable to capture the multi-phonon processes in strong coupling regime, and the heat current based on the nonequilibrium NIBA scheme does not show linear proportion to coupling strength in weak coupling regime tchen2013prb . Recently, the nonequilibrium polaron-transformed Redfield equation(NE-PTRE) was proposed to successfully unify the steady state heat current in the NESB cwang2017pra ; cwang2015sr ; dzxu2016njp .

However, the exploration of the NE-PTRE to study more nonequilibrium spin systems is lacking, but is urgently required for the spin-based quantum heat transfer. This paper aims to fill this gap by applying the NE-PTRE to analyze the heat transfer in the collective-qubit model.

Recently, a superradiant signature of quantum heat transfer has been discovered in the collective-qubit system with weak qubit-bath interaction mvogl2011aop ; mvogl2012pra . The steady state superradiance describes the effect that qubits collectively exchange energy with thermal baths, resulting in the current scaling as $J{\sim}N^{2}_{s}$ with $N_{s}$ the number of qubits. In sharp contrast, the superradiant heat transfer vanishes in strong coupling regime, based on the nonequilibrium NIBA under the Marcus approximation cwang2015aop . Hence, the steady state behavior of the collective-qubit system are significant distinct from each other in limiting interaction regimes. It is demanding to analyze the heat transfer feature from weak to strong couplings from a unified perspective

To address these problems, we extend the NE-PTRE combined with full counting statistics(FCS) mesposito2009rmp ; mcampisi2011rmp to investigate quantum heat transfer in the collective-qubit system. The counting-field dependent NE-PTRE successfully unifies the current and fluctuations(e.g., heat current, current noise, skewness), with the incoherent picture in the weak quit-bath coupling limit and multi-phonon excited transfer picture in the strong coupling limit, respectively. Moreover, the superradiant heat transfer is investigated at large temperature bias mvogl2011aop , and the disappearance of the superradiant signature is explained by enlarging the number of qubits beyond the weak qubit-bath interaction.

This paper is organized as follows: In Sec. II, the collective-qubit model and the NE-PTRE are described. In Sec. III, FCS is briefly introduced and the counting-field dependent NE-PTRE is derived, which enables us to analyze the steady state heat transfer. In Sec. IV, heat current, current noise and skewness are all found to hold the unified features, by extending the application of the NE-PTRE. In Sec. V, we study the transition of superradiant heat transfer from weak to strong couplings, and explain the vanishing mechanism of superradiant signature with large number of qubits. In the final section, we present a concise summary.

II Nonequilibrium collective-qubit system

II.1 Model

The nonequilibrium energy transfer in the collective-qubit system, which interacts with two thermal baths, is modeled as $\hat{H}=\hat{H}_{s}+\sum_{u=L,R}(\hat{H}^{u}_{b}+\hat{V}_{u})$. The collective-qubit model is described as

[TABLE]

where the collective angular-momentum operators are $\hat{J}_{a}=\frac{1}{2}\sum^{N_{s}}_{i}\hat{\sigma}^{i}_{a}~{}(a=x,y,z)$ with $\hat{\sigma}^{i}_{a}$ the Pauli operator of the $i$th qubit, $N_{s}$ is the number of qubits, $\varepsilon$ and $\Delta$ are the Zeeman splitting energy and coherent tunneling strength of the angular momentum operator. In the limit of $N_{s}=1$, the Hamiltonian at Eq. (1) is reduced to the seminal nonequilibrium spin-boson model dsegal2005prl ; jren2010prl . The $u$th thermal bath is composed by noninteracting bosons, shown as $\hat{H}^{u}_{b}=\sum_{k}\omega_{k}\hat{b}^{{\dagger}}_{k,u}\hat{b}_{k,u}$, where $\hat{b}^{{\dagger}}_{k,u}(\hat{b}_{k,u})$ creates (annihilates) one phonon with the frequency $\omega_{k}$. The system-bath interaction is given by

[TABLE]

where $g_{k,u}$ is the coupling strength between the angular-momentum and the $u$th thermal bath. It is characterized by the spectral function $G_{u}(\omega)=4\pi\sum_{k}|g_{k,u}|^{2}\delta(\omega-\omega_{k})$. In this paper, we select the spectral function as the super-Ohmic form $G_{u}(\omega)=\pi\alpha_{u}\omega^{3}/\omega^{2}_{c}\exp(-\omega/\omega_{c})$, where $\alpha_{u}$ is the coupling strength and $\omega_{c}$ is the cutoff frequency of thermal baths. The super-Ohmic spectrum has been extensively considered to investigate the quantum dissipative dynamics and transport in molecular electronics anazir2009prl ; sjjang2013njp ; cklee2015jcp , solid-state devices dpsm2013prl and light-harvesting complexes sjjang2011jcp ; dpsm2011prb .

To analyze the multi-phonon involved energy transfer processes, we apply the canonical transformation $\hat{U}=\exp[i\hat{J}_{z}\sum_{u}\hat{B}_{u}]$ to the Hamiltonian $\hat{H}$ as $\hat{H}^{\prime}=\hat{U}^{{\dagger}}\hat{H}\hat{U}=\hat{H}^{\prime}_{s}+\sum_{u=L,R}(\hat{H}^{u}_{b}+\hat{V}^{\prime}_{u})$ rsilbey1984jcp ; rharris1985jcp , where the collective phonon momentum is $\hat{B}_{u}=2i\sum_{k}(\frac{g_{k,u}}{\omega_{k}}\hat{b}^{{\dagger}}_{k,u}-H.c.)$. After the transformation, the modified system Hamiltonian is given by

[TABLE]

where the factor $\eta={\langle}\cos{\hat{B}}{\rangle}$ is specified as

[TABLE]

and the renormalization energy is $\xi=4\sum_{k,u}|g_{k,u}|^{2}/\omega_{k}$. In weak system-bath interaction limit $|g_{k,u}/\omega_{k}|^{2}{\approx}0$, the renormalization factor becomes $\eta=1$ and $\xi{\approx}0$. While in strong system-bath interaction region $|g_{k,u}/\omega_{k}|^{2}{\gg}1$, the factor $\eta$ approaches zero, and the renormalization energy $\xi{\gg}\{\varepsilon,\Delta\}$. The generation of large $\xi$ will dramatically change the energy structure of the modified qubits system, compared to the counterpart in the weak coupling limit. Moreover, the modified system-bath interaction is given by

[TABLE]

II.2 Nonequilibrium polaron-transformed Redfield equation

We apply the NE-PTRE to investigate the dynamics of the collective-qubit model. The NE-PTRE, one type of the quantum master equation, has bee successfully applied to unify the nonequilibrium energy transfer in the seminal spin-boson model cwang2017pra ; dzxu2016njp . It is known that the modified system-bath interaction disappears under the thermal average (${\langle}\hat{V}^{\prime}_{u}{\rangle}$=0). Hence, it may be safe to perturb ${\langle}\hat{V}^{\prime}_{u}{\rangle}$ up to the second order to obtain the quantum master equation. Under the Born approximation, the density operator of the whole system can be decomposed by $\hat{\rho}(t)=\hat{\rho}_{s}(t){\otimes}\hat{\rho}_{b}$ with $\hat{\rho}_{s}(t)$ the density operator of the qubits and $\hat{\rho}_{b}=e^{-\sum_{u}\hat{H}^{u}_{b}/k_{B}T_{u}}/\textrm{Tr}_{b}\{e^{-\sum_{u}\hat{H}^{u}_{b}/k_{B}T_{u}}\}$ the density operator of the baths at equilibrium. The quantum master equation in the Markovian limit is obtained as

[TABLE]

where the correlation functions are

[TABLE]

with the phonon propagator in $u$th thermal bath $Q_{u}(\tau)=4\sum_{k}|g_{k,u}/\omega_{k}|^{2}[\cos\omega_{k}\tau(2n_{k,u}+1)-i\sin\omega_{k}\tau]$. Furthermore, in the eigen-basis $\hat{H}^{\prime}_{s}|\phi_{n}{\rangle}=E_{n}|\phi_{n}{\rangle}$, the dynamical equation can be re-expressed as

[TABLE]

where the transition rate is $\Gamma_{a}(\omega)=\int^{\infty}_{0}{d\tau}C_{a}(\tau)e^{-i\omega\tau}$ and the element is ${\rho}_{nn^{\prime}}={\langle}\phi_{n}|\hat{\rho}_{s}(t)|\phi_{n^{\prime}}{\rangle}$. The rate $\Gamma_{x(y)}(E_{nm})$ describes that even(odd) number phonons are involved in the transfer process between the states $|\phi_{n}{\rangle}$ and $|\phi_{m}{\rangle}$.

In the weak qubit-bath coupling limit, the heat transfer is dominated by the sequential process and $\eta{\approx}1$. Thus, the correlation function $C_{y}(\tau)$ is reduced to $C_{y}(\tau){\approx}\Delta^{2}[\sum_{u}Q_{u}(\tau)]$, and $C_{x}(\tau){\approx}0$ by ignoring high-order correlations. Accordingly, the lowest order of the transition rate $\Gamma^{(1)}_{y}(\omega)$ includes the term $\Delta^{2}(Q_{L}(\omega)+Q_{R}(\omega))$, with $Q_{v}(\omega)=\int{d\tau}e^{-i\omega\tau}Q_{v}(\tau)$. Moreover, all of off-diagonal elements of the qubits density matrix approach zero at steady state (not shown here, which is quite similar to the result at Fig. 3). Considering the commutating relation $[\hat{H}^{\prime}_{s},\hat{J}_{z}]=-i{\Delta}\hat{J}_{y}$, the quantum master equation at Eq. (8) after long time evolution is simplified as

[TABLE]

which is identical with the dynamical equation based on the Redfield scheme (see appendix A for the detail).

While in strong qubit-bath interaction limit, the coherent tunneling of qubits at Eq. (3) is dramatically suppressed ($\eta{\approx}0$), and the correlation factor $\eta^{2}e^{-\sum_{u}Q_{u}(\tau)}$ vanishes. However, the other factor $\eta^{2}e^{\sum_{u}Q_{u}(\tau)}$ is kept finite, which contributes to the quantum heat transfer. Hence, the correlation functions at Eq. (7) are reduced to $C_{x}(\tau)=C_{y}(\tau)=(\eta\Delta)^{2}[\exp{\sum_{u}Q_{u}(\tau)}]$. Consequently, the master equation at Eq. (8) in the local basis $\{|\phi_{n}{\rangle}\}$ with $(\varepsilon\hat{J}_{z}-\xi\hat{J}^{2}_{z})|\phi_{n}{\rangle}=(\varepsilon{n}-\xi{n}^{2})|\phi_{n}{\rangle}$ is changed to

[TABLE]

where the rate is $\kappa^{\pm}_{n}=\frac{(j^{+}_{n}\Delta)^{2}}{8\pi}\int{d\omega}C_{L}(\pm\omega\mp\Delta_{n})C_{R}(\mp\omega)$, with $j^{+}_{n}=\sqrt{\frac{N_{s}}{2}(\frac{N_{s}}{2}+1)-n(n+1)}$, the energy gap $\Delta_{n}=\varepsilon-\xi(2n+1)$ and the correlation function in the frequency domain $C_{v}(\omega)=\int{dt}e^{i\omega{t}}\eta^{2}_{v}e^{Q_{v}(t)}$. Hence, the dynamical equation is completely reduced to the nonequilibrium NIBA scheme (see Eq.(B) at appendix B for the detail).

III Full counting statistics of heat current

III.1 The general theory

Full counting statistics is a two-time projection protocol to measure the current and fluctuations mesposito2009rmp ; mcampisi2011rmp . For the energy transfer in the multi-terminal setup, the generating function is generally given by hmf2018njp

[TABLE]

where $\chi_{u}$ is the counting-field parameter to count the energy flow into the $u$th bath with the Hamiltonian $\hat{H}_{u}$, $\hat{H}_{u}(t)=\hat{U}^{{\dagger}}\hat{H}_{u}\hat{U}$ with the propagator $\hat{U}=e^{-i\hat{H}t}$, and $\hat{\rho}_{tot}(0)$ is the initial density matrix of the whole system. Moreover, considering the modified propagator $\hat{U}_{\{\chi_{u}\}}(t)=e^{i\sum_{u}\chi_{u}\hat{H}_{u}/2}\hat{U}(t)e^{-i\sum_{u}\chi_{u}\hat{H}_{u}/2}$, it can be re-expressed as $\hat{U}_{\{\chi_{u}\}}(t)=\exp{(-i\hat{H}_{\{\chi_{u}\}}t)}$ with $\hat{H}_{\{\chi_{u}\}}=e^{i\sum_{u}\chi_{u}\hat{H}_{u}/2}\hat{H}e^{-i\sum_{u}\chi_{u}\hat{H}_{u}/2}$. Hence, the generating function is re-expressed as

[TABLE]

After the long-time evolution, the cumulant generating function is obtained as

[TABLE]

Therefore, the $n$th cumulant of heat current fluctuations into the $u$th bath is given by

[TABLE]

Specifically, the lowest cumulant is the steady state heat current J_{u}={\frac{{\partial}\Theta(\{\chi_{u}\})}{{\partial}(i\chi_{u})}}\Big{|}_{\{\chi_{u}=0\}}, the second cumulant is the zero-frequency current noise J^{(2)}_{u}={\frac{{\partial}^{2}\Theta(\{\chi_{u}\})}{{\partial}(i\chi_{u})^{2}}}\Big{|}_{\{\chi_{u}=0\}} and the third cumulant is the current skewness J^{(3)}_{u}={\frac{{\partial}^{3}\Theta(\{\chi_{u}\})}{{\partial}(i\chi_{u})^{3}}}\Big{|}_{\{\chi_{u}=0\}}. In the following, we apply the theory of the FCS to obtain the counting-field dependent NE-PTRE of the nonequilibrium collective-qubit model.

III.2 Counting-field dependent NE-PTRE

To count the heat flow into the right thermal bath, we introduce the counting field parameters as $\hat{H}_{\chi}=e^{i\hat{H}^{R}_{b}\chi/2}\hat{H}e^{-i\hat{H}^{R}_{b}\chi/2}$ hmf2018njp , which results in

[TABLE]

where the system-bath interaction with the counting field parameter is expressed as

[TABLE]

where $\delta_{R,R}=1$ and $\delta_{L,R}=0$. Then, after the canonical transformation $\hat{H}^{\prime}_{\chi}=U^{{\dagger}}_{\chi}\hat{H}_{\chi}{U_{\chi}}$ with transformation operator $U=\exp[i\hat{J}_{z}\sum_{u}\hat{B}_{u}(\chi)]$ and $\hat{B}_{u}(\chi)=2i\sum_{k}(\frac{g_{k,u}}{\omega_{k}}e^{i\frac{\omega_{k}\chi}{2}\delta_{u,R}}\hat{b}^{{\dagger}}_{k,u}-H.c.)$, the modified Hamiltonian is given by

[TABLE]

where the modified system-bath interaction with counting field parameter is given by

[TABLE]

By perturbing $\hat{V}^{\prime}_{u}(\chi)$ under the Born-Markov approximation, we obtain the second-order quantum master equation as

[TABLE]

where the correlation functions with counting field parameter are

[TABLE]

In absence of the counting field parameter, the correlation functions are reduced to the standard version at Eq. (7). In the eigen-basis, the dynamics of density matrix elements can be specified as

[TABLE]

where the transition rates are

[TABLE]

If we reorganize the reduced density matrix of collective-qubit from the matrix form to the vector expression, the dynamical equation at Eq. (21) is expressed as

[TABLE]

where $\hat{\mathcal{L}}(\chi)$ is the superoperator to dominate the system evolution. Then, the cumulant generating function at $t$-time is given by

[TABLE]

where $|\rho(0){\rangle}$ is the vector form of the initial system density matrix, and ${\langle}I|$ is the left-eigenvector of $\hat{\mathcal{L}}$ as ${\langle}I|\hat{\mathcal{L}}=0$, with the normalization relation ${\langle}I|\rho_{\chi=0}(t){\rangle}=1$. Hence, heat current fluctuations at steady state can be straightforwardly obtained by following the Eq. (14).

IV Unified steady state heat current

Quantum heat transfer in the NESB has been successfully investigated based on the Redfield and noninteracting-blip approximation schemes in the weak and strong qubit-bath coupling limits, separately. However, steady state heat current was found to be distinct from each other in a broad coupling regime tchen2013prb . Recently, the nonequilibrium polaron-transformed Redfield equation combined with full counting statistics was proposed to unify heat current between these two limiting approaches cwang2017pra ; dzxu2016njp .

Here, we try to extend the counting-field dependent NE-PTRE to unify the heat transfer in the nonequilibrium collective-qubit model in Fig. 1. At resonance($\varepsilon=0$), we first analyze the steady state heat current in Figs. 1(a) and 1(d). It clearly exhibits the turnover behavior, which unifies the counterparts within the Redfield and NIBA methods as the qubit-bath coupling strength approaches the weak and strong coupling limits. We admit that to gain an analytical expression of heat current at arbitrary qubit-bath coupling is rather difficult. However, it can be obtained in limiting regimes.

Here, we study the analytical expression of the steady state heat current with $N_{s}=2$. Specifically, in the weak interaction limit the heat current is analytically expressed as(see Eq. (39) at appendix A)

[TABLE]

where the coefficient $x=[\sum_{v=L,R}G_{v}(\Delta)n_{v}(\Delta)]/[\sum_{v=L,R}G_{v}(\Delta)(1+n_{v}(\Delta))]$, which is linearly proportional to the qubit-bath coupling strength. The current at Eq. (25) is the special case of the general expression at Eq. (39). It is found in Fig.1(a) that $J_{\textrm{weak}}$ becomes identical with the counterpart from NE-PTRE in weak coupling limit (e.g., $\alpha=0.001$). Moreover, it should be noted that the current $J_{\textrm{weak}}$ at Eq. (25) with $N_{s}=2$ has the similar expression to the case in the standard NESB($N_{s}=1$) in the weak coupling limit, which both proportional to the thermodynamic bias (i.e., $J{\propto}[n_{L}(\Delta)-n_{R}(\Delta)]$) dsegal2006prb .

While in the strong coupling limit, the dynamical equation at Eq. (19) with the number of qubits $N_{s}=2$ is reduced to the kinetic form as (see Eq. (B) at appendix B)

[TABLE]

where the population is $P^{\chi}_{n}={\langle}1,n|\hat{\rho}_{\chi}(t)|1,n{\rangle}$ with $\hat{J}_{z}|1,n{\rangle}=n|1,n{\rangle}$. The transition rate is

[TABLE]

where the coefficient is $j^{+}_{n}=\sqrt{2-n(n{+}1)}$, the energy gap $\Delta_{n}=E_{n+1}-E_{n}={\varepsilon}-\xi(2n+1)$ and the correlation function in the frequency domain is $C_{v}(\omega)=\int^{\infty}_{-\infty}dte^{i\omega{t}}\eta^{2}_{v}e^{Q_{v}(t)}$ with the factor $\eta_{v}=\exp{[-2\sum_{k}|g_{k,v}/\omega_{k}|^{2}(2n_{k,v}+1)]}$. If we select $N_{s}=1$($j^{+}_{-1/2}=1$), the rate $\kappa^{\pm}_{n}(\chi)$ is reduced to the standard NESB case(see Eq. (20) in Ref. lnicolin2011jcp ). Thus, the cumulant generating function at steady state is given by

[TABLE]

where $\kappa^{\pm}_{-1}(\chi)=\kappa^{\mp}_{0}(\chi)$. Consequently, the heat current is given by

[TABLE]

where the first(second) term describes the process that as the qubits release(gain) energy $\xi$, the right bath absorbs(emits) phonon energy $\omega$ and the left bath obtains(provides) the remaining energy $\xi-\omega$. $J_{\textrm{strong}}$ shows similar structure with the counterpart in the standard NESB, which is jointly contributed by two thermal baths lnicolin2011jcp . Moreover, the expression of $J_{\textrm{strong}}$ captures the turnover feature of heat current in Fig. 1(a) and 1(d) with dashed blue lines with squares. Hence, we conclude that the steady state heat flux is clearly unified in the nonequilibrium collective-qubit model.

Next, we analyze the zero-frequency current noise and skewness in Figs. 1(b-c) and 1(e-f). It is interesting to find that the results based on the NE-PTRE also perfectly bridge the limiting counterpart at the weak and strong coupling regimes, which may demonstrate the unification of current fluctuations in extended spin-boson systems (e.g., collective-qubit model). Moreover, the turnover behavior which is presented in the current, are unraveled for the second and third cumulants. Though not shown here, higher cumulants of current fluctuations also show such unified features. We should note that all above results are valid both at resonant and off-resonant conditions, which clearly exhibits that full counting statistics of the heat current is generally unified within the NE-PTRE scheme.

V Suppression of superradiant heat transfer

Superradiant effect, which describes that the system exhibits collective response under the modulation of the external field or thermal bath, has been extensively investigated in quantum phase transition nlambert2004prl ; qhchen2008pra , critical heat engine auch2015sr and quantum transport mvogl2011aop ; mvogl2012pra . In particular, quantum heat transfer in the nonequilibrium large-spin system shows the superradiant signature in the weak spin-bath coupling regime mvogl2011aop . Specifically, under large temperature bias (e.g., $T_{L}{\gg}T_{R}$) the steady state heat current is expressed as

[TABLE]

which follows the condition $x/(1-x){\gg}N_{s}$. It should be noted that Eq. (30) is a special case of the current at Eq. (39) at appendix A. However, with the strong spin-bath interaction based on the NIBA scheme combined with the Marcus approximation, it is found that such superradiant transfer vanishes cwang2015aop . Hence, we apply the NE-PTRE to clarify such seemingly paradox.

We first study the effect of the temperature bias $\Delta{T}$ on the heat current in Fig. 2(a) with spin-bath coupling strength $\alpha=0.01{\ll}\{\Delta,\omega_{c}\}$, which is considered weak for the seminal spin-boson model($\eta{\approx}1$). It is found that the current shows monotonic enhancement by increasing both ${\Delta}T$ and $N_{s}$. Moreover, at large temperature bias (e.g., $T_{L}=8$ and $T_{R}=0.4$) the current becomes nearly stable with large qubits number, which is clearly exhibited in Fig. 2(b) (e.g., $N_{s}=32$). While for relatively small number of qubits (e.g., $N_{s}{<}20$) in Fig. 2(b), the superradiant signature of heat current is numerically obtained as $J{\propto}N^{\gamma}_{s}$ with $\gamma=2.0{\pm}0.1$ (both for resonance and off-resonance). Thus, it is known that the expression of superradiant heat current at Eq. (30) becomes invalid at large $N_{s}$.

To exploit the origin, we plot the current as a function of the qubit-bath coupling strength $\alpha$ for different $N_{s}$ with $T_{L}=8$, $T_{R}=0.4$ and $\varepsilon=0$, as shown in Fig. 2(c). It is shown that for small $N_{s}$(e.g., $N_{s}=8$), the heat flux exhibits approximately linear increasing behavior in the weak coupling regime (e.g., $\alpha{\leq}0.01$). The heat flux demonstrates the sequential transfer process, where a superradiant heat transfer is observed accordingly, which could be described by the Redfield scheme. While the current shows distinct behavior for large $N_{s}$ from that for small $N_{s}$ in the same coupling regime. For example, the current for $N_{s}=32$ with the coupling strength $\alpha=0.01$ has already surpassed the turnover point of the current, while the current for $N_{s}=8$ is almost linearly dependent on the coupling strength in same coupling regime. It is known that the appearance of the turnover point of the current is the significant signature of the multi-phonons involved coherent transfer, as exploited in the nonequilibrium spin-boson model cwang2017pra . Since $\eta$ is nearly equal to $1$ in above two cases and the heat current shows different features in same coupling regime, we conclude that $\eta{\approx}1$ does not necessarily correspond to the weak coupling condition. Furthermore, the behavior of the current should be properly described by the NIBA scheme. It results in the absence of nonequilibrium superradiant signature. Therefore, the effect of superradiant heat transfer in the collective-qubit model will be dramatically suppressed in the large $N_{s}$ regime, i.e. it is a $N_{s}$-dependent phenomenon.

VI Conclusion

To summarize, we investigate the quantum heat transfer in the nonequilibrium collective-qubit system, by applying the nonequilibrium polaron-transformed Redfield equation combined with full counting statistics. We first analyze the effect of qubit-bath coupling on the steady state heat current, which resulting in a turnover behavior, which is similar to the counterpart in the nonequilibrium spin-boson model. Interestingly, the current consistently bridges the results in the weak and strong coupling limits, which clearly demonstrates that heat current can be unified in the multi-qubits case. It should be admitted that the general solution of the heat current is difficult to obtain even for $N_{s}=2$. However, the analytical expression can be still obtained in the weak and strong couplings based on the Redfield(Eq. (25)) and NIBA(Eq. (29)) schemes, respectively. Moreover, the current noise and skewness are shown to be unified accordingly. Though not shown in the present paper, the unification of higher cumulants of heat current can also be observed. Therefore, we propose that full counting statistics of heat current at steady state can be unified by using the NE-PTRE.

Next, we study the superradiant heat transfer at high temperature bias regime. It is found that with small number of qubits, the heat transfer is described by the sequential process under the Redfield scheme. The heat current exhibits $J{\sim}N^{2.0{\pm}0.1}_{s}$, an apparent signature of the steady state superradiance. While with the large number of qubits, the superradiant signature of the heat flux vanishes. The corresponding physical process is described by the NIBA scheme, for multi-phonons should be involved to contribute to the heat transfer. Therefore, we conclude that the superradiant transfer in the collective-qubit model is a size-dependent phenomenon, and it will be strongly suppressed by increasing the qubits number.

We believe that based on the counting-field dependent NE-PTRE, the unified feature of steady state heat transfer may be realized in a much bigger family of the quantum nonequilibrium system, e.g., quantum spin-boson network gtcraven2017prl .

VII Acknowledgement

X.M.C acknowledges the support by the National Natural Science Foundation of China(No. 11874011), and C.W. is supported by the National Natural Science Foundation of China under Grant No. 11704093.

Appendix A Quantum thermal transfer within the Redfield scheme

The nonequilibrium collective-qubit system is expressed as $\hat{H}=\hat{H}_{s}+\sum_{v=L,R}(\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 the 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 directly perturb the system-bath interaction at Eq. (33) up to the second order. Then, based on the Born-Markov approximation, the Redfield equation is given by

[TABLE]

with the commutating relation $[\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]

where the counting-field dependent density matrix element is $\rho^{\chi}_{nn^{\prime}}={\langle}\phi_{n}|\hat{\rho}_{\chi}(t)|\phi_{n^{\prime}}{\rangle}$, the transition rate is $G_{v}(\omega)=\pi\alpha_{v}\omega^{3}/\omega^{2}_{c}\exp{(-|\omega|/\omega_{c})}$, and the Bose-Einstein distribution function is $n_{v}(\omega)=1/[\exp(\omega/k_{B}T_{v})-1]$.

At resonance ($\varepsilon=0$), it is found that off-diagonal elements of the qubits system in eigen-space vanish at steady state in Figs. 3(a), (b) and (c), in absence of the counting field ($\chi=0$). Hence, the steady state equation is given by

[TABLE]

Actually, the system Hamiltonian at resonance is $\hat{H}_{s}=\Delta\hat{J}_{x}$, with the eigen-solution $\hat{H_{s}}|j,m{\rangle}_{x}=\Delta{m}|j,m{\rangle}_{x}~{}(m=-N_{s}/2,-N_{s}/2+1,\cdots,N_{s}/2)$. Thus, the coefficient $J^{nm}_{z}$ can be specified as $J^{nm}_{z}=\frac{1}{2}(j^{+}_{m}\delta_{n,m+1}+j^{-}_{m}\delta_{n,m-1})$, with $j^{\pm}_{m}=\sqrt{N_{s}/2(N_{s}/2+1)-m(m{\pm}1)}$. Consequently, the steady state population can be analytically obtained as

[TABLE]

with the coefficient

[TABLE]

Finally, the steady state heat flux is given by

[TABLE]

where the current factor is expressed as

[TABLE]

While at the off-resonant condition (e.g., $\varepsilon=1$), the off-diagonal elements after long time evolution also vanish, shown at Figs. 3(d), (e) and (f). Hence, the heat current into the right thermal bath is generally expressed as

[TABLE]

Appendix B Quantum thermal transfer within the NIBA scheme

Starting from the counting-field dependent Hamiltonian $\hat{H}_{\chi}=\hat{H}_{s}+\sum_{u=L,R}(\hat{H}^{u}_{b}+\hat{V}_{u}(\chi))$ at Eq. (15), we apply a canonical transformation $\hat{U}=\exp[i\hat{J}_{z}\sum_{u=L,R}\hat{B}_{u}(\chi)]$ with $\hat{B}_{u}(\chi)=2i\sum_{k}(\frac{g_{k,u}}{\omega_{k}}e^{i\frac{\omega_{k}\chi}{2}\delta_{u,R}}\hat{b}^{{\dagger}}_{k,u}-H.c.)$, to obtain the modified Hamiltonian of the whole system $\hat{H}^{\prime}_{\chi}=\hat{U}^{{\dagger}}_{\chi}\hat{H}_{\chi}\hat{U}_{\chi}$ as

[TABLE]

The modified system Hamiltonian is given by

[TABLE]

where the renormalization energy is $\xi=4\sum_{k,u=L,R}|g_{k,u}|^{2}/\omega_{k}$. The eigen-solution is given by $\hat{H}^{\textrm{NIBA}}_{s}|\phi_{n}{\rangle}=E_{n}|\phi_{n}{\rangle}$ with $E_{n}=\varepsilon{n}-\xi{n}^{2}$ and $n=-N_{s}/2,{\cdots},N_{s}/2$. The modified system-bath interaction is given by

[TABLE]

Hence, perturbing $\hat{V}^{\chi}_{sb}$ up to the second order, we obtain the quantum kinetic equation

[TABLE]

where the counting-field dependent population element is $P^{\chi}_{n}={\langle}\phi_{n}|\hat{\rho}_{\chi}(t)|\phi_{n}{\rangle}$ and the modified transition rates are

[TABLE]

where the coefficient is $j^{+}_{n}={{N_{s}}/{2}({N_{s}}/{2}+1)-n(n{+}1)}$, the energy gap $\Delta_{n}=E_{n+1}-E_{n}={\varepsilon}-\xi(2n+1)$, the correlation function in the frequency domain is

[TABLE]

with the renormalization factor $\eta_{v}=\exp[-2\sum_{k}|\frac{g_{k,u}}{\omega_{k}}|^{2}(2n_{k,u}+1)]$, and the correlation phase

[TABLE]

In absence of the counting-field parameter ($\chi=0$), this modified kinetic equation is identical with the dynamical equation at Eq. (10). Hence, the heat current at steady state is obtained as

[TABLE]

where the steady state population is given by $P^{ss}_{n}={\langle}N_{s}/2,n|\hat{\rho}_{s}(t{\rightarrow}\infty)|N_{s}/2,n{\rangle}$.