A Nonequilibrium Variational Polaron Theory to Study Quantum Heat Transport

TL;DR

This paper introduces a nonequilibrium variational polaron approach to accurately analyze quantum heat transport at the nanoscale, overcoming limitations of previous models and capturing complex effects like cross-bath correlations.

Contribution

It develops a variational polaron framework for nonequilibrium steady states that handles infrared divergence and cross-bath effects, extending beyond super-Ohmic bath models.

Findings

More accurate heat current calculations than NIBA for Ohmic baths

Demonstration of cross-bath correlation effects

Observation of current turnover and rectification phenomena

Abstract

We propose a nonequilibrium variational polaron transformation, based on an ansatz for nonequilibrium steady state (NESS) with an effective temperature, to study quantum heat transport at the nanoscale. By combining the variational polaron transformed master equation with the full counting statistics, we have extended the applicability of the polaron-based framework to study nonequilibrium process beyond the super-Ohmic bath models. Previously, the polaron-based framework for quantum heat transport reduces exactly to the non-interacting blip approximation (NIBA) formalism for Ohmic bath models due to the issue of the infrared divergence associated with the full polaron transformation. The nonequilibrium variational method allows us to appropriately treat the infrared divergence in the low-frequency bath modes and explicitly include cross-bath correlation effects. These improvements provide more accurate calculation of heat current than the NIBA formalism for Ohmic bath models. We illustrate the aforementioned improvements with the nonequilibrium spin-boson model in this work and quantitatively demonstrate the cross-bath correlation, current turnover, and rectification effects in quantum heat transfer.

Equations

$$H$$

$${\rho }_{NESS}$$

$$Y=\frac{1}{2}\Omega \left(H_B^L-H_B^R\right){\Omega }^{-1},$$

$$H$$

$${\gamma }_v(\omega )=2\pi \sum _k{g}_{k,v}^2\delta (\omega -{\omega }_{k,u}).$$

$${\gamma }_v(\omega )=\pi {\alpha }_v{\omega }^s{\omega }_{c,v}^{1-s}{\psi }_c\left(\omega /{\omega }_{c,v}\right)$$

$${\rho }_{NESS}$$

$$T_{eff}=\frac{({\alpha }_L{T}_L+{\alpha }_R{T}_R)}{({\alpha }_L+{\alpha }_R)}$$

$$C_v(t)=\frac{2}{\pi }\int {\gamma }_v(\omega )[\coth (\omega T_v)\cos (t\omega )-i\sin (t\omega )]d\omega$$

$$C_R(t)+C_L(t)\approx \frac{2}{\pi }\int [{\gamma }_R(\omega )+{\gamma }_L(\omega )][\coth (\omega T_{eff})\cos (\omega t)-i\sin (\omega t)]d\omega$$

$${\dot{\rho }}_{12}(t)$$

$${\dot{\rho }}_1(t)$$

$${\dot{\rho }}_2(t)$$

$${\rho }_{12}(t)$$

$${\dot{P}}_1(t)$$

$$k_{\pm }$$

$$g_v(t)$$

$$k_{\pm }=2\pi V^2\sqrt{\frac{\bar{\beta }}{4\pi \alpha }}\exp \left(-\bar{\beta }\frac{(\alpha \pm ϵ{)}^2}{4\alpha }\right),$$

$$U$$

$${\tilde{H}}_s$$

$${\tilde{H}}_I$$

$$V_x$$

$$V_y$$

$$V_z$$

$$A\le A_B=-\frac{1}{\beta }\ln \text{Tr}\left(e^{-\beta {\tilde{H}}_s}\right)+⟨{\tilde{H}}_I{⟩}_0+\mathcal{O}(⟨{\tilde{H}}_I^2{⟩}_0),$$

$$F_v(\omega )$$

$$\eta =\exp \left(-\sum _v\int d\omega \frac{{\gamma }_v(\omega )}{\pi {\omega }^2}{F}_v^2(\omega )\coth \left(\frac{{\beta }_v}{2}\omega \right)\right).$$

$$\frac{d{\rho }_s}{dt}=-i\left[H_s,{\rho }_s\right]+$$

$$\sum _{\begin{array}{c}\alpha ,\beta =\\ \{x,y,z\}\end{array}}\sum _{\begin{array}{c}\omega ,{\omega }^{\prime }=\\ 0,\pm \Lambda \end{array}}\left({\Gamma }_{\alpha \beta ,+}(\omega )+{\Gamma }_{\alpha \beta ,-}({\omega }^{\prime })\right)P_{\beta }(\omega ){\rho }_s{P}_{\alpha }({\omega }^{\prime })$$

