Heat transfer statistics in mixed quantum-classical systems

TL;DR

This paper develops a mixed quantum-classical framework to analyze heat transfer at the nanoscale, deriving general statistical expressions and demonstrating their application to a spin boson model.

Contribution

It introduces a novel mixed quantum-classical approach using QCLE and full counting statistics to study heat transfer without simplifying bath assumptions.

Findings

Derived a general expression for the heat moment generating function.

Showed the steady state fluctuation symmetry holds up to order ℏ.

Simulated heat transfer in the nonequilibrium spin boson model.

Abstract

The modelling of quantum heat transfer processes at the nanoscale is crucial for the development of energy harvesting and molecular electronics devices. Herein, we adopt a mixed quantum-classical description of a device, in which the open subsystem of interest is treated quantum mechanically and the surrounding heat baths are treated in a classical-like fashion. By introducing such a mixed quantum-classical description of the composite system, one is able to study the heat transfer between the subsystem and bath from a closed system point of view, thereby avoiding simplifying assumptions related to the bath time scale and subsystem-bath coupling strength. In particular, we adopt the full counting statistics approach to derive a general expression for the moment generating function of heat in systems whose dynamics are described by the quantum-classical Liouville equation (QCLE). From this expression, one can deduce expressions for the dynamics of the average heat and heat current, which may be evaluated using numerical simulations. Due to the approximate nature of the QCLE, we also find that the steady state fluctuation symmetry holds up to order $\hbar$ for systems whose subsystembath couplings and baths go beyond bilinear and harmonic, respectively. To demonstrate the approach, we consider the nonequilibrium spin boson model and simulate its time-dependent average heat and heat current under various conditions.

Equations

$$\hat{H}={\hat{H}}_S(\hat{\mathbf{x}})+\sum _{v=1}^K{\hat{H}}_B^v({\hat{\mathbf{X}}}_v)+{\hat{H}}_I(\hat{\mathbf{x}},\hat{\mathbf{X}}),$$

$$⟨Q_v(t)⟩=⟨{\hat{H}}_B^v(0)-{\hat{H}}_B^v(t)⟩,$$

$$⟨J_v(t)⟩=-\frac{d}{dt}⟨{\hat{H}}_B^v(t)⟩.$$

$$⟨Q_v(t)⟩$$

$$⟨J_v(t)⟩$$

$$P(a_t,a_0)=\mathrm{Tr}\{{\hat{\mathcal{P}}}_{a_t}\hat{U}(t,0){\hat{\mathcal{P}}}_{a_0}{\hat{\rho }}_0{\hat{\mathcal{P}}}_{a_0}{U}^{†}(t,0){\hat{\mathcal{P}}}_{a_t}\},$$

$$P(a_t,a_0)=\mathrm{Tr}\{{\hat{\rho }}_0{\hat{\mathcal{P}}}_{a_0}{\hat{U}}^{†}(t,0){\hat{\mathcal{P}}}_{a_t}\hat{U}(t,0)\}.$$

$$P_t(\Delta a)=\sum _{a_t,a_0}\delta (\Delta a-(a_0-a_t))P(a_t,a_0).$$

$$Z({\chi }_v,t)\equiv \int d\Delta a{e}^{i{\chi }_v\Delta a}{P}_t(\Delta a),$$

$$Z({\chi }_v,t)=\sum _{a_t,a_0}{e}^{-i{\chi }_v(a_t-a_0)}P(a_t,a_0).$$

$$Z({\chi }_v,t)=\mathrm{Tr}[e^{i{\chi }_v{\hat{H}}_B^v}{e}^{-i{\chi }_v{\hat{H}}_B^v(t)}{\hat{\rho }}_0],$$

$$Z(\{{\chi }_v\},t)=\mathrm{Tr}\left[e^{i\sum _v{\chi }_v{\hat{H}}_B^v}{e}^{-i\sum _v{\chi }_v{\hat{H}}_B^v(t)}{\hat{\rho }}_0\right],$$

$$⟨Q_v^n(t)⟩={\frac{{\partial }^n}{\partial (i{\chi }_v{)}^n}Z(\{{\chi }_v\},t)|}_{\{{\chi }_v\}=0}.$$

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

$$⟨\mathcal{R}-\frac{\mathcal{Z}}{2}\left|\hat{O}\right|\mathcal{R}+\frac{\mathcal{Z}}{2}⟩$$

$$=\frac{1}{(2\pi \hslash {)}^{\mu }}\int d\mathcal{P}{e}^{-(i/\hslash )\mathcal{P}\cdot \mathcal{Z}}{O}_W(\mathcal{R},\mathcal{P}),$$

$$D({\mathcal{X}}_1,{\mathcal{X}}_2,t)$$

$$\frac{\partial }{\partial t}D(t)$$

$$O_W(\mathcal{X})=\int d\mathbf{z}{e}^{(i/\hslash )\mathbf{p}\cdot \mathbf{z}}⟨\mathbf{r}-\frac{\mathbf{z}}{2}\left|{\hat{O}}_W(\mathbf{X})\right|\mathbf{r}+\frac{\mathbf{z}}{2}⟩,$$

$$O_W(\mathcal{X})=O_W(\mathbf{X}),$$

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

$$\overline{D}({\mathbf{X}}_1,{\mathbf{X}}_2,t)=\int d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}d{\mathbf{p}}_{1}d{\mathbf{p}}_{2}D({\mathcal{X}}_1,{\mathcal{X}}_2,t)$$

$$=\frac{1}{(2\pi \hslash {)}^{2{\mu }_B}}\int d{\mathbf{r}}_{1}d{\mathbf{Z}}_{1}d{\mathbf{Z}}_2{e}^{-(i/\hslash )({\mathbf{P}}_1\cdot {\mathbf{Z}}_1+{\mathbf{P}}_2\cdot {\mathbf{Z}}_2)}$$

$$\times ⟨{\mathbf{r}}_1|⟨{\mathbf{R}}_1+\frac{{\mathbf{Z}}_1}{2}\left|e^{(i/\hslash )\hat{H}t}\right|{\mathbf{R}}_2-\frac{{\mathbf{Z}}_2}{2}⟩$$

$$\times ⟨{\mathbf{R}}_2+\frac{{\mathbf{Z}}_2}{2}\left|e^{-(i/\hslash )\hat{H}t}{\hat{\rho }}_0\right|{\mathbf{R}}_1-\frac{{\mathbf{Z}}_1}{2}⟩|{\mathbf{r}}_1⟩.$$

$$\frac{\partial }{\partial t}\overline{D}(t)$$

$$\frac{\partial }{\partial t}{\overline{D}}_{QC}(t)$$

$${\overline{D}}_{QC}({\mathbf{X}}_1,{\mathbf{X}}_2,t)=e^{-i\mathcal{L}({\mathbf{X}}_2)t}\overline{D}({\mathbf{X}}_1,{\mathbf{X}}_2).$$

$$Z_{QC}(\{{\chi }_v\},t)$$

$${\overline{D}}^{{\alpha }_1{\alpha }_2{\alpha }_2^{\prime }}({\mathbf{X}}_1,{\mathbf{X}}_2)=\frac{1}{(2\pi \hslash {)}^{2{\mu }_B}}\int d{\mathbf{Z}}_{1}d{\mathbf{Z}}_2$$

Full text

I Introduction

Owing to the rapid development of nanotechnologies in recent decades, heat transfer at the nanoscale has attracted significant attention. Numerous studies have been dedicated to gaining a deep understanding and precise control of the heat transfer, which has impacts at both the fundamental and practical levels. So far, heat transfer has been studied in small and well-characterized quantum systems. On the experimental side, systems such as molecular junctions can be fabricated in the laboratory Wang et al. (2007); Schwab et al. (2000); Carter et al. (2009); Losego et al. (2012); Meier et al. (2014); Cui et al. (2017), while on the theoretical side, simplified models can be put forward and studied with a host of fully quantum methods Segal and Nitzan (2005); Velizhanin et al. (2008); Ren et al. (2010); Nicolin and Segal (2011a, b); Ruokola and Ojanen (2011); Segal (2013); Saito and Kato (2013); Yang and Wu (2014); Wang et al. (2015a); Carrega et al. (2016); Liu et al. (2017); Wang et al. (2017).

