Thermal power of heat flow through a qubit

TL;DR

This paper investigates the thermal power of heat flow through a qubit coupled to two thermal baths, deriving expressions for power in different interaction regimes using the spin-boson model.

Contribution

It introduces a method to compute thermal power in a qubit-bath system using the NIBA approximation and the polaron picture, connecting correlation functions to heat flow.

Findings

Derived expressions for thermal power in the weak interaction limit.

Connected thermal power to bath correlation functions and energy splitting.

Recovered known results in the weak coupling regime.

Abstract

In this paper we consider thermal power of a heat flow through a qubit between two baths. The baths are modeled as set of harmonic oscillators initially at equilibrium, at two temperatures. Heat is defined as the change of energy of the cold bath, and thermal power is defined as expected heat per unit time, in the long-time limit. The qubit and the baths interact as in the spin-boson model, i.e. through qubit operator $\sigma_z$. We compute thermal power in an approximation analogous to `non-interacting blip' (NIBA) and express it in the polaron picture as products of correlation functions of the two baths, and a time derivative of a correlation function of the cold bath. In the limit of weak interaction we recover known results in terms of a sum of correlation functions of the two baths, a correlation functions of the cold bath only, and the energy split.

Equations

$$H_{TOT}=H_S+H_C+H_H+H_{CS}+H_{HS}$$

$$H_S=-\hslash \frac{\Delta }{2}{\hat{\sigma }}_x+\frac{ϵ}{2}{\hat{\sigma }}_z$$

$$H_C$$

$$H_H$$

$$H_{CS}$$

$$H_{HS}$$

$$\Pi =P\left(up\right){\pi }_{up}+P\left(down\right){\pi }_{down}$$

$$P\left(up\right)$$

$$P\left(down\right)$$

$${\pi }_{up}$$

$${\pi }_{down}$$

$${\Pi }^{weak}$$

$${\Pi }^{strong}$$

$${\overline{p}}_C(\Delta E,f|i)=\sum _{{\mathcal{E}}_i^{(C)}}{p}_C(\Delta E,f|{\mathcal{E}}_i^{(C)},i)\frac{e^{-\beta E({\mathcal{E}}_i^{(C)})}}{Z_C(\beta )}$$

$${\overline{p}}_C(\Delta E,f|i)$$

$${\overline{p}}_L(\Delta E,f|i)=\frac{1}{2\pi }\int e^{-i\nu \Delta E}{G}_{if}(\nu )d\nu$$

$$G_{if}(\nu )$$

$$⟨\Delta E_C⟩$$

$$A=\frac{{\Delta }^2}{2}\int e^{-\frac{1}{\hslash }(S_C+S_H)}\cos \frac{1}{\hslash }(X_C+X_H-ϵ\Delta t)d\Delta t$$

$$D=\frac{{\Delta }^2}{2}\int e^{-\frac{1}{\hslash }(S_C+S_H)}\cos \frac{1}{\hslash }(X_C+X_H+ϵ\Delta t)d\Delta t$$

$$C_{+}^C(t)$$

$$C_{+}^H(t)$$

$$C_{-}^C(t)$$

$$C_{-}^H(t)$$

$$A=\frac{{\Delta }^2}{4}\int \left(C_{+}^C(t)C_{+}^H(t)e^{-i\frac{ϵt}{\hslash }}+C_{-}^C(t)C_{-}^H(t)e^{i\frac{ϵt}{\hslash }}\right)dt$$

$${\pi }_{up}$$

$${\pi }_{down}$$

$$A$$

$${⟨{\hat{B}}_{-}(t){\hat{B}}_{+}(0)⟩}_{eq}\approx 1+\frac{4}{{\hslash }^2}\sum _b\frac{C_b^2}{m_b^2{\omega }_b^4}{⟨{\hat{p}}_b(0){\hat{p}}_b(t)⟩}_{eq}$$

$${\pi }_{up}$$

$${\pi }_{down}$$

$$H_S=-\hslash \frac{\Delta }{2}{\hat{\sigma }}_x+\frac{ϵ}{2}{\hat{\sigma }}_z$$

$$H_C$$

$$H_H$$

Full text

Stochastic thermodynamics, quantum power operators, quantum heat switches

I Introduction

Heat and work in classical thermodynamics are properties of processes, and not states. Heat is further in classical thermodynamics energy transferred from the system to an uncontrolled environment such that it cannot later be retrieved to do useful work Sekimoto (2010); Seifert (2012). The translation of these concepts to the quantum domain is therefore not obvious, as discussed in an early review Esposito et al. (2009). Quantum thermal power is average quantum heat per unit time, and is a centrally important topic for e.g. applications to quantum heat engines Kosloff and Levy (2014); Pekola (2015); Vinjanampathy and Anders (2016). While heat and thermal power at weak coupling has been studied for a long time in the literature Weiss (2012); Alicki and Lendi (1987); Karimi and Pekola (2016), the attention to systems interacting strongly with one or several baths is more recent, see cf Esposito et al. (2015); Newman et al. (2017); Goyal and Kawai (2017); Ronzani et al. (2018); Dou et al. (2018); Kwon et al. (2018); Perarnau-Llobet et al. (2018).

The goal of this paper is to revisit these questions in perhaps the simplest non-trivial scenario: one qubit interacting strongly with two heat baths at different temperatures. We will start from the general and unifying point of view that heat is energy change in a bath. Thermal power is thus expected energy change in a bath per unit time, in the long-term limit. For concreteness we will consider thermal power as heat per time to the cold bath, and thus a quantity that has to be non-negative in the long term limit. We assume that the qubit interacts with the baths and with an external drive as in the spin-boson model which allows to re-use many results developed in that literature Leggett et al. (1987). At strong coupling, and in the approximation known as “non-interacting blip approximation” (NIBA), the stationary state of the qubit is then determined by equilibrium correlation functions of certain bath operators related to a polaron transform. Our main result is that in a similar approximation thermal power is determined by derivatives of the same correlation functions with respect to time.

The paper is organized as follows: in Section II we introduce our model, and in Section III we give dimensional arguments what the results should be, first in a version appropriate for weak coupling, and then in a version appropriate for strong coupling. Section IV contains an overview of the calculations, and states the results in path integral language while Section V states in the language of the correlation functions after the polaron transform. Section VI summarizes and discusses the results.

Some of the background and much of the calculations are presented in appendices. Appendix A is thus a summary of the vast literature on the spin-boson problem, sufficient for our purposes. Appendix B summarizes on the other hand earlier work on quantum heat functionals Aurell and Eichhorn (2015); Aurell (2017, 2018a) adapted to the spin-boson setting, and Appendix C gives details of what these formulas mean for Ohmic baths. Appendix D further translates this theory to when the interaction is through bath momentum. Appendices E-J finally contain details of the calculations presented in Section IV.

II The model

We consider one qubit interacting with two harmonic oscillator baths as in the spin-boson model Leggett et al. (1987). Harmonic oscillator baths model, for instance, resistive elements in electrical circuits, and quantum harmonic oscillator baths hence model how such elements interact with other circuit elements at sufficiently low temperatures Devoret (1995). Circuits with superconducting elements that can be assimilated to qubits are widely investigated in scalable quantum information processing Wendin (2017). The state of one qubit interacting with two baths is hence a toy model of a quantum computer perturbed by a heat flow through the dynamical degrees of freedom of quantum computer itself. Quantum thermal power in this setting is conversely how well such a device can transport energy between two baths in the quantum regime.

The system, the baths and the interactions can thus be written down as a total Hamiltonian

[TABLE]

where “$C$” refers to the cold bath (temperature $T_{C}$) and “$H$” refers to the hot bath (temperature $T_{H}$).

The system Hamiltonian is

[TABLE]

where $\Delta$ is a rate (dimension $\left(\hbox{time}\right)^{-1}$), and $\epsilon$ is the level splitting. The bath Hamiltonian are

[TABLE]

where the parameters $m_{b}$ and $\omega_{b}$ are the mass and angular frequency of each oscillator and $C$ and $H$ also stand for the sets of oscillators in respectively the cold bath and the hot bath.

We will take the system-bath interactions to be described by

[TABLE]

where $C_{b}$ is the interaction coefficient between bath oscillator $b$ and the qubit, $q_{b}$ is the oscillator coordinate, and $\hat{\sigma}_{z}$ operates on the qubit. Pauli matrices are by convention dimension-less, and the coupling coefficients $C_{b}$ hence have dimension $(\hbox{energy})\cdot(\hbox{length})^{-1}$. In Leggett et al. (1987) the length scale (there called $q_{0}$) is taken to be the spatial distance between the minima of two potential wells. For a qubit formed out of a non-linear oscillator the length scale could similarly be the typical spatial scale of the oscillator ground state, $\sqrt{\frac{\hbar}{m\omega}}$.

We consider heat as related to two measurements on the cold bath, one at the beginning of the process and one at the end, which we assume to take values ${\cal E}_{i}$ and ${\cal E}_{f}$. In a quantum bath neither ${\cal E}_{i}$ nor ${\cal E}_{f}$ are known; all we can know is the probability of observing ${\cal E}_{i}$ at the beginning, and the probability of observing ${\cal E}_{f}$ at the end. Thermal power is then the expected change of bath energy per unit time $\frac{\left<{\cal E}_{f}-{\cal E}_{i}\right>}{t_{f}-t_{i}}$.

Four remarks are in order. First, “measurement on the bath” is required in the theory we consider, as without measurement the bath energy does not have a definite value. However, expected heat per unit time can, as we will see, be expressed in terms of system properties alone. Thermal power hence does not make any direct references to measurement, the values of which can hence be taken to be unrecorded. We may thus imagine “measurement on the bath” to actually refer to interaction with a large super-bath which forces the bath states to decohere, without assuming any direct control of the bath states by an experimenter. Second, we do not count any part of the interaction energy in the heat. While this issue is important and has been discussed at length on the classical side in the recent literature Seifert (2016); Talkner and Hänggi (2016); Jarzynski (2017); Miller and Anders (2017); Aurell (2017, 2018b), it is reasonable to assume that the interaction energy between one qubit and a bath does not increase at a non-zero rate for long enough times. Third, in applications to superconducting circuits, the system-bath interaction may often more naturally be taken to be proportional to bath oscillator momentum variable $p_{b}$ Devoret (1995). Since both $q_{b}$ and $p_{b}$ can be expressed in Fourier modes of the oscillator this can be expected to make no essential difference, as was indeed stated in Caldeira and Leggett (1983) for the qubit state. For completeness we outline in Appendix D an argument that this is so also for heat (full distribution function of bath energy change). Lastly, in realistic mesoscopic devices effective temperatures of different parts may differ. Such situations fall outside what is considered here, since the devices would then not be systems in thermal equilibrium that could be modelled as baths.

III Dimensional arguments

The long-time limit of the state of one qubit interacting with any number of baths is given by its density matrix, where the diagonal terms (“the populations”) determine the probability for the qubit to be respectively in the up state and in the down state. Suppose these probabilities are $P\left(\hbox{up}\right)$ and $P\left(\hbox{down}\right)$. Suppose further that the memory of the bath is short enough that when the system is in one state the bath does not remember in which states the system was before. We can then suppose that the expected energy given to the cold bath per unit time takes two values that depend on the system state, call them $\pi_{\hbox{up}}$ and $\pi_{\hbox{down}}$. Thermal power can then be estimated as

[TABLE]

To turn this into a quantitative prediction we can suppose that qubit transitions happen with effective rates describing the interactions with the two baths, and call these rates $\Gamma^{C}_{\downarrow\uparrow}$, $\Gamma^{H}_{\downarrow\uparrow}$, $\Gamma^{C}_{\uparrow\downarrow}$ and $\Gamma^{H}_{\uparrow\downarrow}$. This approach is appropriate when the qubit is weakly coupled to the baths, and one considers sufficiently long time scales Weiss (2012); Alicki and Lendi (1987). The up and down probabilities then depend on the rates as for a classical jump process i.e. as

[TABLE]

Power is dimensionally energy per unit time. When interaction energy is negligible the characteristic scale of energy transferred to the cold bath must be $\epsilon$ in an up-to-down transition, and $-\epsilon$ in an down-to-up transition, and these happen with rates $\Gamma^{C}_{\uparrow\downarrow}$ and $\Gamma^{C}_{\downarrow\uparrow}$. This leads to the estimates of power in the two states as

[TABLE]

and overall expected power as

[TABLE]

Expressions of this form are well known in the literature, e.g. in Karimi and Pekola (2016) (Eq. 5), and essentially hold in weak coupling also without the assumption of a short bath memory time.

At strong coupling the above is however not correct because when the qubit flips there is also a change of interaction energy between qubit and the bath. When this is larger than the level splitting the characteristic scale of energy transferred to the bath can be very different from $\epsilon$. Furthermore, in strong coupling one may assume combined effective mean switching rates $\Gamma_{\uparrow\downarrow}$ and $\Gamma_{\downarrow\uparrow}$, but it is not possible to disentangle the actions of the two baths into separate terms $\Gamma^{C}$ and $\Gamma^{H}$.

A different argument can nevertheless be made using the assumption of short enough bath de-correlation time, or equivalently that $\Delta$ is small enough that the residence time of the qubit in one state is long enough. From one qubit jump to the next qubit jump the baths hence on the average behave as follows. Right after the jump into state $s$ there will be some average interaction energy and some average bath energy, $\left<H_{CS}^{i}(s)\right>$ and $\left<H_{C}^{i}(s)\right>$. Between the jumps, when the qubit does not change its state, the sum of these energies is conserved, but in the same time interval the baths will equilibriate with the qubit. At the end of the interval the average interaction energy should hence vanish. This means that during one residence time in state $s$ the expected energy change of the bath should be the expected initial interaction energy i.e. $\left<H_{CS}^{i}(s)\right>$. By this reasoning one gets

[TABLE]

The main contribution of this paper is to derive an estimate like (11) systematically, and explain how the terms follow from the microscopic parameters of the model.

IV Thermal power at strong coupling

We now describe an approach to thermal power at strong coupling based on the Feynman-Vernon formalism Feynman and Vernon (1963). To calculate heat (energy change in a bath) we follow Aurell and Eichhorn (2015); Aurell (2017, 2018b), related general results can also been found in Carrega et al. (2015, 2016) and Funo and Quan (2018). Adapting the Feynman-Vernon formalism to describe the development of one spin interacting with one bath (the spin-boson problem) is already not trivial Leggett et al. (1987). Here we have the complications that we are interested in heat in a spin interacting with two (or more) baths at different temperatures. Technical background and details have therefore been moved to appendices; here we only outline the main idea of the calculation.

We focus on the energy changes of one bath, for concreteness we assume that is the cold bath. The starting point is to assume that initially the baths are independently at thermal equilibrium (at different temperatures), and the system as well as the energy of the cold bath are measured. After that measurement the state of the system and the baths is $\rho_{H}^{eq}\oplus\outerproduct{{\cal E}_{i}^{(C)},i}{{\cal E}_{i}^{(C)},i}$ where $\rho_{H}^{eq}$ is the equilibrium state of the hot bath (or baths).. $i$ indicates the state of the system after measurement and ${\cal E}_{i}^{(C)}$ the state of the cold bath. We take $p_{C}(\Delta E,f|{\cal E}_{i}^{(C)},i)$ to be the conditional probability of observing a final state $\ket{f}$ of the system and energy change of the cold bath $\Delta E$, conditioned on total initial state.

Next we assume that the measured energy of the cold bath is not recorded. This means that we could also say that the cold bath de-coheres by interacting with an unobserved cold super-bath at the same temperature. The initial state of the cold bath is then a statistical mixture where $\ket{{\cal E}_{i}^{(C)}}$ appears with the Gibbs weight $Z_{C}^{-1}(\beta)\exp\left(-\beta E({\cal E}_{i}^{(C)})\right)$. Here $\beta$ is the inverse temperature of the cold bath, and $Z_{C}$ is the partition function. From here we consider the average distribution

[TABLE]

which can be re-written

[TABLE]

where $\rho_{TOT}({\cal E}_{i},i)$ is the total density operator of the system and the bath at the end of the process, when the system and the cold bath started in the pure state $\ket{{\cal E}_{i},i}$. Resolving the delta function one can write

[TABLE]

where

[TABLE]

By linearity the Gibbs weight and the factor $e^{-i\nu E({\cal E}_{i})}$ can be taken inside the the big unitary transformation defining $\rho_{TOT}({\cal E}_{i},i)$. The above is therefore the same as

[TABLE]

where $\rho_{i}^{TOT}=\rho_{H}^{eq}\oplus\outerproduct{i}{i}\oplus\rho_{C}^{eq}$, and the trace is over the cold and the hot bath(s).

$G_{if}(\nu)$ codifies all the information on the distribution of energy change in a bath (here the cold bath), averaged over an initial equilibrium distribution of the baths at their respective temperatures and conditioned on the system starting in pure state $\ket{i}$ and finishing in pure state $\ket{f}$. Derivatives of $G_{if}(\nu)$ with respect to $\nu$ generate moments of the energy change. Here we are interested in the first derivative

[TABLE]

Furthermore we are only interested in thermal power, the limit $\frac{1}{t}\left<\Delta E_{C}\right>$ when $t$, the duration of the process, is long.

Stepping first back a bit, the calculation of $G_{if}(\nu)$ proceeds by representing $U$ and $U^{\dagger}$ as path integrals. Path integrals for spins are known in general Atland and Simons (2006), and have recently been used by one of us to estimate the errors in quantum computing Aurell (2018c). For the problem at hand a much simpler representation is however sufficient, where the spin paths $X$ and $Y$ representing $U$ and $U^{\dagger}$ are piece-wise constant, taking values $\pm\frac{1}{2}$ Leggett et al. (1987). The baths are composed of sets of harmonic oscillators interacting linearly with the spin, and their terms in $U$ and $U^{\dagger}$ as well as $\rho_{H}^{eq}$, $\rho_{C}^{eq}$ and $e^{\pm i\nu H_{C}}$ can be represented as standard path integrals, which can be integrated out as many Gaussians Feynman and Vernon (1963). The functional $G_{if}(\nu)$ can hence be represented as as a double path integral of the spin paths $X$ and $Y$ weighted by an action, i.e. as $e^{\frac{i}{\hbar}{\cal A}[X,Y]}$. At $\nu=0$ this is the same spin-boson path integral derived in Leggett et al. (1987), which represents the quantum operation of moving the density matrix of the spin at time zero to the density matrix of the spin at time $t$. For non-zero values of $\nu$ additional terms appear in ${\cal A}$, details are summarized in Appendix B.

In practice the spin-boson path integrals are quite cumbersome to do without replying on the “non-interacting blip approximation” (NIBA). The terms in ${\cal A}$ that arise from integrating out the bath(s) are double integrals with kernels, and NIBA means that those kernels should have short enough memory. More precisely, memory should be shorter than the duration of the periods when $X$ and $Y$ take the same value, $(\frac{1}{2},\frac{1}{2})$ or $(-\frac{1}{2},-\frac{1}{2})$, so that the bath can only remember the preceding such period. Since the switching rate of paths in the double path integral is given by the tunneling rate in the system Hamiltonian, NIBA is hence expected to hold when that tunneling rate is small. The same reasoning essentially holds for non-zero values of $\nu$. The set-up is summarized in Appendices A and B.

With caveats discussed in Appendix H the stationary state (for the spin) in the spin-boson problem can then (within NIBA) be determined by almost classical arguments. A transition from the up state $(\frac{1}{2},\frac{1}{2})$ to the down state $(-\frac{1}{2},-\frac{1}{2})$ proceeds through two channels labeled by which spin path goes first ($X$ or $Y$), and the time ($\Delta t$) spent in the intermediate “blip” state ($(\frac{1}{2},-\frac{1}{2})$ or $(-\frac{1}{2},\frac{1}{2})$). The first jump occurs with intrinsic rate $i\frac{\Delta}{2}$ or $-i\frac{\Delta}{2}$ and the second jump with the other rate. Altogether, for both kinds of channels, this gives $\frac{\Delta^{2}}{4}$. The two baths are in equilibrium with respect to the spin before the jump, and integrating them out thus leads to characteristic functions $S_{C}$ and $X_{C}$ for the cold bath and $S_{H}$ and $X_{H}$ for the hot bath. Summing contributions from all channels thus gives an overall transition rate from up to down:

[TABLE]

and a similar overall transition rate from down to up

[TABLE]

The stationary probability to be up is $\frac{D}{A+D}$. This expression is formally identical with the dimensional arguments in Section III: $A$ may be identified with $\Gamma_{\uparrow\downarrow}$; and $D$ with $\Gamma_{\downarrow\uparrow}$ 111The sum $A+D$ is proportional to the quantity called $\tilde{g}$ in Leggett et al. (1987) (at zero Laplace transform parameter), and the difference $D-A$ is proportional to $-\tilde{h}$. The magnetization is $\frac{D-A}{A+D}$, which equals $-\tilde{h}/\tilde{g}$ in the notation of Leggett et al. (1987). .

The calculations of thermal power detailed in Appendices E-J rely crucially on exact relations between the derivative of the action ${\cal A}$ with respect to the parameter $\nu$ at $\nu=0$, and the derivatives of the two functions $S$ and $X$ with respect to the time argument. It is then convenient to introduce additional characteristic functions of the hot and the cold baths 222Equivalent functions have been introduced in the previous literature, but not exactly for these quantities

[TABLE]

The quantity $A$ introduced above is then

[TABLE]

and similarly for $D$.

As determined in appendix, the rate of energy change in the cold bath while the system is respectively in the up and the down state can be written, compare (116),

[TABLE]

An interpretation of the above results is that $C_{+}^{C}$, $C_{+}^{H}$, $C_{-}^{C}$ and $C_{-}^{H}$ are the influence functionals from integrating out the baths when the forward and backward paths of the spin are fixed and opposite. These influence functionals are of the form $\hbox{Tr}\left[U\rho^{eq,\uparrow}V^{\dagger}\right]$ with different unitary operators applied to the left and to the right. Differentiating $U$ and $V^{\dagger}$ with respect to time brings down $-\frac{i}{\hbar}\left(H_{B}+H_{I}\right)$ and $\frac{i}{\hbar}\left(H_{B}+H_{I}^{\prime}\right)$ with different interaction Hamiltonians on the two sides because the spin coordinate is different on the two sides. The bath Hamiltonians are however the same and their contributions hence cancel, and the remaining terms are expectation values of the interaction Hamltonians, conditional on which state the spin started from, which path jumped first, and the blip duration. In this way (25) and (26) can be seen to give an estimate of the type of (11).

V The polaron transform picture

Another interpretation of the results in (25) and (26) is based on the polaron transform. Changing $\hat{\sigma}_{z}$ from up to down has the same effect on the bath energy as instantaneously shifting the position of every bath oscillator $q_{b}$ by an amount $2\frac{C_{b}}{m_{b}\omega_{b}^{2}}$. Such a shift is generated by $\hat{B}_{+}=\exp\left(i2\sum_{b}\frac{C_{b}}{\hbar m_{b}\omega_{b}^{2}}\hat{p}_{b}\right)$ where $\hat{p}_{b}$ is the momentum operator of oscillator $b$. Similarly $\hat{B}_{-}=\exp\left(-i2\sum_{b}\frac{C_{b}}{\hbar m_{b}\omega_{b}^{2}}\hat{p}_{b}\right)$ has the same effect on the bath energy as changing $\hat{\sigma}_{z}$ from down to up.

The function $C_{-}(t)=e^{-\frac{1}{\hbar}S-\frac{i}{\hbar}X}$ for the cold or hot bath ($C$ or $H$) is therefore the same as $\left<\hat{B}_{-}(t)\hat{B}_{+}(0)\right>_{eq}$ where the operators are in Heisenberg picture, and the average is over the bath in equilibrium. Similarly $C_{+}(t)=e^{-\frac{1}{\hbar}S+\frac{i}{\hbar}X}$ is the same as $\left<\hat{B}_{-}(0)\hat{B}_{+}(t)\right>_{eq}$. The effective jump rates are thus

[TABLE]

and similarly for $D$. The above may be used to derive the weak-interaction limit, since then $\hat{B}_{+}\approx\mathbf{1}+i2\sum_{b}\frac{C_{b}}{\hbar m_{b}\omega_{b}^{2}}\hat{p}_{b}$ and $\hat{B}_{-}\approx\mathbf{1}-i2\sum_{b}\frac{C_{b}}{\hbar m_{b}\omega_{b}^{2}}\hat{p}_{b}$, and (linear terms cancel)

[TABLE]

Except for $\epsilon$ very small this gives the effective jump rate proportional to the sum of the spectral powers of the cold and hot bath at frequency $\epsilon/\hbar$, which can be compared e.g. to Karimi and Pekola (2016) (Eq. 3).

In a similar manner one may also consider (25) and (26). The derivatives $\frac{dC_{+}^{C}(t)}{dt}$ and $\frac{dC_{-}^{C}(t)}{dt}$ translate (in weak coupling) to $\frac{4}{\hbar^{2}}\sum_{b}\frac{C_{b}^{2}}{m_{b}^{2}\omega_{b}^{4}}\left<\hat{p}_{b}(0)\frac{d\hat{p}_{b}(t)}{dt}\right>_{C,eq}$ and $\frac{4}{\hbar^{2}}\sum_{b}\frac{C_{b}^{2}}{m_{b}^{2}\omega_{b}^{4}}\left<\frac{d\hat{p}_{b}(t)}{dt}\hat{p}_{b}(0)\right>_{C,eq}$. The dependence on the hot bath is only to higher orders in the interaction coefficients, and therefore drops out. Given that $C_{+}^{C}(0)$ and $C_{-}^{C}(0)$ are both equal to one, one may integrate by parts, which gives

[TABLE]

which is of the same form as (8) and (9).

VI Discussion

In this paper we have considered thermal power (heat per unit time) through a qubit interacting with two or several baths as in the spin-boson problem Leggett et al. (1987). By an extension of the Feynman-Vernon influence functional method it is possible to compute the distribution of energy changes in a bath or baths of harmonic oscillators interacting with a general quantum system Aurell and Eichhorn (2015); Aurell (2018a); Funo and Quan (2018); Carrega et al. (2016). Here we have adapted this approach to the situation where the system in one spin.

The advantage of the Feynman-Vernon method is that while each oscillator in the bath is only perturbed slightly, and the system-bath interaction hence assumed linear in the harmonic oscillator coordinates, the accumulated effect on the system from all the bath oscillators can be large. A Feynman-Vernon theory of energy changes in a bath is thus a way to model quantum heat in a system interacting strongly with its environment. In this paper we have only considered the expected value, but in principle higher moments can also be computed e.g. by the formulae given in Aurell (2018a). Furthermore we have only considered the stationary case (constant drive) and the long-time limit which can be analyzed by Laplace transforms, as was already done in Leggett et al. (1987).

If an assumption analogous to the “non-interacting blip approximation” (NIBA) is made, the general structure of the answer is quite simple, and basically follows by dimensional arguments. It can also be expressed in terms of correlation functions and time derivatives of correlation functions after a polaron transform. While the final result is simple, the intermediate calculations are not, as seems to be the case for most path integral treatments of the spin-boson problem, compare Leggett et al. (1987) as well as the later literature Weiss (2012); Grifoni and Hänggi (1998); Grifoni et al. (1997); Hartmann et al. (2000). For the quantum state a much simpler approach is possible using the polaron transform directly Aslangul, C. et al. (1986); Dekker (1987). Since our result for thermal power can also be expressed in terms of quantities after a polaron transform, it would be interesting to know if it can also be found in a simpler manner. We leave this question to future work, as well as numerical determination terms (25) and (26) in thermal power.

We end by noting that for a qubit interacting with two baths the prediction of NIBA may be not only incorrect, but also physically inadmissable. The limits of validity of NIBA may thus be qualitatively different in non-equilibrium compared to equilibrium. This question deserves further study. We further note that in NIBA the condition that thermal power to the cold bath be positive appears different than the admissibility condition on the state. Conceivably there may hence be situations where NIBA is appropriate, for the quantum state but not for quantum thermodynamics. This issue also deserves further study.

Acknowledgments

This work was supported by ESPCI Chaire Joliot 2018 (EA). EA thanks Jukka Pekola and Bayan Karimi for many discussions on heat flows in superconducting devices, Yuri Galperin for a critical reading of the MS, and Dmitry Golubev for showing results prior to publication. Results equivalent to Eq. (25) and (26) have also been derived independently by Golubev in the case of zero bias. FM was supported by an Erasmus+ Student Mobility for Traineeship (Politecnico di Torino, Italy), and thanks Nordita (Stockholm, Sweden) for hospitality.

Appendix A Summary of spin-boson theory and NIBA

The calculations in Section IV are for the quantum thermal power and two baths what Leggett and collaborators did in the 80ies for the development of the quantum state and one bath Leggett et al. (1987). This Appendix summarizes relevant results from that earlier calculation. For ease of comparison (here and in later related Appendices) we follow the notation of Leggett et al. (1987). We restate the system (qubit) Hamiltonian:

[TABLE]

where $\Delta$ is a rate (dimension $\left(\hbox{time}\right)^{-1}$), and $\epsilon$ is the level splitting. The bath Hamiltonians are, in classical notation,

[TABLE]

where the parameters $m_{b}$ and $\omega_{b}$ are the mass and angular frequency of each oscillator and $C$ and $H$ also stand for the sets of oscillators in respectively the cold bath and the hot bath. The system-bath interactions are similarly

[TABLE]

where $C_{b}$ is the interaction coefficient between bath oscillator $b$ and the qubit, and $\hat{\sigma}_{z}$ operates on the qubit. The coupling coefficients $C_{b}$ have dimension $(\hbox{energy})\cdot(\hbox{length})^{-1}$.

The Feynman-Vernon transition probability of a general quantum system interacting with two baths is

[TABLE]

where the initial state of the baths $\rho_{CH}^{eq}$ is the product state of two thermal states $\rho_{C}^{eq}$ and $\rho_{H}^{eq}$, at two temperatures. $U$ is the big unitary expressing the forward time evolution due to the total Hamiltonian given by (31), (32) (33) (34) and (35), and $U^{\dagger}$ (the adjoint) is the backward time evolution.

The bath coordinates in (36) can be integrated out to yield

[TABLE]

where $X$ is the system coordinate in the forward system path (part of the representation of $U$), $Y$ is the system coordinate in the backward system path (part of the representation of $U^{\dagger}$), and $\int_{if}$ denotes the projection on initial and final states (integrals over initial and final positions of the system in the forward and backward path). The result of integrating out the cold bath is $\frac{i}{\hbar}S^{C}_{i}[X,Y]-\frac{1}{\hbar}S^{C}_{r}[X,Y]$, and the result of integrating out the hot bath is $\frac{i}{\hbar}S^{H}_{i}[X,Y]-\frac{1}{\hbar}S^{H}_{r}[X,Y]$. The real terms ($S_{r}$) depend on the difference $X-Y$ at two different times while the imaginary terms ($S_{i}$) depend on the difference $X-Y$ at a later time, and the sum $X+Y$ at an earlier time.

For the system and bath interaction described by (31), (32) (33) (34) and (35) the system paths $X$ and $Y$ can be represented as piece-wise constant, taking value $\frac{1}{2}$ when the spin is up, and $-\frac{1}{2}$ when the spin is down. This means that at any one one time the forward-backward system path pair can take only four values $(\frac{1}{2},\frac{1}{2})$, $(-\frac{1}{2},-\frac{1}{2})$, $(\frac{1}{2},-\frac{1}{2})$ and $(-\frac{1}{2},\frac{1}{2})$. The two first are in the terminology of Leggett et al. (1987) called sojourns and correspond to populations, the diagonal elements of the density matrix, up and down. The last first are in the terminology of Leggett et al. (1987) called blips and correspond to coherences, the off-diagonal elements of the density matrix. The kind of sojourn and blip can be indicated by variables $\chi=X+Y$ and $\xi=X-Y$, both taking values $\pm 1$. A given double path in $X$ and $Y$, conventionally starting from the up sojourn, can therefore be represented as

[TABLE]

where $\Delta t_{0},\Delta t_{2},\ldots$ are the durations of the sojourns and $\Delta t_{1},\Delta t_{3},\ldots$ are the durations of the blips. The first sojourn starts at time $t_{0}$ and the $n$’th sojourn at time $t_{2n}=t_{0}+\sum_{j=0}^{2n-1}\Delta t_{j}$; the $n$’th blip starts at time $t_{2n+1}=t_{0}+\sum_{j=0}^{2n}\Delta t_{j}$.

The $\hat{\sigma}_{x}$ terms in (31) translate to weights in the integrations ${\cal D}X$ and ${\cal D}Y$ in (37) which are $+i(\frac{\Delta}{2})$ if the forward path ($X$) jumps, and $-i(\frac{\Delta}{2})$ if the backward path ($Y$) jumps. Everything else is included in the total exponent in (37) which one can write as

[TABLE]

where all terms are integrals over time of the terms in the exponent in (37). The first line in above hence represent the terms $\frac{i}{\hbar}S_{S}[X]-\frac{i}{\hbar}S_{S}[Y]$ which have only one time integral, and which are non-zero only for blips, the terms $-\frac{1}{\hbar}S^{C}_{r}[X,Y]-\frac{1}{\hbar}S^{H}_{r}[X,Y]$, with both terms in the same blip, and $\frac{i}{\hbar}S^{C}_{i}[X,Y]+\frac{i}{\hbar}S^{H}_{i}[X,Y]$ with the sojourn immediately before the blip. The second and third line in (39) correspond to times separated by at least one sojourn.

The Non-interacting blip approximation (NIBA) of Leggett et al. (1987) is to ignore the second and third line of (39), and to assume that $X^{C}_{j,j-1}$ and $X^{H}_{j,j-1}$ only depend on the associated blip duration $\Delta t_{2j-1}$. The validity of this approximation was discussed in depth in Leggett et al. (1987) and in the later literature, see e.g Weiss (2012); Grifoni and Hänggi (1998); Grifoni et al. (1997); Hartmann et al. (2000). Here we only note that it is essentially an expansion in small tunneling rates $\Delta$, as lucidly explained in Aslangul, C. et al. (1986) and Dekker (1987), with long blip durations suppressed as a result of the interaction between the system and the baths.

The content of NIBA is thus expressed in the following two characteristic functions of the baths, which we write for the cold bath as

[TABLE]

In above the sums are over oscillators in the cold bath and $\beta_{C}$ is the inverse temperature of the cold bath. The formulas for the contributions from the hot bath are analogous.

It is customary to also write the above functions as $Q_{1}$ and $Q_{2}$ as these are equivalent in NIBA. If one does not assume NIBA, $X_{j,j-1}$ would however be the sum of three terms $Q_{1}$ with different arguments, where the one above is the shortest time.

Appendix B Heat and NIBA

The starting point is the generating function of energy changes in the cold bath

[TABLE]

This equation is the same as (36) above, except that exponentials of the Hamiltonian of the cold bath have been inserted at the initial and final time. It is assumed in (42) that $e^{-i\nu H_{C}}$ commutes with the initial density matrix of the baths $\rho_{CH}^{eq}$; this issue, related to strong coupling, will be discussed below.

As for (36) we can introduce path integral representations of $U$ and $U^{\dagger}$ and integrate out the two baths. The result must analogously to (37) look like

[TABLE]

where the two new functionals ${\cal J}_{\nu}$ and $\tilde{\cal J}_{\nu}$, which represent the distribution of energy changes in the cold bath, are quadratic in $X$ and $Y$. The two terms are for later convenience separated as to and respectively depending anti-symmetric and symmetric combinations in the exchange of times. In earlier contributions the same two functionals and their kernels were distinguished by superscripts $(2)$ and $(3)$ Aurell and Eichhorn (2015); Aurell (2017, 2018b). Here we choose to streamline the formalism, additionally because the similar functional with superscript $(1)$ does not appear; for a discussion, see Aurell (2017).

When $\nu$ is equal to zero $G_{if}(\nu)$ is equal to $P_{if}$, and the two functionals ${\cal J}_{\nu}$ and $\tilde{\cal J}_{\nu}$ must vanish. In this paper we are concerned with the terms linear in $\nu$ which are given by

[TABLE]

with two kernels

[TABLE]

These two kernels are the same as $h^{(2)}$ and $h^{(3)}$ in Aurell and Eichhorn (2015), except for a factor $\hbar$.

It is a non-trivial fact Aurell and Eichhorn (2015) that $h$ and $\tilde{h}$ are proportional to time derivatives of the Feynman-Vernon kernels

[TABLE]

Similar relations between second integrals of these kernels will be crucial in the following.

We can now represent $G_{if}(\nu)$ in a similar way to (39) with new terms stemming from $J$ and $\tilde{J}$. We can write these as

[TABLE]

In above $X^{(1)}_{jk}$ are the first-order terms in $\nu$ from the kernels anti-symmetric in the time exchange. In contract to the imaginary Feynman-Vernon kernel, both the blip-sojourn and sojourn-blip terms appear. Furthermore $S^{(1)}_{j}$ and $S^{(1^{\prime})}_{j}$ are the first-order terms in $\nu$ from the kernels symmetric in the time exchange where both times fall in the same time interval. In contrast to the real Feynman-Vernon kernel, there are such terms from both blips and sojourns. Finally $\Lambda^{(1)}_{jk}$ and $\Lambda^{(1^{\prime})}_{jk}$ are terms from two intervals of the same kind, either two blips or two sojourns.

A NIBA-like approximation to (50) means to include only the terms from an adjacent blip and sojourn. These are on the one hand terms like $-\frac{1}{2}\nu\xi_{j}\chi_{j-1}X^{(1)}_{j,j-1}$, and on the other $\frac{1}{2}\nu\xi_{j}\chi_{j}X^{(1)}_{j,j}$ both of which depend on time increments as discussed for $X_{j,j-1}$ above. Only one of these time increments is for a blip interval (the same blip interval), and we are therefore led to

[TABLE]

From this we have the NIBA-like approximation

[TABLE]

Comparing to (49) and (41) we see that

[TABLE]

A NIBA-like approximation to (51) is a bit more involved, for two reasons. First the two terms $S^{(1)}_{j}$ and $S^{(1^{\prime})}_{j}$ both need to be included, and they are both diverging in the bath cut-off frequency. This requires a separate discussion which we give below in Appendix E. Second, the terms on the second line of (51) cannot be neglected entirely. This is so because the interaction of two neighboring sojourns ($\frac{1}{2}\nu\chi_{j}\chi_{j-1}\Lambda^{(1^{\prime})}_{j,j-1}$) has one terms which depends on the intervening blip time, and which hence gives

[TABLE]

Comparing to (40) we see that

[TABLE]

Appendix C Ohmic baths

Ohmic baths have spectra (density of states) that are continuous up to some very large upper cut-off $\Omega$ and increase quadratically with frequency. The number of oscillators with frequencies in the interval $[\omega,\omega+d\omega]$ is $f(\omega)d\omega$ can then be taken to be

[TABLE]

where $\omega_{c}$ is some characteristic frequency less than $\Omega$. The total number of oscillators is then $\frac{2}{3\pi}\left(\frac{\Omega}{\omega_{c}}\right)^{3}$.

An alternative version is to take a smooth cut-off:

[TABLE]

In this case the number of bath oscillators is $\frac{12}{\pi}\left(\frac{\Omega}{\omega_{c}}\right)^{3}$.

The system-bath interactions are characterized by two parameters $\eta_{C}$ and $\eta_{H}$ such that for an oscillator in the cold bath

[TABLE]

and for an oscillator in the hot bath

[TABLE]

For the spin-coupling problem the dimensions of $\eta_{C}$ and $\eta_{H}$ are $(\hbox{mass})\cdot(\hbox{length})^{2}\cdot(\hbox{time})^{-1}$ i.e. the action.

The terms $X_{j,j-1}(\tau)$ and $S_{j}(\tau)$ in (40) and (41) were computed in Leggett et al. (1987) as $\eta\tan^{-1}(\Omega\tau)$, and $\frac{1}{2}\eta\log(1+\Omega^{2}\tau^{2})+\eta\log\left(\frac{\hbar\beta}{\pi\tau}\sinh\frac{\pi\tau}{\hbar\beta}\right)$. The first is essentially a sign function. The second starts as $\frac{\eta}{2}\Omega^{2}\tau^{2}$ in the interval $\tau<<\Omega^{-1}$, then grows as $\eta\log\Omega\tau+\frac{\eta}{2\Omega^{2}\tau^{2}}$ in the interval $\Omega^{-1}<<\tau<<\hbar\beta$ and finally behaves as $\eta\log\frac{\Omega\hbar\beta}{2\pi}+\frac{\eta\pi}{\hbar\beta}|\tau|$ when $\tau>>\hbar\beta$. The derivative $\partial_{\beta}S_{j}(\tau)$ evaluates to $\eta/\beta(1-\frac{\pi\tau}{\hbar\beta}\coth\frac{\pi\tau}{\hbar\beta})$. which is always negative. $S_{j}(\tau)$ is hence an increasing function of bath temperature. The second derivative $\partial_{\beta\tau}S_{j}(\tau)$ evaluates to $\frac{\pi\eta}{\hbar\beta^{2}}(-\coth\frac{\pi\tau}{\hbar\beta}+\frac{\pi\tau}{\hbar\beta}\sinh^{-2}\frac{\pi\tau}{\hbar\beta})$. which is also always negative. $\partial_{t}S_{j}(\tau)$ is hence also an increasing function of bath temperature.

$K$ and $\tilde{K}$ can be computed from (54) and (56): $\tilde{K}$ is essentially a delta function on the bath cut-off frequency scale $\Omega^{-1}$, while $K$ is basically a delta function on the time scale $\hbar\beta$, and for large $\tau$ a constant.

Appendix D Interaction through bath momentum

Theorem D.1

Let a system described by coordinate $X$ interact with by a bath of harmonic oscillators described by coordinate and momenta $(q_{b},p_{b})$ through a combined bath and interaction Hamiltonian $\sum_{b}\frac{1}{2m_{b}}\left(p_{b}+m_{b}C_{b}X\right)^{2}+\frac{1}{2}m_{b}\omega_{b}q_{q}^{2}$. The coupling coefficients $C_{b}$ vanish at the beginning and the end of the process. Then the generating function of the change of bath energy is the same is if the combined bath and interaction Hamiltonian would have been $\sum_{b}\frac{1}{2m_{b}}p_{b}^{2}+\frac{1}{2}m_{b}\omega_{b}^{2}\left(q_{q}-\frac{C_{b}}{\omega_{b}}X\right)^{2}$.

The proof proceeds by adapting the calculation in Aurell (2018a), in the following steps.

The action corresponding to the Hamiltonian coupled through momentum is $\int\frac{1}{2}m_{b}\dot{q}^{2}_{b}-m_{b}C_{b}X\dot{q}_{b}-\frac{1}{2}m_{b}\omega_{b}q_{q}^{2}$. By an integration by parts the term linear in $\dot{q}_{b}$ is changed to $\hbox{{boundary terms}}+\int m_{b}\frac{d}{dt}(C_{b}X)q_{b}$. 2. 2.

The path integral of the bath oscillator with fixed initial and final positions can then be considered to be that of a Lagrangian $\int\frac{1}{2}m_{b}\dot{q}^{2}_{b}-\frac{1}{2}m_{b}\omega_{b}^{2}q_{b}^{2}+m_{b}\frac{d}{dt}(C_{b}X)q_{b}$. This path integral can then be done as in Feynman-Vernon theory giving integrals of the external drive (here $m_{b}\frac{d}{dt}(C_{b}X)$) multiplying the initial and final positions of the oscillator, and a constant. 3. 3.

The integrals are of the type ($u$ in the notation of Aurell (2018a), Appendix A) $\frac{1}{\sin\omega_{b}t}\int^{t}_{0}\sin\omega_{b}(t-s)m_{b}\frac{d}{ds}(C_{b}X)(s)ds$. By a partial integration they can be combined with the boundary terms to give $\frac{m_{b}\omega_{b}}{\sin\omega_{b}t}\int^{t}_{0}\cos\omega_{b}(t-s)C_{b}Xds$, multiplying the initial position of the bath oscillator in the forward path. There are four terms of this type with two sign changes compared to Aurell (2018a), Appendix A. 4. 4.

The constant ($B$ in the notation of Aurell (2018a), Appendix A) is two terms of the type $\frac{1}{m_{b}\omega_{b}\sin\omega_{b}t}\int^{t}_{0}\int^{s}\sin\omega_{b}(t-s)\sin\omega_{b}s^{\prime}m_{b}\frac{d}{ds}(C_{b}X)(s)m_{b}\frac{d}{ds^{\prime}}(C_{b}X)(s^{\prime})ds^{\prime}ds$. By two integrals by parts the sines are turned into cosines multiplying $(C_{b}X)(s)(C_{b}X)(s^{\prime})$, and there is a change of sign. Additionally there is a boundary term $-\frac{m_{b}}{2}\int C_{b}^{2}X^{2}ds$, the same as appears in the complete square $-\frac{1}{2}m_{b}\omega_{b}^{2}\left(q_{q}-\frac{C_{b}}{\omega_{b}}X\right)^{2}$. 5. 5.

The integration over the initial and final coordinates of the bath oscillator proceeds as in Aurell (2018a), Appendix A, and gives in fact the same result, with $m_{b}C_{b}\omega_{b}$ appearing instead of $C_{b}$. One of the authors (E.A.) points out that there is an error in Eqs (25) and (A14) in Aurell (2018a): the constant appearing in the kernel ${\cal J}^{(2)}$ should read $(yz^{\prime}-y^{\prime}z)/\Delta$ (instead of $(y^{\prime}z^{\prime}-yz)/\Delta$). To linear order in the parameter $\nu$ these two quantities are however the same, hence there is no difference to the present paper.

In summary, the only difference to coupling through coordinate is hence that if the coupling coefficient to bath momentum is $C$, then the equivalent coupling coefficient to bath coordinate is $m\omega C$, as is also required dimensionally.

Appendix E The singular NIBA heat terms

In this appendix we estimate the contributions $S^{(1)}_{j}$ and $S^{(1^{\prime})}_{j}$ to (51). Both these terms are second integrals of the kernel $\tilde{h}$ in (47) over one blip or one sojourn interval, hence proportional to

[TABLE]

For an Ohmic bath with sharp cut-off this expression is $\frac{2\eta}{\pi}\left(\Omega-\delta_{\Omega}(\Delta t)\right)$ where $\delta_{\Omega}(\Delta t)$ a delta-function smoothened at time scale $\Omega^{-1}$. The contribution to $G_{if}(\nu)$ from $n+1$ sojourns and $n$ blips is hence

[TABLE]

While the first line sums to a large number it does not scale with the time, and there will hence not be any contribution to thermal power from these terms.

The large terms are in fact an artifact from assuming that the baths are in equilibrium at the start and the end of the process while still interacting strongly with the system. It has been known for quite some time that this leads to problems already for the open quantum system state Grabert et al. (1984); Ford et al. (1985); Rosenau da Costa et al. (2000); Ingold et al. (2009). One way to resolve the problem for heat is to assume that the interaction coefficients $C_{b}$ depend on time, and vanish in the beginning of the process Aurell (2017). Assuming as in Aurell (2017) and in analogy with (59) above that $C_{b}(s)=\sqrt{\omega_{c}^{3}m_{\omega}\eta(s)}$ we have instead of above

[TABLE]

In above the bath cut-off frequency has been taken to infinity. Clearly if the function $\eta(s)$ is constant except at the boundaries this does not give anything proportional to the duration of the process.

Appendix F The non-singular NIBA heat terms: general formalism

The main idea is to write the sum $G_{i}(\nu)=\sum_{f}G_{if}(\nu)$ as a matrix product (transfer matrix formalism). The formulation is as follows:

Starting state $i$ is by convention “up”. The starting vector is therefore \chi_{0}=\left(\begin{array}[]{c}1\\ 0\end{array}\right)=(\uparrow,\uparrow). 2. 2.

End vector, when we sum over the final state of the system, is \chi_{n}=\left(\begin{array}[]{c}1\\ 1\end{array}\right)=(\uparrow,\uparrow)+(\downarrow,\downarrow). 3. 3.

The phase terms at the jumps are determined by the translation tables

[TABLE]

and

[TABLE] 4. 4.

To every transition sojourn $\to$ blip are associated terms $e^{\frac{i}{\hbar}\chi_{j-1}\xi_{j}\left(X_{j.j-1}+\frac{1}{2}\nu K\right)}$. Combine this and the phase factor to a matrix $\frac{i\Delta}{2}\mathbf{T}$. 5. 5.

To every blip interval is associated the terms $e^{-\frac{1}{\hbar}S_{j}-\frac{i}{\hbar}\epsilon(t_{2j}-t_{2j-1})}$. Call this diagonal matrix $\mathbf{\Lambda}$. 6. 6.

To every transition blip $\to$ sojourn is associated a term $e^{-\frac{i}{\hbar}\chi_{j}\xi_{j}\nu\frac{1}{2}K}$ Combine this and the phase factors to a matrix $\frac{i\Delta}{2}\mathbf{S}$. 7. 7.

To every transition sojourn $\to$ sojourn is additionally associated as term $e^{\frac{i}{\hbar}\chi_{j}\chi_{j-1}\nu\frac{1}{2}\tilde{K}}$. This is the same for both signs of the blip in between. 8. 8.

The transition sojourn $\to$ sojourn is then given by a matrix $\mathbf{M}$ formed by $\mathbf{S}\mathbf{\Lambda}\mathbf{T}$ and the modifications due to $\tilde{K}$. By matrix multiplication one finds

[TABLE]

For simplicity the blip interval is written $t$. 9. 9.

The whole generating function can hence, within NIBA, be written as

[TABLE]

where all the blip times are implicit in the matrices $\mathbf{M}$ on the right-hand side.

To analyze (65) in a stationary setting (the bias $\epsilon$ and all other parameters are constant in time) one takes a Laplace transform. Every sojourn interval then yields a factor $\lambda^{-1}$, and the $n$’th term in (65) hence a factor $\lambda^{-1-n}$. For the Laplace transform of the matrix it is convenient to write

[TABLE]

where

[TABLE]

All $S$, $X$, $K$ and $\tilde{K}$ depend on the blip time $t$ (at least in principle).

The Laplace transform of the generating function is

[TABLE]

Appendix G The generating function at $\nu=0$

The special case of $\nu=0$ is an important check, because that should give the quantity computed by Leggett in Leggett et al. (1987): $\tilde{P}(\lambda)=\int dte^{-\lambda t}\left<\sigma_{z}\right>(t)$. The relation is $\left<\sigma_{z}\right>(t)=2\cdot\hbox{Prob}(\hbox{"up"},t)-1$ and hence $\tilde{P}(\lambda)=\tilde{G}_{if}(\nu=0,\lambda)-\lambda^{-1}$ where $i$ and $f$ are both “up”. The formula found by Leggett is

[TABLE]

where

[TABLE]

We hence consider (74) at $\nu=0$. We have the simplification that $C=A$ and $B=D$, and the Laplace transform matrix is hence

[TABLE]

The eigenvalues of this matrix are [math] and $2(A+D)$. Positive powers of this matrix ($n\geq 1$) are thus simply

[TABLE]

which means that

[TABLE]

We may identify $\frac{\Delta^{2}}{2}(A+D)=\tilde{g}$ and $\frac{\Delta^{2}}{2}A=\frac{1}{2}(\tilde{g}+\tilde{h})$ and so

[TABLE]

This means that

[TABLE]

which is (75), as required. The result $G_{i}(\nu=0,t)=1$ (normalization of the system state) follows from \left(\begin{array}[]{ll}1&1\end{array}\right)\tilde{\mathbf{M}}=0, which means that ${G}_{i}(\nu=0,\lambda)=\lambda^{-1}$ (only $n=0$ term survives).

Appendix H The long term limit of the generating function at $\nu=0$

On physical grounds it is reasonable to assume that for long times the generating function is

[TABLE]

where $p$ is the long term limit of the probability to be up, and $q_{k}$ and $r_{k}$ are some constants. The Laplace transform is then

[TABLE]

from which follows

[TABLE]

Inserting (80) we have

[TABLE]

where in the integrals defining $A$ and $D$ the Laplace transform parameter $\lambda$ is zero.

A physical density matrix of the qubit must lie inside the Bloch sphere. A necessary condition for $\frac{D}{A+D}$ and $\frac{A}{A+D}$ to be the diagonal elements of a stationary density matrix in the long-time limit is hence that they fall between zero and one. For a qubit interacting with one bath at one temperature this was shown to be always the case in Leggett et al. (1987), even when the density matrix computed under these assumption of NIBA is not correct.

For our case of one qubit interacting with two baths the situation is more involved, and we state it as

Theorem H.1

Consider $S=S_{C}+S_{H}$ and $X=X_{C}+X_{H}$ as an even and an odd function on the whole line. Let $\hat{F}(\omega)$ be the Fourier transform of $e^{-\frac{1}{\hbar}S+\frac{i}{\hbar}X}$ and $\hat{F}^{*}(\omega)=\hat{F}(-\omega)$ the Fourier transform of $e^{-\frac{1}{\hbar}S-\frac{i}{\hbar}X}$. Then $\frac{D}{A+D}$ and $\frac{A}{A+D}$ are possible diagonal elements of a density matrix if $|{\cal I}\hat{F}(\frac{\epsilon}{\hbar})|<|{\cal R}\hat{F}(\frac{\epsilon}{\hbar})|$.

The proof is by simple translation. We may write

[TABLE]

and the condition

[TABLE]

is hence the same as

[TABLE]

Multiplying out and identifying terms says that the imaginary part of the Fourier transform should be smaller in absolute value than the real part, at the frequency of the level splitting. Note that the theorem does not give a condition for NIBA with two baths to be correct, only a condition for it to give physically admissible populations.

With two caveats one may interpret (86) in an almost classical manner. First we can (trivially) rewrite it as

[TABLE]

where (at $\lambda=0$)

[TABLE]

The two terms in $\frac{\Delta^{2}}{2}A$ are the integrals over time $t$ of the influence functionals of two particular spin histories, where the state is $(\uparrow,\uparrow)$ before time zero, then at time zero either the forward or the backward path jumps to down, and then at time $t$ the other path follows. The two terms $\left(i\frac{\Delta}{2}\right)$ and $\left(-i\frac{\Delta}{2}\right)$ are the jump rate amplitudes (dimension $\left(\hbox{time}\right)^{-1}$) for the two paths. These combined with the integral over time $t$ hence gives a quantitity analogous to the probability that the state transits from $(\uparrow,\uparrow)$ to $(\downarrow,\downarrow)$ per unit time. The two terms in $\frac{\Delta^{2}}{2}D$ may similarly be taken to represent the total rate of the state transiting from $(\downarrow,\downarrow)$ to $(\uparrow,\uparrow)$.

The first of the two caveat is that by the above $A$ and $D$ may have different signs so that one of $\frac{A}{A+D}$ and $\frac{D}{A+D}$ is negative, and the other is larger than one. If so, NIBA would not give a physically admissable state. The second is that even when $\frac{A}{A+D}$ and $\frac{D}{A+D}$ are both between zero and one, both $A$ and $D$ could be negative. NIBA would in that case give a physically admissable state, but not one that can be described as from a classical jump process.

Appendix I Derivatives of generating function formula at $\nu=0$

The expected energy change of the bath is given by the derivative of the generating function (74) with respect to $i\nu$ taken at $\nu=0$. At any $\nu$ this quantity is

[TABLE]

At $\nu=0$ the sums on the left and the right simplify as above. On the left only the zeroth order term ($l=0$) survives, while on the right we have

[TABLE]

The dependence on $\nu$ comes either through the function $K$, or the function $\tilde{K}$. In the first case only the off-diagonal elements ($B$ and $C$) depend on $\nu$, and the total expression is

[TABLE]

where \dot{C}=\frac{dC}{d(i\nu)}|_{\nu=0\,\hbox{through K}} and \dot{B}=\frac{dB}{d(i\nu)}|_{\nu=0\,\hbox{through K}}. These derivatives follow from (68) and (69) and are

[TABLE]

Following (54) we can rewrite this as

[TABLE]

In the second case of dependence through $\tilde{K}$ the derivative matrix is

[TABLE]

where we have used (56). Together with (98) we have hence also

[TABLE]

where

[TABLE]

Appendix J Long-time limit of the derivative

On physical grounds it is reasonable to assume that the derivative of the generating function with respect to its argument is for long times

[TABLE]

where $\Pi$ is the long time limit of the power (heat per unit time), and $b$, $c_{k}$ and $\lambda_{k}$ are some constants. The Laplace transform is then

[TABLE]

from which follows

[TABLE]

Inserting the various formulas above we have

[TABLE]

where in the integrals defining $A$ and $D$ the Laplace transform parameter $\lambda$ is zero, and where the subscript $C$ indicates that only the quantities for the cold bath are considered. Clearly we now have an expression for power similar to the dimensional formula (7). For the case of only one bath we can integrate by parts in (110) to get

[TABLE]

which is the expected result. In the long term limit the thermal power from one qubit equilibrating with one bath must vanish. If we were to consider heat to the hot bath, all that would change (110) is that the time derivatives would be $\partial_{t}S_{H}$ and $\partial_{t}X_{H}$. By adding the same argument as in (111) shows that the the sum of thermal power to the cold bath and the hot bath cancel.

In the case of two baths and heat to one bath it is on the other hand more convenient to write $S=S_{C}+S_{H}$ and $X=X_{C}+X_{H}$ and to introduce the kernels 333Similar kernels have been introduced in the literature before, but not exactly these ones; hence the new notation.

[TABLE]

in terms of which (110) can be written

[TABLE]

This is the formulation used in Section IV and Section V in the main text.

Physically, thermal power to the cold bath must be positive. Referring to the discussion at the end of Appendix H we may identify $A$ as $\frac{1}{2}\hat{F}^{*}(\frac{\epsilon}{\hbar})$ and $D$ as $\frac{1}{2}\hat{F}(\frac{\epsilon}{\hbar})$ and the terms in parentheses in (116) as Fourier components of the function $H(t)=i\hbar\frac{dC_{+}^{C}(t)}{dt}C_{+}^{H}(t)$. Thermal power would then be $-\frac{\Delta^{2}}{2{\cal R}[\hat{F}(\frac{\epsilon}{\hbar})]}{\cal R}[\hat{F}(\frac{\epsilon}{\hbar})\hat{H}^{*}(\frac{\epsilon}{\hbar})]$.