$$-\sum _{\begin{array}{c}\alpha ,\beta =\\ \{x,y,z\}\end{array}}\sum _{\begin{array}{c}\omega ,{\omega }^{\prime }=\\ 0,\pm \Lambda \end{array}}\left({\Gamma }_{\alpha \beta ,+}(\omega )P_{\alpha }({\omega }^{\prime })P_{\beta }(\omega ){\rho }_s+h.c.\right),$$

$$P[q_t,q_0]\text{leavevmode}\equiv \text{leavevmode}\mathrm{Tr}\{K_{q_t}U(t,0)K_{q_0}\rho (0)K_{q_0}{U}^{†}(t,0)K_{q_t}\},$$

$$p(Q,t)=\sum _{q_t,q_0}\delta (Q(t)-(q_t-q_0))P[q_t,q_0],$$

$$G(\chi ,t)\text{leavevmode}=\text{leavevmode}\ln \int dQ(t)p(Q,t)e^{i\chi Q(t)}.$$

$$J^{(n)}(t)={\frac{{\partial }^{n}G(\chi ,t)}{\partial (i\chi {)}^n}|}_{\chi =0}$$

Full text

1 Introduction

As modern electrical, optical and mechanical devices[1, 2] continue to shrink in size, nanoscale heat transfer has become an increasingly important research direction. A thorough understanding and characterization of heat dissipation and fluctuation will be critical to maintain the stability of nanoscale devices. For instance, better design of nanostructures to avoid joule heating[3] could be realized with theoretical insights gained from studies of heat transfer in model systems. Furthermore, the advancement of novel technologies such as phononics[4] and quantum heat engines[5, 6, 7] have spurred further interests in acquiring a precise control[8, 9, 10] of heat flows at nanoscale. Beyond the application driven needs, the quantum transport also provides an experimentally accessible platform to explore the rich set of nonequilibrium physical phenomena[11, 12, 13] in the quantum regime. For instance, the fluctuation theorem for charge transport have been experimentally verified in quantum dot systems [14, 15, 16].

Motivated by the aforementioned interests, our group has recently developed a novel approach, based on a combination of the nonequilibrium polaron transformed Redfield equation[17, 18] (NE-PTRE) and the full counting statistics, to study heat transfer and its higher order moments for a finite-size quantum system simultaneously coupled to two heat reservoirs held at different temperatures. The NE-PTRE approach inherits the physically transparent structure of the Redfield equation that facilitates the analyses of the underlying heat-conduction physics. For instance, the parity classified transfer processes[17] through a two-level spin junction has been unraveled within the polaron framework. More importantly, the NE-PTRE addresses the shortcomings of the standard Redfield equation approach which perturbatively treats the system-reservoir interactions in the weak coupling limit. In Ref. 17, it has been demonstrated that the NE-PTRE provides an analytical expression for heat current interpolating accurately from the weak to strong coupling regimes for a super-Ohmic bath model. This successful unification of heat current calculations builds upon the physical picture that bath modes displace to new stable configurations under coupling to the system. However, the polaron technique[19, 20, 21] performs rather poorly with the slower bath modes as they are sluggish and fail to dress the system. Hence, potential problems arise when the polaron technique is applied to Ohmic and sub-Ohmic bath models which feature more prominent contributions from the low-frequency modes (i.e., the infrared divergence). Indeed, it can be easily shown that the NE-PTRE reduces to the Nonequilibrium Non-Interacting Blip Approximation[22] (NE-NIBA) results and captures only the incoherent part of the heat transfer processes for Ohmic baths.

Inspired by earlier studies on open quantum systems, the variational polaron transformation[23, 24] is adopted in this work to extend the applicability of the NE-PTRE beyond the super-Ohmic bath models. Under the variational treatment, slow modes are displaced with reduced displacements. Technically, this variational modulation avoids the infrared divergence of a full polaron transformation when applied to Ohmic baths. The optimized dispaclements are determined by minimizing the upper bound of an effective free energy based on the Feynman-Bogoliubov inequality. This variational ansatz only holds strictly for the equilibrium systems. In this work, we extend this equilibrium technique to the nonequilibrium domain. In subsequent discussions, we denote the improved method nonequilibrium variational polaron transformed Redfield equation (NE-VPTRE). As confirmed by the numerical study, the NE-VPTRE method provides significantly more accurate results than that of the original NE-PTRE. In particular, the polaron-based calculation of heat current is now pushed beyond the NIBA accuracy for Ohmic and sub-Hhmic bath models. This achievement enables a unified heat current calculation (from weak to strong couplings) for all spectral densities of the thermal bath.