When the heat transfer occurs in a complex, many-body system such as a molecular aggregate May and Kühn (2011) or a self-assembled monolayer junction Majumdar et al. (2015), which may not be well described in terms of a simplified model containing a small number of degrees of freedom (DOF), a fully quantum approach to modelling the heat transfer dynamics will be computationally intractable. In this case, an approximate treatment of the dynamics is required to gain insight into the system under study. Mixed quantum-classical dynamics methods, which treat a set of light particles of interest (i.e., subsystem) quantum mechanically and the remaining particles in the system (i.e., bath or environment) in a classical-like fashion, provide tremendous computational advantages over fully quantum methods Tully (1990); Prezhdo and Kisil (1997); Martens and Fang (1997); Tully (1998); Kapral and Ciccotti (1999); Wan and Schofield (2000); Horenko et al. (2002); Kelly and Markland (2013); Bai et al. (2014); Kim and Rhee (2014); Wang et al. (2015b); Martens (2016); Wang et al. (2016); Agostini et al. (2016); Subotnik et al. (2016).

In this work, we adopt a mixed quantum-classical approach to modelling heat transfer dynamics that is based on the quantum-classical Liouville equation (QCLE) Aleksandrov (1981); Gerasimenko (1982); Zhang and Balescu (1988); Kapral and Ciccotti (1999), which stems from a linearization of the quantum Liouville equation expressed in the partial Wigner representation Wigner (1932), viz., a description of the subsystem and bath DOF in terms of operators and phase space variables, respectively. The QCLE is chosen as the starting point for our work because (i) several of the popular mixed quantum-classical methods may be derived from this equation Kapral (2015, 2016), and (ii) it yields the exact quantum dynamics for quantum subsystems that are bilinearly coupled to harmonic environments MacKernan et al. (2002), which are frequently used as models for studying energy transfer at the nanoscale. In particular, we combine the QCLE and full counting statistics (FCS) Levitov and Lesovik (1993); Levitov et al. (1996); Belzig and Nazarov (2001); Klich (2003); Bagrets and Nazarov (2003); Pilgram et al. (2003); Saito and Utsumi (2008); Gutman et al. (2010) approaches to derive a general expression for the moment generating function (MGF) of heat, which may then be used to compute the time-dependent average heat and its fluctuations in a system. As the QCLE treats the dynamics of the heat baths explicitly, one can start from the exact definition of the MGF in FCS and does not need to impose any constraints on the bath timescale and subsystem-bath coupling strength, in contrast to the conventional Redfield master equation Segal and Nitzan (2005) and nonequilibrium Green’s function method Liu et al. (2017). Thus, one can apply this combined approach to a wide range of parameter regimes.

Because heat fluctuates at the nanoscale, its average is insufficient to fully characterize a heat transfer process. For a fully quantum system at steady state, heat fluctuations are governed by the steady state fluctuation symmetry (SSFS) of the MGF Esposito et al. (2009); Campisi et al. (2011); Nicolin and Segal (2011a). However, when the dynamics of a fully quantum system is approximated, the behavior of the heat fluctuations may be altered and, as a result, the SSFS may not be satisfied. A direct consequence of this is the breakdown of the fluctuation-dissipation theorem in the linear response regime. Thus, it is of vital importance to assess to what extent the SSFS holds in systems whose dynamics are described by the QCLE. In the case of systems for which QCLE dynamics is exact (e.g., subsystems that are bilinearly coupled to harmonic environments), one expects the SSFS to be strictly preserved, while in the case of systems for which QCLE dynamics is approximate, one expects to reach an approximate nonequilibrium steady state. Nevertheless, in the limit of high temperature and a very small mass ratio between the subsystem and bath particles, the approximations introduced by the QCLE dynamics are expected to be minor.

To illustrate the utility of our approach, we consider the nonequilibrium spin-boson (NESB) model, a prototypical model in the study of quantum energy transfer over the last decade Boudjada and Segal (2014). In particular, we compute the time-dependent average heat and heat current using a recently proposed method for solving the QCLE Liu and Hanna (2018). This method deterministically propagates the dynamics of the system by numerically solving a set of coupled first-order differential equations for the subsystem and bath coordinates. Given its demonstrated accuracy and efficiency in several prototype systems, we believe that a QCLE-based approach to heat transfer statistics will provide a viable way of studying more realistic models of many-body systems.

The paper is organized as follows. We describe the model and MGF of heat in section II. In section III, we derive a general expression for the MGF of heat in the quantum-classical limit. In section IV, we address the question of the extent to which the SSFS holds in systems whose dynamics are described by the QCLE. In section V, we apply our formalism to the NESB model and present and discuss our numerical results for the time-dependent heat and heat current. We summarize our findings in section VI.

II General background

II.1 Model

We consider a composite quantum system in which a subsystem is in contact with $K$ ($K\geq 2$) bosonic heat baths at different temperatures and whose Hamiltonian is given by

[TABLE]

where $\hat{H}_{S}$ is the subsystem Hamiltonian; $\hat{H}_{B}^{v}=\sum_{j=1}^{N_{v}}[\hat{P}_{j,v}^{2}/2+\omega_{j,v}^{2}\hat{R}_{j,v}^{2}/2]$ is the Hamiltonian of the $v$th heat bath at inverse temperature $\beta_{v}$ with $\hat{P}_{j,v}$, $\hat{R}_{j,v}$, and $\omega_{j,v}$ the mass-weighted momentum, position, and frequency of the $j$th oscillator, respectively; and $\hat{H}_{I}$ is the subsystem-bath interaction Hamiltonian with $\boldsymbol{\hat{X}}=(\boldsymbol{\hat{X}}_{1},\boldsymbol{\hat{X}}_{2},\ldots,\boldsymbol{\hat{X}}_{K})$. In the above equation, $\boldsymbol{\hat{x}}=(\boldsymbol{\hat{r}},\boldsymbol{\hat{p}})$ and $\boldsymbol{\hat{X}}_{v}=(\boldsymbol{\hat{R}}_{v},\boldsymbol{\hat{P}}_{v})$ with $\boldsymbol{\hat{R}}_{v}=(\hat{R}_{1,v},\hat{R}_{2,v},\ldots,\hat{R}_{N_{v},v})$ and $\boldsymbol{\hat{P}}_{v}=(\hat{P}_{1,v},\hat{P}_{2,v},\ldots,\hat{P}_{N_{v},v})$, where $N_{v}$ is the number of harmonic oscillators in the $v$th heat bath. In what follows, we assume factorized initial density operators $\hat{\rho}_{0}=\hat{\rho}_{S}(0)\otimes\hat{\rho}_{B}(0)$, where $\hat{\rho}_{S}(0)$ is the initial subsystem density operator and $\hat{\rho}_{B}(0)=\hat{\rho}_{B}^{1}(0)\otimes\cdots\otimes\hat{\rho}_{B}^{K}(0)$ is the initial bath density operator with each $\hat{\rho}_{B}^{v}(0)\propto e^{-\beta_{v}\hat{H}_{B}^{v}}$ assuming a canonical form.

To quantify the heat transfer between the subsystem and its heat baths, we define the average of the heat transferred from the $v$th heat bath to the subsystem as the average change in the bath energy during a time interval $[0,t]$ Agarwalla et al. (2012)

[TABLE]

where the time dependence should be understood in the Heisenberg picture. It follows that the average heat current from the $v$th heat bath to the subsystem may be obtained by taking the time derivative of the above equation, i.e.,

[TABLE]