The paper is organized as follows. In Sec. II, we discuss the formal structure of nonequilibrium steady state and define the effective temperature for quantum transport, and then we introduce the non-equlibrium spin-boston (NESB) model and outline the derivation of NE-VPTRE. In Sec. III, we analyze the proposed nonequilibrium variational method, bechmark the theoretical predictions with the numerical exact results for NESB models, and discuss the considerable advantages of the variational polaron transformation over the usual polaron transformation. Finally, we present a brief summary in Sec. IV.

2 Method

2.1 Nonequilibrium Steady State

We illustrate the general non-equilibrium transport setup with the example of heat transfer, as illustrated in Fig. 1. Here, an open system couples directly to two thermodynamic reservoirs held at different temperatures, and the long-time limit of the system defines the nonequilibrium steady state (NESS), which is characterized by the heat flux from the hot to cold baths. This nonequilibrium phenomenon can be described by a general Hamiltonian,

[TABLE]

where $H_{s}$ is the open system, $H^{R}_{B}$ and $H^{L}_{B}$ denote the left and right bosonic reservoirs, and $H_{I}$ describes the system-reservoir interaction. The two heat reservoirs are held fixed at the inverse temperature $\beta_{R}$ and $\beta_{L}$ respectively. This setup encompasses a broad range of dissipative and transport settings. In this work, we set $\hbar=1$ and $k_{B}=1$.

Our approach is built upon the theoretical concept[25, 26] that claims the formally exact nonequilibrium steady state (NESS), such as the heat transfer model with a temperature bias, can be cast in a Gibbs-like expression. This effective global equilibrium state (for the open system and its two heat baths) is augmented with an extra term, often referred to as the $Y$ operator, which corresponds to the entropy production associated with the irreversible transport such as the heat transfer driven by a temperature gradient. This formal structure reads

[TABLE]

where $\bar{\beta}$ is an effective temperature and

[TABLE]

$\Omega=\lim_{t\rightarrow\infty}\exp(itH)\exp(-itH_{0})$ and $\delta\beta=\beta_{L}-\beta_{R}$.

Equation (2) is a joint density matrix for the open system and its environment. The close resemblance of Eq. (7) to a thermal equilibrium state implies that various nonequilibrium quantities can be computed by adopting the standard equilibrium techniques for observable calculations, provided one has access to the exact form of the $Y$ operator. In general, this is an extremely difficult task except for the simplest models. In this work, we simply exploitH the knowledge of the formal structure[27, 28, 26] of NESS and devise a variational polaron transformation technique for quantum transport problems.

2.2 Mapping of Non-equilibrium Spin-Boson Model

For convenience, we first introduce the non-equlibrium spin-boson (NESB) model, which is used throughout the paper to formulate and calibrate the NE-VPTRE approach. The model Hamiltonian[29, 30, 31, 22, 32] is given by

[TABLE]

On the first line, the Hamiltonians $H_{0}=H_{s}+H_{B}$ and $H_{I}$ describe the uncoupled system $(s)$ and bath $(B)$ and their coupling, respectively. On the third line, we specify the system as a two-level spin with $\sigma^{z/x}$ referring to the standard Pauli matrices. $H_{B}=\sum_{v=L,R}H_{B}^{v}=\sum_{k,v}\omega_{k,v}b^{\dagger}_{k,v}b_{k,v}$ denotes the left and right bosonic reservoirs with $b_{k,v}^{\dagger}$ and $b_{k,v}$ the creation and annihilation operators for the $k$-th mode in the $v$-th bath. The two heat reservoirs are held fixed at the inverse temperature $\beta_{v}$, respectively. For the NESB model, the reservoirs are linearly coupled to the system as $B_{v}=\sum_{k}g_{k,v}(b_{k,v}^{\dagger}+b_{k,v})$. The model is summarized in Fig. 1. Though formulated in the context of the NESB model, we shall emphasize that the subsequent derivations of the effective temperature in Eq. (7) and NE-VPTRE are completely general and are not limited to the model system.

The influences of bosonic reservoirs on the system are succinctly encoded in the spectral density

[TABLE]

In the continuum limit, the spectral density can be assumed to take on the standard form

[TABLE]

with $\alpha_{v}$, the dimensionless system-reservoir coupling strength, and $\omega_{c,v}$, the cutoff frequency of the $v$–th bosonic reservoir. In this work, we will focus on the Ohmic case, $s=1$, which features more prominent contributions from the low–frequency modes in comparison to the super-Ohmic models with $s=3$ analyzed in earlier works using NE-PTRE method. We will also consider different cutoff functions $\psi_{c}\left(\left.\omega\right/\omega_{c,v}\right)$ to illustrate the robustness of NE-VPTRE in providing consistently improved results on heat current calculations under a variety of model specifications.