In anticipation for a mixed quantum-classical description of the system’s dynamics, we express Eqs. (2) and (3) in the partial Wigner representation by taking the Wigner transform Wigner (1932) of these equations over the bath degrees of freedom (DOFs). For a general operator $\hat{A}(t)$, its expectation value in this representation is given by $\langle\hat{A}(t)\rangle=\sum_{\alpha\alpha^{\prime}}\int d\boldsymbol{X}A_{W}^{\alpha\alpha^{\prime}}(\boldsymbol{X},t)\rho_{W}^{\alpha^{\prime}\alpha}(\boldsymbol{X},0)$, where $\{|\alpha\rangle\}$ denotes a complete set of basis states that span the Hilbert space of the subsystem, $(\cdot)_{W}^{\alpha\alpha^{\prime}}\equiv\langle\alpha|(\cdot)_{W}|\alpha^{\prime}\rangle$, and $\hat{\rho}_{W}(0)$ is the partial Wigner transform of $\hat{\rho}_{0}$ Sergi et al. (2003). Using this result, one can directly write Eqs. (2) and (3) in the partial Wigner representation as

[TABLE]

In Eq. (4), the delta function results from the fact that the subsystem and bath are uncorrelated initially. However, at finite times, one must consider the matrix elements of the bath Hamiltonian in the subsystem basis because they depend on the subsystem operators due to the subsystem-bath interaction.

II.2 Moment generating function of heat

To fully characterize a heat transfer process at the nanoscale, not only is information about the average heat and heat current important, but one should also consider the higher order heat fluctuations. The FCS approach Levitov and Lesovik (1993); Levitov et al. (1996); Belzig and Nazarov (2001); Klich (2003); Bagrets and Nazarov (2003); Pilgram et al. (2003); Saito and Utsumi (2008); Gutman et al. (2010) provides a general route for obtaining such statistics of heat in open quantum systems. Recalling that the heat transferred from the $v$th bath to the subsystem in a time interval $t$ is given by $\langle Q_{v}(t)\rangle=\langle\hat{H}^{v}_{B}(0)-\hat{H}^{v}_{B}(t)\rangle$, one may evaluate $\langle Q_{v}(t)\rangle$ using a two-time measurement Campisi et al. (2011) in which the instantaneous eigenvalues (eigenvectors) of $\hat{H}_{B}^{v}$ at time $t$ are $a_{t}$ ($|a_{t}\rangle$). This two-time measurement can be described in terms of the joint probability of measuring $a_{0}$ at time zero and $a_{t}$ at time $t$, i.e.,

[TABLE]

where $\hat{\mathcal{P}}_{a_{t}}=|a_{t}\rangle\langle a_{t}|$, $\hat{U}(t,0)$ is the time evolution operator governed by the total Hamiltonian $\hat{H}$, and $\hat{\rho}_{0}$ is the initial total density operator. Since it has been previously shown that only the part of $\hat{\rho}_{0}$ that commutes with $\hat{H}_{B}^{v}$ determines the moment generating function Esposito et al. (2009), for convenience, we choose $\hat{\rho}_{0}$ such that $[\hat{\rho}_{0},\hat{H}_{B}^{v}]=0$ at time zero (which is the case for factorized initial states). Furthermore, since $[\hat{\mathcal{P}}_{a_{0}},\hat{\rho}_{0}]=0$ (as a result of $[\hat{\rho}_{0},\hat{H}_{B}^{v}]=0$), Eq. (6) becomes

[TABLE]

The probability distribution for the difference $\Delta a$ (i.e., the amount of heat transferred from the measured bath to the subsystem) between the output of the two aforementioned measurements is given by

[TABLE]

The corresponding MGF may be defined as

[TABLE]

where $\chi_{v}$ is the counting field associated with the measurement on the $v$th bath. Upon substituting Eq. (8) into Eq. (9), the MGF becomes

[TABLE]

Noting that $f(\hat{B})=\sum_{b}\hat{\mathcal{P}}_{b}f(b)$, where $f$ is an arbitrary function of an arbitrary operator $\hat{B}$ with $\hat{B}|b\rangle=b|b\rangle$ and $\hat{\mathcal{P}}_{b}=|b\rangle\langle b|$, and substituting Eq. (7) into Eq. (10), the MGF simplifies to

[TABLE]

where $\hat{H}_{B}^{v}(t)=\hat{U}^{\dagger}(t,0)\hat{H}_{B}^{v}\hat{U}(t,0)$. Generalizing this expression to the multiple bath case leads to

[TABLE]

where $\chi_{v}$ is the counting field for the $v$th heat bath, $\{\chi_{v}\}\equiv\{\chi_{1},\chi_{2},\ldots,\chi_{K}\}$, and the trace is performed over all DOFs. It should be noted that, to arrive at this expression, one requires that $[\hat{\rho}_{0},\hat{H}_{B}^{v}]=0$ Esposito et al. (2009), which is the case for the factorized initial state $\hat{\rho}_{0}$.

By differentiating the MGF with respect to the counting field and evaluating the result at $\chi_{v}=0$, one obtains the $n$th moment of heat for the $v$th bath, i.e.,

[TABLE]

As $\langle Q_{v}(t)\rangle$ corresponds to the transferred energy from the $v$th bath to the subsystem during the time interval $[0,t]$, the time derivative of the first moment will give rise to the time-dependent energy current. Higher moments contain information about higher order correlations of the transferred energy.

III Quantum-classical limit of MGF

III.1 Derivation of the MGF

To obtain the quantum-classical limit of the MGF, we start by introducing a coordinate representation $\{\boldsymbol{\mathcal{Q}}\}=\{\boldsymbol{r},\boldsymbol{R}\}$ (calligraphic symbols are used to denote variables for the entire system) into Eq. (12)

[TABLE]

where $\hat{U}(t,0)$ is the time evolution operator governed by the total Hamiltonian $\hat{H}$. We next make a change of variables, $\boldsymbol{\mathcal{Q}}_{1}=\boldsymbol{\mathcal{R}}_{1}-\boldsymbol{\mathcal{Z}}_{1}/2$, $\boldsymbol{\mathcal{Q}}_{2}=\boldsymbol{\mathcal{R}}_{1}+\boldsymbol{\mathcal{Z}}_{1}/2$, $\boldsymbol{\mathcal{Q}}_{3}=\boldsymbol{\mathcal{R}}_{2}-\boldsymbol{\mathcal{Z}}_{2}/2$, and $\boldsymbol{\mathcal{Q}}_{4}=\boldsymbol{\mathcal{R}}_{2}+\boldsymbol{\mathcal{Z}}_{2}/2$, and rewrite the above equation as

[TABLE]

where we have used the notation $\boldsymbol{\mathcal{R}}=(\boldsymbol{r},\boldsymbol{R})$ and $\boldsymbol{\mathcal{P}}=(\boldsymbol{p},\boldsymbol{P})$ (the lowercase and uppercase symbols refer to the subsystem and bath variables, respectively). To arrive at the above equation, we used the fact that the matrix element of an arbitrary operator $\hat{O}$ may be expressed in terms of its Wigner transform $O_{W}$ as follows

[TABLE]

where $\mu=\mu_{S}+\mu_{B}$ is the coordinate-space dimension of the total system and $\boldsymbol{\mathcal{Z}}=(\boldsymbol{z},\boldsymbol{Z})$. Finally, the time-dependent weight function $D$ has the following form

[TABLE]

where $\boldsymbol{\mathcal{X}}=(\boldsymbol{\mathcal{R}},\boldsymbol{\mathcal{P}})$. It is interesting to note that $D(\boldsymbol{\mathcal{X}}_{1},\boldsymbol{\mathcal{X}}_{2},t)$ has the same structure as the spectral density appearing in previous derivations of transport coefficients for mixed quantum-classical systems Sergi and Kapral (2004); Kim and Kapral (2005a, b), with the only difference being that we consider a factorized initial density operator as opposed to a thermal equilibrium state of the total system. Taking into consideration that $\hat{\rho}_{0}$ and $\hat{H}$ do not commute in general, one can show that $D(\boldsymbol{\mathcal{X}}_{1},\boldsymbol{\mathcal{X}}_{2},t)$ obeys the following equation of motion (EOM) Kim and Kapral (2005a)

[TABLE]

where $\widetilde{\Lambda}_{2}=\overleftarrow{\nabla}_{\boldsymbol{\mathcal{P}}_{2}}\overrightarrow{\nabla}_{\boldsymbol{\mathcal{R}}_{2}}-\overleftarrow{\nabla}_{\boldsymbol{\mathcal{R}}_{2}}\overrightarrow{\nabla}_{\boldsymbol{\mathcal{P}}_{2}}$ is the Poisson bracket operator (with the direction of an arrow indicating the direction in which the operator acts).

The MGF in Eq. (15) is exact but computationally intractable in general because it involves a fully quantum mechanical treatment of the total system. By taking the quantum-classical limit of Eq. (15), one can obtain an expression that is amenable to numerical simulations. To take this limit, we first note that the full Wigner transform of an operator, $O_{W}(\boldsymbol{\mathcal{X}})$, may be written as

[TABLE]

where $\hat{O}_{W}(\boldsymbol{X})$ is the partially Wigner-transformed operator. For a quantity that depends only on the variables of the baths, one further has

[TABLE]

i.e., its full Wigner transform is equivalent to its partial Wigner transform. This is the case for the exponential functions $(e^{i\sum_{v}\chi_{v}H_{B}^{v}})_{W}$ and $(e^{-i\sum_{v}\chi_{v}H_{B}^{v}})_{W}$ in Eq. (15). Thus, the MGF in Eq. (15) reduces to

[TABLE]

where

[TABLE]

Similar to Eq. (18), the EOM for the weight function $\overline{D}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)$ reads

[TABLE]

where $\Lambda_{2}=\overleftarrow{\nabla}_{\boldsymbol{P}_{2}}\overrightarrow{\nabla}_{\boldsymbol{R}_{2}}-\overleftarrow{\nabla}_{\boldsymbol{R}_{2}}\overrightarrow{\nabla}_{\boldsymbol{P}_{2}}$ is the Poisson bracket operator that acts in the phase space of the heat baths. Up to this point, no approximations were employed, so Eq. (21) is equivalent to the exact quantum MGF in Eq. (12). However, the difficulties associated with solving Eq. (23) are formidable.

The quantum-classical limit of the MGF is then taken by replacing the evolution equation for $\overline{D}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)$ with its quantum-classical limit Sergi and Kapral (2004); Kim and Kapral (2005a, b), i.e., replacing the exponential operator in Eq. (23) with its expansion to first order in $\hbar$ Kapral and Ciccotti (1999):

[TABLE]

where the subscript “QC”denotes the quantum-classical limit, $\{\hat{H}_{W},\cdot\}_{a}=\frac{1}{2}\{\hat{H}_{W},\cdot\}-\frac{1}{2}\{\cdot,\hat{H}_{W}\}$ is the anti-symmetrized Poisson bracket, and $i\mathcal{L}$ is the quantum-classical Liouville operator. One can formally solve Eq. (24) to obtain

[TABLE]

In the above equation, we note that $\overline{D}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$, the zero-time limit of $\overline{D}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)$ in Eq. (III.1), is the initial condition for $\overline{D}_{QC}(t)$. As such, all of the quantum information is retained at the initial time. Thus, the basis-independent quantum-classical MGF is

[TABLE]

Noting that the quantum-classical Liouville operator in Eq. (25) depends only on $\boldsymbol{X}_{2}$, one can move the action of the evolution operator onto the term $\left(e^{-i\sum_{v}\chi_{v}H_{B}^{v}}\right)_{W}(\boldsymbol{X}_{2})$ Kim and Kapral (2005a) to obtain an equivalent expression for the MGF

[TABLE]

This equation serves as a convenient starting point for computations.

To evaluate the quantum-classical MGF in Eq. (27), one must insert complete sets of basis states $\{|\alpha\rangle\}$ that span the Hilbert space of the quantum subsystem

[TABLE]

where

[TABLE]

It should be noted that the MGF in Eq. (28) has a similar structure to that of a quantum correlation function in the quantum-classical limit Sergi and Kapral (2004); Kim and Kapral (2005a, b).

III.2 Average heat and heat current

In this subsection, we show how one can work out the expected expressions for the average heat and heat current from the quantum-classical MGF in Eq. (28). From the formal expression in Eq. (13), the average heat is given by

[TABLE]

The first term on the right-hand-side (RHS) of Eq. (30) may be simplified by first summing over $\alpha_{2}$ (using the completeness of the basis) and integrating over $\boldsymbol{P}_{2}$, which results in the delta function $\delta(\boldsymbol{Z}_{2})$. The integrals with respect to $\boldsymbol{Z}_{2}$ and $\boldsymbol{R}_{2}$ may then be evaluated to yield

[TABLE]

where we have used the fact that $\frac{1}{(2\pi\hbar)^{\mu_{B}}}\int d\boldsymbol{Z}_{1}e^{-\frac{i}{\hbar}\boldsymbol{P}_{1}\cdot\boldsymbol{Z}_{1}}\langle\boldsymbol{R}_{1}+\boldsymbol{Z}_{1}/2|\hat{\rho}_{0}|\boldsymbol{R}_{1}-\boldsymbol{Z}_{1}/2\rangle$ is the partially Wigner-transformed initial density matrix $\hat{\rho}_{W}(0)$. Similarly, for the second term on the RHS of Eq. (30), we have

[TABLE]

Summing these two terms together leads to the expression for the average heat in Eq. (4). By taking the time derivative of this expression, one obtains the expression in Eq. (5) for the heat current from the $v$th bath to the quantum subsystem. However, the time evolution in both expressions is now dictated by the quantum-classical Liouville operator.

IV Fluctuation symmetry in the quantum-classical limit

It has been previously shown that fully quantum composite systems exhibiting microscopic reversibility obey the following SSFS Nicolin and Segal (2011a) (in appendix A, we provide a proof of this from a closed system point of view)

[TABLE]

where $Z_{ss}(\{\chi_{v}\},t)\equiv\lim\limits_{t\to\infty}Z(\{\chi_{v}\},t)$. Physically speaking, the SSFS determines the heat fluctuations at steady state. If one now considers a two-heat bath (left and right) setup, the two counting fields $\chi_{L}$ and $\chi_{R}$ simply measure the same amount of energy in the steady state. Thus, if one introduces a new counting field $\chi=\chi_{R}-\chi_{L}$ (for the case in which the left bath has a higher temperature than the right bath), then the well-known heat exchange fluctuation symmetry is recovered from the SSFS above Esposito et al. (2009); Campisi et al. (2011)

[TABLE]

where $\Delta\beta=\beta_{R}-\beta_{L}$ is the thermodynamic affinity associated with the steady state heat current.

In the case of an arbitrary quantum subsystem bilinearly coupled to harmonic baths, the SSFS in Eq. (33) should exactly hold for the MGF $Z_{QC}(\vec{\chi},t)$ because the QCLE yields the exact quantum dynamics. On the other hand, for systems that go beyond the scope of bilinear interactions and harmonic environments, the approximations inherent to QCLE dynamics may alter the behavior of the heat fluctuations at steady state. It is therefore important to determine to what extent the SSFS holds in such systems. The MGF in Eq. (21) is equivalent to the exact quantum MGF and therefore its long-time limit satisfies the SSFS in Eq. (33). However, under QCLE dynamics, one can show that the long-time limit of the quantum-classical weight function in Eq. (26) has the following approximate form (see appendix B for details)

[TABLE]

where $\overline{D}_{ss}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)\equiv\lim\limits_{t\to\infty}\overline{D}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)$ and $\overline{D}_{QC,ss}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)\equiv\lim\limits_{t\to\infty}\overline{D}_{QC}(\boldsymbol{X}_{1},\boldsymbol{X}_{2},t)$, and the subscript in $\mathcal{O}_{t}(\hbar^{2})$ indicates that the correction term is time-dependent in the long-time limit. Given the above relation, it follows that the SSFS holds only up to order $\hbar$ under QCLE dynamics (see appendix C for details), i.e.,

[TABLE]