A challenging question in studying non-equilibrium transport is the definition of the effective temperature of the system. In the non-equilibrium set-up, the temperatures of the reservoirs are well-defined, but the temperature of the system is not directly determined or even defined. This posts a conceptual difficulty in extending variational polaron method to a non-equlibrium setup. To address this question, the effective temperature is first introduced within the framework of the nonequilibrium steady state and then justified by means of a simple mapping procedure.

To proceed, we make an ansatz that the NESS can be approximated by the zero-th order term (with respect to $H_{I}$) of Eq. (2)[25],

[TABLE]

where the effective temperature is given

[TABLE]

where $\alpha_{\rm L}$ and $\alpha_{\rm R}$ are the coupling strengths to the left and right baths, defined in Eq. (6). The choice of the effective temperature $T_{eff}$ is derived next in Sec. 2.3 and also justified in the next paragraph. We note that Eq. (7) preserves the factorized form of the density matrix due to the absence of $H_{I}$. At this point, we have turned the nonequilibrium setup to an approximate but much more familiar equilibrium setting in order to derive variational polaron transformed master equation under the effective temperature $K_{B}T_{eff}=1/\bar{\beta}$ with the canonical form of equilibrium density matrices.

Alternatively , we can adopt a simple mapping procedure to determine the effective temperature. In the non-equilibrium setup illustrated in Fig. 1, the thermal effect of a reservoir on the system is characterized by the influence functional

[TABLE]

where $\gamma_{v}(\omega)$ is the spectral density. We now assume the two thermal reservoirs couple to the system with the same functional form of spectral density but different coupling strengths and temperatures. Then, the two reservoirs can be combined to a single effective reservoir characterized by the total influence functional $C_{eff}(t)=C_{R}(t)+C_{L}(t)$. When the temperature difference is small ($\delta T/\tilde{T}\ll 1$) and/or the temperatures are high ($\beta h\omega\ll 1$), we can approximately write

[TABLE]

where the effective temperature is given exactly as in Eq. (8). Thus, we can map the non-equilibrium setup to an effective equilibrium setup such that the system relaxes to the equilibrium distribution characterized by the effective temperature. It is also possible to generalize the above discussion to a more general case where the two coupling operators are not identical, i.e., non-communitive transport.

While we use the NESB model to illustrate our newly introduced approach to the nonequilibrium quantum transport, we emphasize the essential assumption, such as Eq. (7), is based on a rigorous theoretical framework[25] and can be adapted to other related transport models.

For a general non-equilibrium transport set-up, the basic principle of non-equilibrium variational polaron transformation remains applicable, though the calculation can be more involved. To go beyond the NESB model, we can consider a multi-site system, such as an exciton chain or a spin chain.[33] Then the effective temperature is not a constant, but becomes a function of the coordinate along the chain. Often, we can find an analytical solution of the temperature profile in some limiting cases, such as the strong coupling limit and/or Markov limit. Then, we use the effective temperature profile to define the non-equilibrium Gibbs state and apply variational polaron evaluation. Though the choice of the reference system is flexible and can affect the accuracy of the theoretical prediction, the NE-VPTRE approach remains valid and general.

2.3 Non-equilibrium Fermi’s Golden-Rule Rate and Effective Temperature

For the NESB model considered in this work, one can write down the exact Liouville equation[34] for the composite system projected in the system’s local basis $\{\left|1\right\rangle,\left|2\right\rangle\}$,

[TABLE]

where the Liouville superoperators are defined as $\mathcal{L}_{1/2}A\equiv[H_{1/2},A]$, $\mathcal{L}_{12}A=H_{1}A-AH_{2}$ and $\rho_{12}=\rho_{21}^{*}$. In this subsection, we note the notation changes: $V=\Delta/2$, $H_{1}=\epsilon/2+\sum_{v}H_{1}^{v}$ with $H_{1}^{v}\equiv H_{B}^{v}+\left|1\right\rangle\left\langle 1\right|B_{v}$ and $H_{2}\equiv\sum_{v}H_{B}^{v}-\left|2\right\rangle\left\langle 2\right|B_{v}$.

Assuming $\rho_{12}(0)=0$, one obtains a formal solution,

[TABLE]

Next we assume the density matrix elements remain factorized as $\rho_{1}(t)=P_{1}(t)\rho_{b,1}$ and $\rho_{2}(t)=P_{2}(t)\rho_{b,2}$ in which the two baths are separately equilibrated with the system such that $\rho_{b,i}\propto\exp(-\beta_{L}H_{i}^{L})\exp(-\beta_{R}H_{i}^{R})$ with $i=1,2$. Substituting Eq. (12) into the last two equations in Eq. (2.3) and use the factorization assumption to trace out the bath, we obtain a general Fermi golden rule rate equation

[TABLE]