where $Z_{QC,ss}(\{\chi_{v}\},t)\equiv\lim\limits_{t\to\infty}Z_{QC}(\{\chi_{v}\},t)$ and $\mathcal{S}_{QC}(\{\chi_{v}\})\equiv\lim\limits_{t\to\infty}\frac{1}{t}\ln Z_{QC}(\{\chi_{v}\},t)$ is the quantum-classical scaled cumulant generating function of the heat current. It should be noted that $\mathcal{S}_{QC}(\{\chi_{v}\})$ is time-independent because we are considering situations where the cumulants of heat grow linearly with time. This is usually the case when measuring the statistics of quantities associated with nonequilibrium energy fluxes Esposito et al. (2009) (anomalous heat statistics have been observed in exceptional cases Esposito and Lindenberg (2008)).

To understand the physical implications of not strictly satisfying the SSFS, we focus on systems with two heat baths such that the scaled cumulant generating function in Eq. (36) reduces to

[TABLE]

By introducing the heat transport coefficients

[TABLE]

one can obtain the following Saito-Utsumi (SU) relations Saito and Utsumi (2008) in the quantum-classical limit

[TABLE]

If we now consider the case with $n=0$ and $m=2$, we find that

[TABLE]

where we have used the fact that $L_{2}^{0}(\Delta\beta)=0$ due to the normalization condition on the density matrix Saito and Utsumi (2008). In the above equation, the left-hand-side (LHS) equals $\frac{\partial}{\partial\Delta\beta}\langle J\rangle_{ss}$ with $\langle J\rangle_{ss}$ the steady state heat current from the hot bath to the cold one and the RHS is the variance of the steady state heat current. In the linear response regime, the LHS of Eq. (40) is proportional to the heat conductance. Thus, Eq. (40) reveals that the fluctuation-dissipation relation is satisfied up to order $\hbar$ in the quantum-classical limit.

We conclude this section by noting that, although the SSFS and fluctuation dissipation theorem are not strictly preserved under QCLE dynamics, the quantum-classical approximation becomes more accurate in the limit that the bath DOF are much heavier than the subsystem DOF and/or under high temperature conditions.

V Nonequilibrium spin-boson model

V.1 Model

To illustrate the formalism, we consider the NESB model, a prototypical model in the study of quantum energy transfer at the nanoscale Boudjada and Segal (2014). This model consists of an unbiased two-level subsystem in contact with two bosonic heat baths at different temperatures. The QCLE dynamics of the NESB model is dictated by the following Weyl-ordered, partially Wigner-transformed Hamiltonian

[TABLE]

where $\hat{\sigma}_{x/z}$ are the Pauli spin matrices, $\Delta$ is the tunneling frequency between the two states, and $C_{j,v}$ is the coupling coefficient between the spin and the $j$th harmonic oscillator in the $v$th heat bath. The bilinear coupling between the subsystem and $v$th heat bath is characterized by an Ohmic spectral density with an exponential cutoff, namely $I_{v}(\omega)=\frac{\xi_{v}}{2}\pi\omega e^{-\omega/\omega_{c,v}}$ with $\xi_{v}$ the Kondo parameter characterizing the subsystem-bath coupling strength and $\omega_{c,v}$ the cutoff frequency. In our simulations, we use dimensionless variables and parameters with time scaled by $\omega_{c}$.

The initial state of the system is chosen to be the product state $\hat{\rho}_{W}(0)=\hat{\rho}_{S}(0)\rho_{B,W}(0)$, where $\hat{\rho}_{S}(0)=|+\rangle\langle+|$ ($|+\rangle$ is the spin-up state of $\hat{\sigma}_{z}$) and $\rho_{B,W}(0)=\prod_{v}\rho_{B,W}^{v}(0)$ with

[TABLE]

the partially Wigner-transformed canonical distribution.

Sampling from the weight function $\overline{D}^{\alpha_{1}\alpha_{2}\alpha_{2}^{\prime}}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$ in Eq. (III.1) remains a challenging numerical task, so here we focus on simulations of the average heat and heat current. In this case, $\overline{D}^{\alpha_{1}\alpha_{2}\alpha_{2}^{\prime}}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$ reduces to the initial total density matrix (see Sec. III.2), which can be readily sampled from. The expression for the average heat transferred through the system (obtained from Eq. (4)) is

[TABLE]

where, for example, $\left(P_{j,v}^{2}(t)\right)^{\alpha\alpha^{\prime}}=\sum_{\beta}P_{j,v}^{\alpha\beta}(t)P_{j,v}^{\beta\alpha^{\prime}}(t)$. From this expression, we see that the time-dependent heat is determined by the time dependence of the matrix elements $(\boldsymbol{P}_{j,v}^{2})^{\alpha\alpha^{\prime}}(t)$ and $(\boldsymbol{R}_{j,v}^{2})^{\alpha\alpha^{\prime}}(t)$ originating from the bath Hamiltonian $H_{B,W}^{v}$. It should be noted that $P_{j,v}^{\alpha\alpha^{\prime}}(t)\neq P_{j,v}(t)\delta_{\alpha\alpha^{\prime}}$ because its time evolution depends on the subsystem’s operators due to the subsystem-bath coupling. The heat current $\langle J_{v}(t)\rangle$ is defined as the negative of the time derivative of $\langle Q_{v}(t)\rangle$:

[TABLE]

In our simulations, the heat current is obtained by simply calculating the derivative of $\langle Q_{v}(t)\rangle$ at each molecular dynamics (MD) time step.

V.2 Numerical simulations

V.2.1 DECIDE solution of QCLE

In order to evaluate the expressions for the average heat and heat current, we used a recently developed approximate solution of the QCLE known as the DECIDE (Deterministic Evolution of Coordinates with Initial Decoupled Equations) method Liu and Hanna (2018). Instead of propagating the observables directly as in the previous QCLE-based methods, DECIDE evolves the coordinates corresponding to the subsystem and bath (viz., $\boldsymbol{\hat{x}}(t)$ and $\boldsymbol{X}(t)$, respectively) according to the following set of equations of motion (EOMs)

[TABLE]

In the above equations, the time arguments are placed outside of their respective brackets to indicate that one should first evaluate the commutator and Poisson brackets with respect to the initial bath coordinates and then apply the time dependence to the coordinates in the resulting expressions.

As the DECIDE algorithm provides an approximate solution of the QCLE, it is worthwhile to discuss the core approximations that enter into the method. To arrive at Eq. (V.2.1), one starts with the partially Wigner-transformed (with respect to the initial bath coordinates) quantum Heisenberg equations for $\boldsymbol{\hat{x}}(t)$ and $\boldsymbol{\hat{X}}(t)$, and then truncates them by applying the following approximation for an arbitrary time-dependent operator $(\hat{B}(\boldsymbol{\hat{x}}(t),\boldsymbol{\hat{X}}(t)))_{W}\equiv(e^{i\hat{\mathcal{K}}t}\hat{B}(\boldsymbol{\hat{x}},\boldsymbol{\hat{X}}))_{W}$

[TABLE]

where $\mathcal{\hat{K}}$ is the quantum Liouville operator and $\Lambda$ is the Poisson bracket operator. To arrive at the second line of this equation, the quantum Liouville operator is replaced with the quantum-classical Liouville operator $i\mathcal{L}$ and only zeroth-order terms in $\hbar$ in the Moyal product expansion are retained. For the full details of the derivation of Eq. (V.2.1), we refer the readers to Ref. Liu and Hanna (2018).

The replacement of the quantum Liouville operator with the quantum-classical Liouville operator in Eq. (46) is exact if one considers harmonic environments and bilinear subsystem-bath interactions. However, by neglecting the higher order terms in $\hbar$ in the Moyal product, one may underestimate the back-action from the heat baths to the subsystem even in cases with harmonic environments and bilinear subsystem-bath interactions. To illustrate this, we focus on the first-order correction term to Eq. (46), namely $e^{i\mathcal{L}t}\left(\hbar\Lambda/2i\right)\hat{B}_{W}(\boldsymbol{\hat{x}},\boldsymbol{X})=-i\hbar\{e^{i\mathcal{L}t},\hat{B}_{W}(\boldsymbol{\hat{x}},\boldsymbol{X})\}$. For demonstration purposes, let us assume that $\hat{B}_{W}(\boldsymbol{\hat{x}},\boldsymbol{X})\propto\hat{x}\boldsymbol{R}$, which arises from a bilinear subsystem-bath interaction. In this case, the first-order correction term is

[TABLE]

To evaluate the derivative of $\hat{x}(t)$ with respect to the initial momenta, one must know the complete history of the dynamics from the initial time to time $t$. Thus, if the subsystem dynamics is highly non-Markovian, such correction terms cannot be ignored.

In light of its inherent approximations, the DECIDE solution can give rise to inaccurate results in the long-time limit in parameter regimes where non-Markovian effects are pronounced. We note that strong memory effects can be induced by strong subsystem-bath coupling, slow heat baths characterized by $\omega_{c}\ll\Delta$, and very low temperatures $k_{B}T_{v}\ll\hbar\Delta$. Therefore, DECIDE should be used with caution in such regimes. However, in regimes with weak memory effects, the contributions to the dynamics from the dropped terms in the EOMs for $\boldsymbol{\hat{x}}$ and $\boldsymbol{X}$ are negligible and DECIDE is expected to perform very well (as seen in Ref. Liu and Hanna (2018) and as will be shown below).

For the NESB model, the generalized coordinates of the spin subsystem are taken to be the Pauli matrices $(i.e.,\boldsymbol{\hat{x}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}))$. Before solving Eq. (V.2.1), one must cast the EOMs in an arbitrary basis $\{|\alpha\rangle\}$ that spans the $2\times 2$ Hilbert space of the two-level subsystem, namely

[TABLE]

where the dot denotes a time derivative. In total, there are $4\times(3+2N)$ (with $N=N_{L}+N_{R}$) coupled first-order differential equations (FODEs) for the matrix elements $(\sigma_{x}^{\{\alpha\alpha^{\prime}\}},\sigma_{y}^{\{\alpha\alpha^{\prime}\}},\sigma_{z}^{\{\alpha\alpha^{\prime}\}},\boldsymbol{X}^{\{\alpha\alpha^{\prime}\}})$, where $\{\alpha\alpha^{\prime}\}$ denotes all the combinations of basis indices. We remark that the superscript in $\boldsymbol{X}^{\alpha\alpha^{\prime}}(t)$ serves as a label to distinguish the various matrix elements that arise due to the subsystem-bath coupling.

V.2.2 Simulation details

To solve the FODEs in Eq. (V.2.1), we must first specify the nature of the $\{|\alpha\rangle\}$ basis set. In this work, we consider two basis sets that are frequently used in studies of the spin-boson model.

The first basis set is a subsystem basis, consisting of the eigenstates of $\hat{\sigma}_{z}$, i.e., $\{|\alpha\rangle\}=\{|+\rangle,|-\rangle\}$. In this basis, the initial values of the matrix elements of the subsystem coordinates are $\sigma_{x}^{+-}(0)=\sigma_{x}^{-+}(0)=1$, $\sigma_{x}^{++}(0)=\sigma_{x}^{--}(0)=0$, $\sigma_{y}^{++}(0)=\sigma_{y}^{--}(0)=0$, $\sigma_{y}^{-+}=i,\sigma_{y}^{+-}=-i$, $\sigma_{z}^{+-}(0)=\sigma_{z}^{-+}(0)=0$, $\sigma_{z}^{++}(0)=1,\sigma_{z}^{--}(0)=-1$; and the initial values of the matrix elements of the bath coordinates are $\boldsymbol{X}^{\alpha\alpha^{\prime}}(0)=\boldsymbol{X}(0)\delta_{\alpha\alpha^{\prime}}$ (due to the initial product state), with $\boldsymbol{X}(0)$ sampled from Eq. (42). In this basis, the expression for the average transferred heat in Eq. (43) reduces to

[TABLE]

using the fact that $\rho_{S}^{++}(0)=1$.

The second basis set is the adiabatic basis $\{|\alpha\rangle\}=\{|1\rangle,|2\rangle\}$, which can be expressed in terms of $|\pm\rangle$ as follows MacKernan et al. (2002)

[TABLE]

where $G=\frac{1}{\gamma(\boldsymbol{R})}[-\Delta+\sqrt{\Delta^{2}+\gamma(\boldsymbol{R})^{2}}]$ with $\gamma(\boldsymbol{R})=-\sum_{v}\sum_{j=1}^{N_{v}}C_{j,v}R_{j,v}$. In this basis, the initial conditions for the subsystem coordinates are $\sigma_{x}^{11}(0)=\frac{1-G^{2}}{1+G^{2}},\sigma_{x}^{12}(0)=\sigma_{x}^{21}(0)=\frac{2G}{1+G^{2}},\sigma_{x}^{22}(0)=-\frac{1-G^{2}}{1+G^{2}}$, $\sigma_{y}^{12}(0)=-i,\sigma_{y}^{21}(0)=i$, $\sigma_{y}^{11}=\sigma_{y}^{22}=0$, $\sigma_{z}^{11}(0)=\frac{2G}{1+G^{2}},\sigma_{z}^{12}(0)=\sigma_{z}^{21}(0)=-\frac{1-G^{2}}{1+G^{2}},\textrm{and}~{}\sigma_{z}^{22}(0)=-\frac{2G}{1+G^{2}}$ (see Appendix D for details), where $G$ is determined by the initial bath coordinates $\boldsymbol{R}(0)$; the initial values of the bath coordinates are again $\boldsymbol{X}^{\alpha\alpha^{\prime}}(0)=\boldsymbol{X}(0)\delta_{\alpha\alpha^{\prime}}$ (due to the initial product state), with $\boldsymbol{X}(0)$ sampled from Eq. (42). It should be noted that, because the adiabatic basis states $\{|1\rangle,|2\rangle\}$ are only used to set the initial values of the coordinates, one does not need to diagonalize the Hamiltonian matrix on-the-fly, in contrast to surface-hopping approaches. After setting the initial values of the coordinates, one just updates them by integrating Eq. (V.2.1). Using the adiabatic basis, the density operator corresponding to the initial spin-up state is given by

[TABLE]

Since the four matrix elements of $\hat{\rho}_{S}(0)$ are nonzero, there will be four non-zero components in the expression for the average transferred heat in Eq. (43).

To simulate the Ohmic spectral density with the exponential cutoff, we adopt a discretization scheme Thompson and Makri (1999); Wang et al. (2001) with

[TABLE]

where $j$ runs from 1 to $N_{v}$, $\omega_{0,v}=\frac{\omega_{c,v}}{N_{v}}(1-e^{-\omega_{m,v}/\omega_{c,v}})$, and $\omega_{m,v}$ is the maximum frequency of the $v$th heat bath. In our simulations, we take $\omega_{m,L}=\omega_{m,R}=\omega_{m}$ and $\omega_{c,L}=\omega_{c,R}=\omega_{c}$. Although we employ an Ohmic spectral density in this study, it should be emphasized that this approach, just like any other mixed quantum-classical dynamics method, can handle arbitrary bath spectral densities.

Finally, to integrate Eq. (V.2.1), we adopt the standard fourth-order Runge-Kutta scheme Dormand and Prince (1980). Noting that $(\boldsymbol{P}_{v}^{2})^{\alpha\alpha^{\prime}}(t)=\sum_{\beta}\boldsymbol{P}_{v}^{\alpha\beta}(t)\boldsymbol{P}_{v}^{\beta\alpha^{\prime}}(t)$ and $(\boldsymbol{R}_{v}^{2})^{\alpha\alpha^{\prime}}(t)=\sum_{\beta}\boldsymbol{R}_{v}^{\alpha\beta}(t)\boldsymbol{R}_{v}^{\beta\alpha^{\prime}}(t)$, the time evolution of the heat and heat current can then be constructed in terms of the time-dependent coordinates by averaging over an ensemble of trajectories according to Eqs. (43) and (44).

V.2.3 Equilibrium condition

Before considering a temperature gap between the two baths, it is instructive to first consider the equilibrium case ($\beta_{L}=\beta_{R}$) and investigate whether the heat currents of the left and right bath vanish at steady state. In doing so, we also consider symmetric subsystem-bath couplings (i.e., $\xi_{L}=\xi_{R}$) and asymmetric ones (i.e., $\xi_{L}\neq\xi_{R}$), because numerical methods may predict vanishing heat currents in the long-time limit in one case and fail in the other.