Imposing the Markov approximation and keep only the first order expansion of the rate kernel $K_{1/2}$, one obtains the standard rate equation,

[TABLE]

with the lineshape function given by

[TABLE]

The above derivation follows closely the formulation of non-Markov quantum rate equation in Ref. [34].

If one further takes a short time and high temperature expansion on all trigonometric functions and $\coth(\beta_{v}\omega/2)$ for $g_{v}(t)$, the rate constants can be obtained in generalized Marcus form[35] after performing a Gaussian integration.

[TABLE]

with $T_{eff}=\frac{\alpha_{L}}{\alpha}T_{L}+\frac{\alpha_{R}}{\alpha}T_{R}$, and $\bar{\beta}=1/T_{eff}$. Since the rate constant is controlled by the effective temperature, this motivates us to propose the ansatz that the system is characterized by the effective temperature in Eq. (8).

2.4 The Nonequilibrium Variational Polaron Transform

We now apply the variational polaron transformation to the NESB model. The generalized polaron displacement operator is given by

[TABLE]

The transformed Hamiltonian $\tilde{H}=UHU^{\dagger}=\tilde{H}_{s}+\tilde{H}_{I}+H_{b}$ with

[TABLE]

The rotated system-bath interactions take on the form,

[TABLE]

where the displacement operator $D$ was defined in Eq. (17). The renormalized tunneling matrix element $\Delta_{R}=\Delta\eta$ reflects the polaron dressing effects and $\eta$ is given in Eq. (22).

The displacement parameters, $\{f_{k,v}\}$, are determined by minimizing the “effective” free energy upper bound $A_{B}$ given by the Feynman-Bogoliubov inequality,

[TABLE]

where the second term $\langle\tilde{H}_{I}\rangle_{0}=0$ by construction and higher-order terms at the end of right hand side are ignored. If we write $f_{k,v}=g_{k,v}F_{v}(\omega_{k,v})$, then the minimization conditions, $\partial A_{B}/\partial f_{k,v}=0$, leads to the set of self-consistent equations consisting of,

[TABLE]

where $\Lambda=\sqrt{\epsilon^{2}+\Delta_{R}^{2}}$ and the tunneling renormalization factor reads

[TABLE]

While our discussion uses the nonequilibrium spin-boson (NESB) model[29, 30] for illustration, we emphasize Eqs. (7)-(21) can be generalized for other models.

Following the derivation[36] of the Born-Markovian Redfield equation in the polaron picture[37, 38, 20, 17, 21], one obtains

[TABLE]

where the transition rates $\Gamma_{\alpha\beta,\pm}(\omega)$ are provided explicitly in Supplementary Materials and the eigen-basis decomposition of Pauli matrices in the interaction picture gives $\sigma^{\alpha}(-\tau)=\sum_{\omega=0,\pm\Lambda}P_{\alpha}(\omega)e^{i\omega\tau}$.

It has been observed that the variational polaron transformation may suffer from sharp changes in the renormalized tunneling constant,[23, 37, 19] which can lead to difficulties in the prediction of density matrix propagation. However, the numerical results reported in the next section suggest that the VPTRE prediction of steady-state heat current does not suffer from the discontinuity problem. In essence, heat current a non-equilibrium steady-state (NESS) solution and is thus more related to thermal equilibrium in the long-time limit than to dynamic coherence at finite times. In this context, the thermodynamic consistency of PTRE or VPTRE established early[21] can help explain the accuracy of the heat current prediction reported next.

2.5 Heat Current and Full Counting Statistics

The definition of heat current is formalized through a two-time measurement protocol[11, 12]: at time $t=0$, an initial measurement implemented via a projector $K_{q_{0}}=|q_{0}\rangle\langle q_{0}|$ to determine the energy content of $H_{B}^{R}=\sum_{k}\omega_{k,R}b^{\dagger}_{k,R}b_{k,R}$ and get an outcome $q_{0}$. A second measurement is performed at a later time $t>0$ with another projector $K_{q_{t}}=|q_{t}\rangle\langle q_{t}|$ that gives an outcome $q_{t}$. Hence, the net transferred heat is determined by $Q(t)=q_{t}-q_{0}$. The joint probability to measure $q_{0}$ and $q_{t}$ at two specified time points reads

[TABLE]

where $U(t,0)$ the unitary time evolution operator of the total system and $\rho(0)$ is the initial density matrix. One can further define the probability distribution for the net transferred quantity, $Q(t)=q_{t}-q_{0}$ over a period of time $t$,

[TABLE]

where $\delta(x)$ denotes the Dirac distribution. The corresponding cumulant generating function (CGF) for $p(Q,t)$ reads

[TABLE]

with $\chi$ the counting-field parameter. From CGF $G(\chi,t)$, one obtains an arbitrary $n$-th order cumulant of heat transfer via

[TABLE]

Notably, the first and second cumulants correspond to the heat current and its noise power.

If the reduced density matrix (RDM) $\rho_{s}$ is decomposed into different subspaces of net transferred $Q$, i.e. $\rho_{s}(t)=\int dQ\rho_{s}(Q,t)$, then $p(Q,t)=\text{Tr}_{s}\rho_{s}(Q,t)$. Substituting the $Q$-resolved RDM into Eq. 26, one may further define $\chi$-resolved RDM,

[TABLE]

The statistics of heat transfer can be extracted from $\rho_{s}(\chi,t)$ via the relation $G(\chi,t)=\ln\text{Tr}_{s}\rho_{s}^{\chi}(t)$.

The dynamics of $\rho^{\chi}_{s}(t)$ is obtained after a counting field $\chi$ is introduced via a unitary transformation $U_{\chi}=\exp(-iH_{b}^{R}\chi/2)$ to count the net amount of transferred energy into the right bath. Without delving into further derivations, which is a straightforward generalization of earlier works, we get a $\chi$-dependent NE-VPTRE,

[TABLE]

By adopting the Liouville-space notation, Eq. (2.5) could be succinctly expressed as $d\rho_{s}^{\chi}(t)/dt=-\mathcal{L}_{\chi}\rho_{s}^{\chi}(t)$. The formal solution assumes the simple form $\rho_{s}^{\chi}(t)=e^{-\mathcal{L}_{\chi}t}\rho_{s}^{\chi}(0)$ when the Hamiltonian is time-independent. The stationary CGF at the steady state is obtained via $G(\chi)=\lim_{t\rightarrow\infty}\frac{1}{t}\ln\text{Tr}_{s}\rho_{s}^{\chi}(t)$. In the asymptotic limit[39, 17], one can show that $G(\chi)=-E_{0}(\chi)$ where $E_{0}(\chi)$ is the ground state of the superoperator $\mathcal{L}_{\chi}$. Hence, the heat current can be conveniently calculated via

[TABLE]

In the case of the unbiased NESB model $\epsilon=0$, we derive the following expression for the heat current

[TABLE]

where the renormalized energy gap defined earlier reduces to $\Lambda=\Delta\eta$ for $\epsilon=0$, and

[TABLE]

with $\Sigma_{d}=\sum_{\nu=y,z}\phi_{\nu\nu}(\Lambda)+\phi_{\nu\nu}(-\Lambda)$, $\Sigma_{od}=\sum_{\nu=y,z}\phi_{\nu\bar{\nu}}(\Lambda)+\phi_{\nu\bar{\nu}}(-\Lambda)$, $\Delta_{d}=\sum_{\nu=y,z}\phi_{\nu\nu}(\Lambda)-\phi_{\nu\nu}(-\Lambda)$ and $\Delta_{od}=\sum_{\nu=y,z}\phi_{\nu\bar{\nu}}(\Lambda)-\phi_{\nu\bar{\nu}}(-\Lambda)$. Note the subscript satisfies $\bar{\nu}=z/y$ if $\nu=y/z$. The correlation functions read,

[TABLE]

where $C_{\nu\mu}(t)$ are explicitly given in the supplementary information, and $C_{\nu\mu}(\omega,\omega^{\prime})$ are defined on the second line of this equation.

In the weak and strong coupling regimes, Eq. (31) reduces smoothly back to the Redfield and NIBA results, respectively. In the weak coupling limit, we note the polaron displacement $F_{v}(\omega)\rightarrow 0$,which subsequently leads to the correlation functions $C_{xx}(0,\omega)\rightarrow 0$ and $C_{yy}(\pm\Lambda,\omega)\rightarrow 0$. In this case, only the two terms, $C_{zz}(\pm\Lambda,\omega)$, survive and approach the Redfield result. On the other hand, in the strong coupling limit, $F_{v}(\omega)\rightarrow 1$ yield a full polaron displacement while the last two terms of Eq. (31) vanish. In this case, Eq. (31) reduces to NIBA results in the strong coupling limit. In the next section, we will numerically demonstrate the limiting behaviors of Eq. (31) in various examples.

3 Results

3.1 Nonequilibrium Variational Polaron Displacement: Cross-bath Correlation

We first take a closer look at the nonequilibrium variational polaron displacements proposed in the Method section. The functional form of $F_{v}(\omega)$ in Eq. (21) is almost identical to the standard equilibrium forms. Hence, many of the equilibrium results carry over to the nonequilibrium set-ups. For instance, $F_{v}(\omega)$ prescribes small displacements, $f_{vk}\approx 0$, for slow modes and close to the full polaron displacements, $f_{vk}=g_{vk}$, for fast modes. Due to additivity nature of the exponents representing contributions from different baths in Eq. (22), the renormalized tunneling factor $\eta$ can be approximately factorized into two components to reflect interaction with two baths, respectively. Because of this resembleance to equilibrium results, we can clearly see that the original NE-PTRE works well for high-temperature and fast-relaxing heat reservoirs when $\eta\ll 1$.