In the symmetric coupling case shown in Fig. 1, we see that both $\langle Q_{L}(t)\rangle$ and $\langle Q_{R}(t)\rangle$ become constant (to within numerical error) in the long-time limit, as expected (see Figs. 1 (a) and (b) for the results obtained using the adiabatic and subsystem bases, respectively). In principle, $\langle Q_{L}(t)\rangle$ and $\langle Q_{R}(t)\rangle$ should be identical in the symmetric coupling case, but minor deviations between them are observed due to the fact that Eq. (V.2.1) is not exact.

Nevertheless, the resulting $\langle J_{L}(t)\rangle$ and $\langle J_{R}(t)\rangle$ vanish in the long-time limit, implying that the energy conservation condition is satisfied by our simulations. In the transient regime, we find that the left and right bath heat currents are negative, which (according to our sign convention) means that heat is flowing into the baths. This behaviour has been previously observed and is due to the sudden switch-on of the subsystem-bath couplings at $t=0$ Cuansing and Wang (2010).

On the other hand, in the asymmetric coupling case shown in Fig. 2, we see that $\langle Q_{R}(t)\rangle$ (see Figs. 2 (a) and (b) for the results obtained using the adiabatic and subsystem bases, respectively) has a larger absolute steady state value due to a larger coupling strength between the subsystem and right heat bath. The ratio between the two steady state heat values is about 2, which is consistent with the ratio of the coupling strengths.

The time dependences of $\langle J_{L}(t)\rangle$ and $\langle J_{R}(t)\rangle$ at short times are also different, the latter having a larger drop. However, in the long-time limit, both currents still vanish identically. Furthermore, by comparing the results from the two bases in Figs. 1 and 2, we find that the results are indeed basis-independent in both the symmetric and asymmetric coupling cases, pointing to the utility of the DECIDE method for simulating heat transfer processes.

V.2.4 Nonequilibrium condition

We now consider the nonequilibrium case where the temperatures of the two baths are not equal ($\beta_{L}\neq\beta_{R}$). In light of the results of the previous subsection, in this case, we only use the subsystem basis to carry out our calculations and only consider symmetric subsystem-bath couplings ($\xi_{L}=\xi_{R}=\xi$).

The results for the time-dependent average heat and heat current under different subsystem-bath coupling and temperature conditions are shown in Figs. 3, 4, and 5.

We first focus on the results obtained with high bath temperatures in Fig. 3. From Figs. 3 (a) and (b), we see that, at very short times, both $\langle Q_{L}(t)\rangle$ and $\langle Q_{R}(t)\rangle$ are the same because the total system starts from an initial product state and it takes time for the system to adjust to the temperature difference. However, at longer times, $\langle Q_{L}(t)\rangle$ and $\langle Q_{R}(t)\rangle$ begin to exhibit differences and ultimately grow linearly with time with opposite slopes, resulting in stationary heat currents. For $\langle Q_{L}(t)\rangle$, the slope is positive, so heat is leaving the left (higher temperature) heat bath, while for $\langle Q_{R}(t)\rangle$, the slope is negative, so heat is entering the right (lower temperature) heat bath. This behavior was also observed in open quantum linear systems (where a quantum harmonic oscillator is coupled to two harmonic heat baths) by using the nonequilibrium Green’s function method Agarwalla et al. (2012). As for the left and right bath heat currents (see Figs. 3 (c) and (d)), their short-time behaviors are similar to those in the equilibrium case. At later times, the currents ultimately plateau with positive and negative values for $\langle J_{L}(t)\rangle$ and $\langle J_{R}(t)\rangle$, respectively, i.e., heat is flowing from the hot to the cold bath as expected. In comparing the left and right panels, we see that, in the strong subsystem-bath coupling case, the oscillations in the heat current are more pronounced (as the back-action from the heat baths becomes stronger) and that the total system takes a longer time to evolve to its steady state. As a benchmark, in Fig. 3 (c), we provide the value of the steady state heat current predicted by the quantum master equation (QME) Segal and Nitzan (2005) for the weak subsystem-bath coupling case. As can be seen, there is a very small deviation between the DECIDE and QME results, which is not surprising as the former is an approximate method and the latter becomes exact in the weak coupling regime at high temperatures. We further note that the behaviors of the heat currents are qualitatively similar to the regularized heat currents obtained by the multilayer multiconfiguration time-dependent Hartree approach in Ref. Velizhanin et al. (2008).

If we further increase the subsystem-bath coupling strength (from $\xi=0.5$ to $\xi=2$),

we observe large fluctuations in the heat currents, even though the heat curves are quite smooth (see inset of Fig. 4). Due to the large magnitude of the heat current fluctuations, one cannot therefore extract a meaningful steady-state heat current (since the actual steady-state heat current is of order $10^{-3}$ according to the non-interacting blip approximation Nicolin and Segal (2011a)). This is expected because the approximation of truncating the exact EOMs for the subsystem and bath coordinates in the DECIDE method deteriorates in the strong coupling regime. Thus, one should be cautious when applying DECIDE to strong subsystem-bath coupling cases.

For heat baths at low temperatures and relatively small subsystem-bath coupling strengths (see Fig. 5), we find that the transient behaviors of both the transferred heat and heat current are similar to those at high temperatures (see Fig. 3). However, we see that the heat current curve in the weaker coupling regime (see Fig. 5 (c)) is smoother than that in the high temperature case (see Fig. 3 (c)) due to the suppression of the thermal noise from the heat baths at lower temperatures. We also notice that the recurrence of the heat current after its initial drop takes longer than in the high temperature case because of the smaller thermodynamic force (resulting from a smaller temperature difference) at the lower temperature.

VI Summary

In this work, we presented a general formalism for studying nonequilibrium heat transfer processes in mixed quantum-classical systems that combines the FCS and QCLE approaches. In particular, starting from its exact definition from FCS, we derived a general expression for the MGF of heat in the partial Wigner representation whose dynamics is prescribed by the QCLE. Using this expression, we obtained explicit expressions for the time-dependent average heat and heat current in a system. Owing to its mixed quantum-classical nature, this formalism offers a computationally efficient way of studying quantum heat transfer in realistic molecular environments at the nanoscale.

Since approximations that lead to mixed quantum-classical treatments are expected to alter the behavior of the heat fluctuations at steady state, we further considered to what extent the SSFS holds under QCLE dynamics. We found that the SSFS is preserved up to order $\hbar$ for systems that are beyond the scope of bilinear subsystem-bath interactions and harmonic baths. Using the SU relations, we also showed that a violation of the SSFS is related to a breakdown of the fluctuation-dissipation theorem in linear response regimes. However, if one considers systems in which the bath DOF are much heavier than the subsystem DOF and/or are at high temperatures, the approximations inherent to QCLE dynamics are not expected to significantly affect the SSFS and fluctuation-dissipation theorem.

We demonstrated the performance of this formalism by computing the time-dependent average heat and heat current for the NESB model using the recently developed DECIDE solution of the QCLE. Under equilibrium conditions (i.e., heat baths at the same temperature), DECIDE yields the expected zero steady state heat currents for both symmetric and asymmetric subsystem-bath couplings. Under nonequilibrium conditions (i.e., heat baths at different temperatures), DECIDE also yields the expected trends in the average heat and heat current, as compared to previous results obtained with fully quantum methods. Therefore, the present formalism together with the DECIDE method provide a valuable approach for simulating energy transfer processes in open quantum systems out of equilibrium.

Future studies will aim at analyzing the steady state heat current in the NESB model over a wide parameter space using the present method. In particular, it is essential to demonstrate whether the present method can reproduce the turn-over behaviour in the steady state heat current as a function of subsystem-bath coupling strength Velizhanin et al. (2008). The generalization of the method to calculate higher order heat fluctuations is also worthwhile. For instance, the noise power of the heat current, obtained from the second order cumulant of the heat, provides rich information beyond what could be inferred from the average heat and heat current Liu et al. (2018). We also anticipate applications of the method to multi-level (i.e., beyond two levels) subsystems with more complex environments. Finally, one could consider applying other mixed quantum-classical and semi-classical methods, such as those previously used in the study of vibrational energy transfer in condensed phases Shi and Geva (2003); Hanna and Geva (2008); Jain and Subotnik (2018), to nonequilibrium heat transfer problems.