The major distinction that sets apart nonequilibrium variational method is the cross-bath correlation effects mediated by the central spin via the hyperbolic tangent factor, $\tanh\left(\frac{\bar{\beta}\Lambda}{2}\right)$, evaluated at the weighted average of the inverse temperature $\bar{\beta}$ in Eq. (21). In Fig. 2, the displacements of the low-frequency modes in the left bath clearly depend on the parameters of the right bath as long as $\Lambda\neq 0$, i.e. the two system eigenstates are not degenerate. When the right bath is more strongly coupled to the system or is elevated to a higher temperature, the displacements of left bath modes also become more pronounced as displayed in the figure.

We first consider an NESB model with two super-Ohmic baths characterized by a rational cutoff function $\psi_{c}(\omega/\omega_{c,v})=\frac{1}{(1+(\omega^{2}/\omega_{c}^{2}))^{2}}$. As demonstrated in Fig. 3, the new NE-VPTRE results are in excellent agreement with the numerically exact hierarchical equation of motion[40, 41, 42, 43] (HEOM) data than the original NE-PTRE method. In the intermediate coupling strength regime presented in Fig. 3, the full polaron displacement overestimates the heat current. This discrepancy is precisely due to the inaccurate treatment of the low-frequency modes in the baths. This observation can be confirmed by increasing the cutoff frequency $\omega_{c,v}$ of the bath, the dissipations induced by the low-frequency modes are diluted and the variational functions $F_{L/R}(\omega)$, Eq. (21), approach a constant unity. The results of the two polaron methods then converge in this limit.

3.2 Unified Heat Current Calculation for Ohmic Spectral Density

Next, we turn to the NESB cases in which heat reservoirs are featured with Ohmic spectral densities. The NE-PTRE reduces exactly to NIBA regardless of the cutoff function $\psi_{c}(\omega/\omega_{c,v})$. Hence, the full polaron displacement can only capture the incoherent part of the heat current and provide an accurate results only in the strong coupling and/or scaling limits.

The variational method extends the applicability of the polaron method beyond super-Ohmic cases. To demonstrate the improvement, we first calculate the heat currents as function of the coupling strength in Fig. 4. Different spectral cutoff functions are used in the panel a: $\psi_{v}(\omega)=\frac{1}{(1+(\omega^{2}/\omega_{c}^{2}))^{4}}$ and panel b: $\psi_{v}(\omega)=e^{-\omega/\omega_{c}}$, respectively. Over the entire range coupling strength considered, the NE-VPTRE agrees well with the exact result (i.e. NEGF in panel a and HEOM in panel b). In the weak-coupling limit, the exact results approach the Redfield result; while the NE-PTRE method (equivalent to NIBA for Ohmic baths) do not fare well. These two results establish the superiority of NE-VPTRE over the NE-PTRE[17] in dealing with Ohmic baths.

We next investigate the heat current as a function of the spin tunneling frequency $\Delta$ in Fig. 5. The two baths are taken to be an identical Ohmic form, characterized by an exponential cutoff function, but held at different temperatures. A weak system-bath coupling strength is used for the two cases presented in Fig. 5 such that the Redfiled results provide reasonably accurate benchmarks within the specified range of $\Delta/\omega_{c}$. Under both small (panel a) and large (panel b) thermal bias, an excellent agreement between Redfield and NE-VPTRE affirm the applicability of the variational polaron-based framework to calculate heat current beyond the scaling limit, i.e. $\Delta/\omega_{c}\ll 1$. The essential need of adopting a variational polaron displacement is further supported by the observation on how the NE-PTRE (or NIBA) fails to capture the turnover behaviors portrayed in Fig. 5 because it only accounts for the spin tunneling up to $\Delta^{2}$ in the heat current calculation. On the other hand, the variational polaron result takes into account of higher order terms of $\Delta$ in the weak system-bath coupling limit.

The turnover can be more transparently explained via the Redfield expression for heat current,

[TABLE]