Acknowledgements.

J. Liu and G. Hanna acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). C.-Y. Hsieh acknowledges support from the Singapore-MIT Alliance for Research and Technology (SMART). D. Segal acknowledges support from an NSERC Discovery Grant and the Canada Research Chair program.

Appendix A Fluctuation symmetry in the long-time limit

In this appendix, we provide an alternative proof to that in Ref. Nicolin and Segal (2011a) that, in quantum composite systems with time-reversal symmetry, the MGF of heat satisfies the symmetry relation $Z_{ss}(\{\chi_{v}\},t)=Z_{ss}(\{i\beta_{v}-\chi_{v}\},t)$ in the long-time limit, where $\beta_{v}$ is the inverse temperature of $v$th heat bath.

To start, by using the facts that the trace in Eq. (12) is invariant to cyclic permutation and that $[\hat{\rho}_{0},\hat{H}_{B}^{v}]=0$, we may re-express the MGF as

[TABLE]

In the absence of any external driving, we can shift the time arguments in Eq. (53) as follows by noting that $\hat{U}(t,0)=\hat{U}(t/2,0)\hat{U}(t/2,0)$:

[TABLE]

with $\hat{\rho}\left(t/2\right)\equiv\hat{U}(t/2,0)\hat{\rho}_{0}\hat{U}^{\dagger}(t/2,0)$. Next, applying the transformation $\chi_{v}\to i\beta_{v}-\chi_{v}$ to Eq. (54), we obtain

[TABLE]

Since $\hat{\rho}_{0}$ commutes with the bath Hamiltonians, we can show that

[TABLE]

In the limit of $t\to\infty$, the RHS of Eq. (A) becomes

[TABLE]

Given the micro-reversibility of the closed system without external driving, one may argue that $\hat{H}_{B}^{v}(\infty)$ and $\hat{H}_{B}^{v}(-\infty)$ are equal and, therefore, Eq. (57) reduces to $\hat{\rho}_{ss}$ in the long-time limit.

To complete the proof, we will need the following relations for an arbitrary operator $\hat{B}$:

[TABLE]

where $\Theta$ is the quantum mechanical time-reversal operator and $\hat{\rho}_{ss}$ is the steady state density operator of the total system. In arriving to the last equality in the above equation, we used the fact that the time evolution operator commutes with $\hat{\rho}_{ss}$. Now, similarly to what is done above, we can shift the time arguments in Eq. (12), take the long-time limit, and use Eq. (58) to yield the following:

[TABLE]

Finally, comparing Eq. (59) with the long-time limit of Eq. (55) (obtained with the aid of Eq. (57)), we see that

[TABLE]

Appendix B Long-time limit of $\overline{D}_{QC}$

Given Eqs. (23) and (24), $\overline{D}$ and $\overline{D}_{QC}$ in the long-time limit should respectively satisfy the following equations:

[TABLE]

In order to analyze the connection between $\overline{D}_{ss}(t)$ and its approximated form $\overline{D}_{QC,ss}(t)$, we first expand these quantities in power series of $\hbar$ Nielsen et al. (2001):

[TABLE]

We then substitute these power series back into Eq. (B) and group by powers of $\hbar$. For $\overline{D}_{ss}(t)$, this leads to

[TABLE]

and so on. For $\overline{D}_{QC,ss}$, this leads to the following recursion relations:

For $\hbar^{0}$ order,

[TABLE]

and for $\hbar^{n}$ order with $n\geqslant 1$,

[TABLE]

Comparing Eqs. (65) and (66) with Eq. (B), we find that $\overline{D}_{QC,ss}(t)$ and $\overline{D}_{ss}(t)$ are identical to order $\hbar$, namely, $\overline{D}_{QC,ss}^{(0)}(t)=\overline{D}_{ss}^{(0)}(t)$, $\overline{D}_{QC,ss}^{(1)}(t)=\overline{D}_{ss}^{(1)}(t)$, and $\overline{D}_{QC,ss}^{(n)}(t)\neq\overline{D}_{ss}^{(n)}(t)$ for $n>1$, thereby proving Eq. (35) in the main text.

Appendix C Fluctuation symmetry in the quantum-classical limit

The long-time limits of the quantum and quantum-classical MGFs in Eqs. (21) and (26), respectively, are

[TABLE]

where $Z_{QC,ss}(\{\chi_{v}\},t)\equiv\lim\limits_{t\to\infty}Z_{QC}(\{\chi_{v}\},t)$. We can expand these MGFs in power series of $\hbar$ (as done in Eqs. (62) and (63)) to yield

[TABLE]

with $Z_{ss}^{(n)}(\{\chi_{v}\},t)$ and $Z_{QC,ss}^{(n)}(\{\chi_{v}\},t)$ solely determined by $\overline{D}_{ss}^{(n)}(t)$ and $\overline{D}_{QC,ss}^{(n)}(t)$, respectively, e.g., $Z_{ss}^{(2)}(\{\chi_{v}\},t)$ does not depend on $\overline{D}_{ss}^{(1)}(t)$ because they are associated with different orders of $\hbar$. Given the analysis in appendix B for the weight functions, one can conclude that $Z_{QC,ss}(\{\chi_{v}\},t)$ and $Z_{ss}(\{\chi_{v}\},t)$ are also identical to order of $\hbar$, i.e.,

[TABLE]

where $\mathcal{O}_{t}(\hbar^{2})$ is time-dependent. Since $Z_{ss}(\{\chi_{v}\},t)$ preserves the SSFS, it immediately follows that

[TABLE]

Based on Eq. (13), we can further rewrite the MGFs in terms of the moments of heat as

[TABLE]

where the other counting fields except for $\chi_{v}$ are zero. If we now expand the moments of heat in power series of $\hbar$, and compare Eqs. (69) and (70) with Eqs. (73) and (74), respectively, we obtain the following relations

[TABLE]

Therefore, according to Eq. (71), we find that the moments of heat $\left\langle Q_{v}^{m}(t)\right\rangle_{ss}$ and $\left\langle Q_{v}^{m}(t)\right\rangle_{QC,ss}$ are identical to order $\hbar$, i.e.,

[TABLE]

Next, we introduce the long-time limits of the cumulant generating functions (CGFs) of heat $G_{ss}(\{\chi_{v}\},t)\equiv\ln Z_{ss}(\{\chi_{v}\},t)$ and $G_{QC,ss}(\{\chi_{v}\},t)\equiv\ln Z_{QC,ss}(\{\chi_{v}\},t)$. Given these definitions, one can establish the following relationship between the quantum and quantum-classical cumulants of heat $\left\langle\langle Q_{v}^{m}(t)\right\rangle\rangle_{ss}$ and $\left\langle\langle Q_{v}^{m}(t)\right\rangle\rangle_{QC,ss}$, respectively:

[TABLE]

In most cases, the cumulants of heat grow linearly with time Esposito et al. (2009), so one may argue that $\mathcal{O}_{t}(\hbar^{2})\sim t\cdot\mathcal{O}(\hbar^{2})$ in the above equation and consequently in Eq. (72), thereby proving the first equation in Eq. (36). Given this linear time dependence of the cumulants in the long-time limit, one may define time-independent scaled CGFs of the heat current as follows

[TABLE]

Since the CGFs can also be expanded in terms of the cumulants of heat (in analogy with Eqs. (73) and (74) for the MGFs), Eq. (78) implies that

[TABLE]

Finally, given the fact that $\mathcal{S}(\{\chi_{v}\})$ preserves the SSFS, one can recover the second equation in Eq. (36).

Appendix D Pauli matrices in the adiabatic basis

Using Eq. (V.2.2), one can express the Pauli matrices in the adiabatic basis as

[TABLE]

[TABLE]

From these expressions, one can determine the initial values of the subsystem coordinates given below Eq. (V.2.2).