where $n_{v}(\omega)$ is the Bose-Einstein distribution. In the Redfield framework, a classical-like sequential energy transfer is at work, and only bath modes in resonance (i.e. $\omega=\Delta$) with the system contribute to the heat conduction. Since we only consider parameter regime $\Delta<\omega_{c}$ for both cases considered in Fig. 5, the Ohmic spectral densities $\gamma_{v}(\Delta)$ are monotonically increasing with respect to $\Delta$. The turnovers of the heat current in Fig. 5 are entirely controlled by the Bose-Einstein distributions in Eq. (34). In short, the current diminishes once $\Delta$ exceeds the thermal energy appreciably such that the resonant modes in the bath are not thermally excited. The shift of the peak current to higher value of $\Delta$ (for the higher temperature case in Fig. 5b) also supports the claim that the heat conduction depends largely on the temperatures in this case.

3.3 Thermal Rectifications

In this section, we consider an NESB model where the two Ohmic baths have asymmetrical coupling strengths, i.e. $\alpha_{L}\neq\alpha_{R}$. Because the spin is an anharmonic system, a rectification can arise when the temperature bias $T_{L}-T_{R}=\delta T$ on the two asymmetric baths are switched to $T_{L}-T_{R}=-\delta T$. More precisely, we define the rectification ratio $R\equiv\left|J_{R}(\delta T)/J_{R}(-\delta T)\right|$ under the constraint of a fixed $\alpha=\alpha_{L}+\alpha_{R}$ to remove the dependence of the rectification ratio $R$ on the overall magnitude of the coupling strengths, $\alpha$. The rectification ratio is a useful indicator to refelect whether the central quantum system is an ideal thermal diode[44, 45].

In Fig. 6, we consider two overall magnitudes of the coupling strengths: $\alpha=0.05$ and $\alpha=0.2$. For both cases, the same Redfield result (green curve in Fig. 6) is obtained. This seemingly universal rectification behavior is an artifact of approximations and can be seen from Eq. (34). Given the linearity of the current with respect to $\alpha$ in Eq. (34), the rectification ratio R only depends on the ratio of $\alpha_{R}/\alpha_{L}$ and not their overall magnitude $\alpha$. Nevertheless, the Redfield result is reliable for the weak value case such as $\alpha=0.05$. Indeed, the NE-VPTRE result (red cross) also agrees well with the Redfield. For larger coupling case ($\alpha=0.2$), the NE-VPTRE predicts an enhanced rectification ratio as typically expected for an anharmonic junction. On the other hand, the two NIBA results (blue circle and dashed line corresponding to $\alpha=0.05$ and $0.20$, respectively) also collapse onto the same line in Fig. 6. This is mainly due to the set of parameters chosen for illustration. In particular, $\omega_{c}=10\Delta$ is at the borderline of the scaling limit. When applying NIBA outside their valid parameter regimes, the deviation from accurate results could be significant as manifested by the rectification calculation shown in Fig. 6 as well as Fig. 5.

In this case, the superiority of NE-VPTRE is attributed to its better handling of the low-frequency modes through reduced polaron displacements. When $\omega_{c}$ is further increased, the overall contributions from the low-frequency modes is diluted, and one would find NIBA results on $R$ to agree better with that of NE-VPTRE.

4 Discussion

In summary, we extend the nonequilibrium polaron transformed Redfield equation (NE-PTRE) framework by variationally tuning the polaron displacements. The generalization of a free-energy-based variational principle and Feynman-Bogoliubov inequality is built upon the ansatz, Eq. (7), that the zero-th order nonequilibrium steady state of the composite system assumes a form resembling an equilibrium density matrix. Similar to the equilibrium cases, the variational method extends the usefulness of a polaron picture beyond the nonadiabatic limit for nonequilibrium processes. This achievement allows us to formulate a heat transfer theory beyond the super-Ohmic models in the polaron picture.

In the Result section, we explicitly demonstrate several improvements of the newly proposed NE-VPTRE method in the aforementioned circumstances. Specifically, we observe that (1) improved numerical accuracy for super-Ohmic bath models over a broader range of model parameters. (2) Correctly recover the coherent tunneling effects on the heat current for the Ohmic bath models beyond the scaling limit as illustrated in Fig. 5. (3) For calculations of rectification, an important indicator for the quality of a thermal diode, NE-VPTRE provides absolute improvements over both Redfield and the original NE-PTRE (equivalent to NIBA) methods in the Ohmic cases.

The present work not only significantly extends the applicability of a transparent and numerically efficient polaron picture to calculate heat current but also builds a foundation to calculate higher order statistics[18, 13] of heat transfer, such as the noise power spectrum relating to the fluctuations of heat transfer. Future work will generalize the present framework to handle non-Markovian effects and to formulate a unified energy transfer theory incorporating the low-temperature regime and to investigate fluctuations.

Acknowledgements

All authors acknowledge the support from Singapore-MIT Alliance for Research and Technology (SMART). J. C. acknowledges the support from the National Science Foundation (NSF) (Grant CHE 1836913 and CHE 1800301).