Quantum Thermodynamics with Degenerate Eigenstate Coherences

Abstract

We establish quantum thermodynamics for open quantum systems weakly coupled to their reservoirs when the system exhibits degeneracies. The first and second law of thermodynamics are derived, as well as a finite-time fluctuation theorem for mechanical work and energy and matter currents. Using a double quantum dot junction model, local eigenbasis coherences are shown to play a crucial role on thermodynamics and on the electron counting statistics.

1. Introduction

The study of nonequilibrium open quantum systems is an active field of research with particular relevance to routinely-devised systems, such as quantum dots, or electronic circuits [1,2,3,4], or assemblies of cold atoms [5,6,7], for example. The possibility to monitor thermodynamically-relevant quantities, such as heat and work, during a single experimental realization has motivated the study of their fluctuating properties and, thereby, the identification of universal laws satisfied by their statistics [8,9,10].

Quantum master equations have been widely used for the study of the thermodynamic properties of open quantum systems [9,11,12,13]. They are usually derived for systems weakly interacting with their reservoirs using the Born–Markov and secular (BMS) approximation [14,15,16]. The resulting quantum master equation can be shown to be of the Lindblad form [17]. In absence of degeneracies in the system Hamiltonian, the density matrix populations in the system energy eigenbasis satisfy a closed stochastic equation, whereas coherences undergo an independent decay in time [18]. For many processes that only depend on populations or for steady state dynamics where eigenstates coherences are always vanishing, a classical stochastic thermodynamics (ST) [19,20,21] can be easily built for the population dynamics. This provides a consistent framework for the study of the thermodynamics of open quantum systems at both the average and the single trajectory level [22,23,24,25,26,27,28]. However, various time-dependent processes do depend on eigenstate coherences. This happens for instance for systems driven by fast periodic time-dependent forces, where what we just said holds at the level of quasi-energies instead of eigenenergies [13,29]. It also happens for multi-stroke machines or for systems undergoing feedback control, where eigenstate coherences can be shown to play an important thermodynamics role (see, e.g., [30], respectively [31]).

Open quantum systems with degenerate system energies constitute another important case in which eigenstate coherences come into play already in the weak coupling limit. In this case, time-dependent driving is not even required, and coherences may survive even at steady state. Such situations are very important in mesoscopic physics, leading to interesting phenomena, such as maximum output power enhancement [32,33]. The aim of this paper is to extend the central results of stochastic thermodynamics to open quantum systems with degeneracies. When applying the BMS approximations, while the dynamics of populations and coherences between non-degenerate states of the system Hamiltonian remains uncoupled, the populations and coherences between degenerate states remain coupled. We propose consistent definitions for energy, work, heat, entropy and entropy production for such dynamics. We further obtain the counting statistics of the mechanical work and energy and particle currents from the aforementioned quantum master equation and derive a finite time fluctuation theorem that extends its classical counterpart [34] to quantum systems with eigenstate coherences. We illustrate our results on a degenerate double quantum dot system, which exhibits a quantum suppression of the particle current due to coherences [35,36,37,38,39,40,41]. We show that coherences cause a bi-modality in the finite time current distribution [38], which is nevertheless compatible with the fluctuation theorem symmetry.

The paper is organized as follows. The BMS master equation for a general open quantum system with exact degeneracies is exposed in Section 2.1. The analysis of the nonequilibrium thermodynamics is presented in Section 2.2, where we establish the energy and entropy balance, as well as the positivity of entropy production. The thermodynamics analysis is exposed in Section 2.3. An expression for the work and currents statistics is derived in Section 2.3.1 using the dressed quantum master equation formalism [9]. In Section 2.3.2, we prove a finite-time fluctuation theorem for systems described by the quantum master Equation (3). Finally, our approach is applied in Section 3 to study the thermodynamics of a degenerate double quantum dot connected to two electronic leads. A summary is given in Section 4.

2. General Formalism

2.1. Microscopic Derivation of Lindblad Master Equations

We consider an open quantum system with Hamiltonian $H=H_S+H_R+H_I$, in terms of system ($H_S$), reservoir ($H_R$) and interaction ($H_I$) Hamiltonians, respectively. We aim to describe the effective dynamics of the system with a master equation of the form:

$$\begin{array}{c}\hfill {\dot{\rho }}_S=\mathcal{L}{\rho }_S\end{array}$$

(1)

for the system density matrix ${\rho }_S={\mathrm{Tr}}_{\mathrm{R}}\left\{\rho \right\}$ only (here and in the following, ${\mathrm{Tr}}_{\mathrm{R}}\left\{\dots \right\}$ denotes the partial trace over the reservoir degrees of freedom). This equation should preserve the density matrix properties (trace, hermiticity and positivity), at least in an approximate sense. The Lindblad master equation [17] is the most general master equation that preserves the density matrix properties exactly. There are multiple ways of obtaining Lindblad master equations from microscopic Hamiltonians for various parameter regimes [42,43]. Here, we will constrain ourselves to the weak-coupling limit between system and reservoir, in which the Born, Markov and secular (BMS) approximations [15] can be applied, the latter often also termed rotating wave approximation. As such, we will be concerned with systems whose relaxation dynamics is much slower than the fast correlation time of the reservoirs.

Under the aforementioned approximations and for a decomposition of the interaction Hamiltonian:

$$\begin{array}{c}\hfill H_I=\sum _{\alpha }{A}_{\alpha }\otimes B_{\alpha }\end{array}$$

(2)

into system operators $A_{\alpha }$ and reservoir operators $B_{\alpha }$, respectively, the BMS Lindblad master equation becomes for a single reservoir [44]:

$$\begin{array}{c}\hfill {\dot{\rho }}_S(t)=-\mathrm{i}\left[H_S+\sum _{ab}{\sigma }_{ab}{L}_{ab},{\rho }_S(t)\right]+\sum _{ab,cd}{\gamma }_{ab,cd}\left[L_{ab}{\rho }_S(t)L_{cd}^{†}-\frac{1}{2}\left\{L_{cd}^{†}{L}_{ab},{\rho }_S(t)\right\}\right]\equiv \mathcal{L}{\rho }_S(t),\end{array}$$

(3)

where we use the fixed eigen-operator basis $L_{ab}=\left|a\right.〉\left.〈b\right|$ of the system Hamiltonian $H_S\left|a\right.〉=E_a\left|a\right.〉$. We note that this basis is unique when the spectrum of $H_S$ is non-degenerate. Here, the matrix elements of the Lamb shift Hamiltonian ${\sigma }_{ab}={\sigma }_{ba}^{*}$ and the positive definite matrix ${\gamma }_{ab,cd}$ are given by:

$$\begin{array}{ccc}{\sigma }_{ab}\hfill & =\hfill & {\delta }_{E_b,E_a}{\sum }_{\alpha \beta }{\sum }_c\frac{{\sigma }_{\alpha \beta }(E_b-E_c)}{2\mathrm{i}}\left.〈c\right|A_{\beta }\left|b\right.〉\left.〈c\right|A_{\alpha }^{†}{\left|a\right.〉}^{*},\hfill \\ {\gamma }_{ab,cd}\hfill & =\hfill & {\delta }_{E_b-E_a,E_d-E_c}{\sum }_{\alpha \beta }{\gamma }_{\alpha \beta }(E_b-E_a)\left.〈a\right|A_{\beta }\left|b\right.〉\left.〈c\right|A_{\alpha }^{†}{\left|d\right.〉}^{*}.\hfill \end{array}$$

(4)

They depend on the matrix elements of the system coupling operators $A_{\alpha }$ and the even (${\gamma }_{\alpha \beta }$) and odd (${\sigma }_{\alpha \beta }$) Fourier transforms:

$$\begin{array}{c}\hfill {\gamma }_{\alpha \beta }(\omega )=\int C_{\alpha \beta }(\tau )e^{+\mathrm{i}\omega \tau }d\omega ,{\sigma }_{\alpha \beta }(\omega )=\int \mathrm{sgn}(\tau )C_{\alpha \beta }(\tau )e^{+\mathrm{i}\omega \tau }d\omega \end{array}$$

(5)

of the reservoir correlation functions (bold symbols denote the interaction picture ${\mathit{B}}_{\mathit{\alpha }}(\tau )=e^{+\mathrm{i}{H}_R\tau }{B}_{\alpha }{e}^{-\mathrm{i}{H}_R\tau }$):

$$\begin{array}{c}\hfill C_{\alpha \beta }(\tau )=〈{\mathit{B}}_{\mathit{\alpha }}(\tau )B_{\beta }〉={\mathrm{Tr}}_{\mathrm{R}}\left\{{\mathit{B}}_{\mathit{\alpha }}(\tau )B_{\beta }{\overline{\rho }}_R\right\},\end{array}$$

(6)

where ${\overline{\rho }}_R$ denotes the stationary state of the reservoir. For a single reservoir, it is usually chosen as a thermal reference state:

$$\begin{array}{c}\hfill {\overline{\rho }}_R=e^{-\beta (H_R-\mu N_R-{\varphi }_R)}\end{array}$$

(7)

in terms of the reservoir thermodynamic grand-potential ${\varphi }_R=-{\beta }^{-1}\ln \mathrm{Tr}\left\{e^{-\beta (H_R-\mu N_R)}\right\}$. It is characterized by the inverse temperature β and chemical potential μ of the reservoir.

We now summarize a few useful properties of the BMS Lindblad master equation beyond preservation of the density matrix properties.

First, we observe that coherences ${\rho }_{ij}\equiv \langle i|{\rho }_S(t)|j\rangle$ of basis states i and j with different energies $E_i\ne E_j$ will evolve decoupled from the populations ${\rho }_{aa}\equiv \langle a|{\rho }_S(t)|a\rangle$

$$\begin{array}{c}\hfill {\dot{\rho }}_{ij}=-\mathrm{i}\left(E_i-E_j+{\sigma }_{ii}-{\sigma }_{jj}\right){\rho }_{ij}+\sum _{ab}{\gamma }_{ia,jb}{\rho }_{ab}-\frac{1}{2}\sum _{ab}{\gamma }_{ab,ai}{\rho }_{bj}-\frac{1}{2}\sum _{ab}{\gamma }_{aj,ab}{\rho }_{ib},\end{array}$$

(8)

which formally results from the Kronecker-delta functions in (4). This implies that for a non-degenerate system (where $E_i\ne E_j$ implies $i\ne j$), one can directly show that in the system energy eigenbasis, the master equation decouples the evolution of all populations and all coherences. Whereas the coherences are damped and will fade away in the long-term limit, the equation governing the dynamics of populations in this case just becomes a simple rate equation with transition rates from b to a given by ${\gamma }_{ab,ab}$:

$$\begin{array}{c}\hfill {\dot{\rho }}_{aa}=\sum _b{\gamma }_{ab,ab}{\rho }_{bb}-\sum _b{\gamma }_{ba,ba}{\rho }_{aa}.\end{array}$$

(9)

Instead, for states with the same energies, the populations of the system density matrix are coupled to the coherences of the states with the same energy. The treatment that disregards all couplings of the populations to the coherences will in this paper be denoted the rotating wave approximation (RWA). In contrast, the treatment that preserves the couplings to the degenerate coherences will be denoted the secular approximation (BMS).

Second, for a single reservoir in thermal equilibrium (7), the correlation functions acquire additional analytic properties, so-called Kubo–Martin–Schwinger relations (KMS), which enable a thermodynamically-consistent description even for degenerate systems. In the absence of chemical potentials, the KMS relations read:

$$\begin{array}{c}\hfill C_{\alpha \beta }(\tau )=C_{\beta \alpha }(-\tau -\mathrm{i}\beta ),\end{array}$$

(10)

and transfer to the even Fourier transforms as ${\gamma }_{\alpha \beta }(+\omega )={\gamma }_{\beta \alpha }(-\omega )e^{+\beta \omega }$. In the master equation, these eventually lead to a detailed balance, and the system thermalizes in the long run with the temperature of the reservoir, i.e., ${\overline{\rho }}_S\propto e^{-\beta H_S}$ is a stationary state of the master equation [15]. With a chemical potential and an interaction that conserves the total particle number, i.e., under the assumption that $[H_S,N_S]=[H_R,N_R]=[H_I,N_S+N_R]=0$, the KMS relations can be generalized to:

$$\begin{array}{c}\hfill \sum _{\overline{\alpha }}{A}_{\overline{\alpha }}{C}_{\alpha \overline{\alpha }}(\tau )=\sum _{\overline{\alpha }}{e}^{+\beta \mu N_S}{A}_{\overline{\alpha }}{e}^{-\beta \mu N_S}{C}_{\overline{\alpha }\alpha }(-\tau -\mathrm{i}\beta ),\end{array}$$

(11)

which is explicitly shown in Appendix A. This leads to the local detailed detailed balance (LDB) relation among the coefficients:

$$\begin{array}{c}\hfill \frac{{\gamma }_{ab,cd}}{{\gamma }_{dc,ba}}=e^{\beta \left[(E_b-E_a)-\mu (N_b-N_a)\right]},\end{array}$$

(12)

where $E_a$ and $N_a$ denote the energy and particle number of state a, respectively. Eventually, these relations imply equilibration of both the temperature and chemical potential [45], i.e., ${\overline{\rho }}_S=e^{-\beta (H_S-\mu N_S-{\varphi }_S)}$, where ${\varphi }_S=-{\beta }^{-1}\ln \mathrm{Tr}\left\{\exp -\beta (H_S-\mu N_S)\right\}$, is one stationary state of the BMS master equation, even in presence of degeneracies.

We extend our setup by admitting two kinds of drivings.

First, we may allow for a slow external driving of the system Hamiltonian $H_S\to H_S(t)$. However, this driving must be significantly slower than the decay time of the reservoir correlation functions. Furthermore, the driving should not lift the degeneracy in the energy spectrum, and the non-degenerate states should not cross at any time. In other words, the driving should only operate on well-separated eigenenergies. Under these assumptions, our approximations remain applicable, and we still arrive at the same microscopically-derived master equation. The only difference is that the previously constant Hamiltonian and all associated quantities become time-dependent: $H_S\to H_S(t)$, $E_a\to E_a(t)$, $L_{ab}\to L_{ab}(t)$, and $\left|a\right.〉\to \left|a(t)\right.〉$ in Equations (3) and (4).

Second, we can consider N multiple reservoirs $H_R={\sum }_{\nu }{H}_R^{(\nu )}$ held at different equilibrium states:

$$\begin{array}{c}\hfill {\overline{\rho }}_R=\underset{\nu }{⨂}{e}^{-{\beta }_{\nu }(H_R^{(\nu )}-{\mu }_{\nu }{N}_R^{(\nu )}-{\varphi }_R^{(\nu )})}.\end{array}$$

(13)

where we introduced the inverse temperatures ${\beta }_{\nu }$, chemical potentials ${\mu }_{\nu }$, particle number operators $N_R^{(\nu )}$ and thermodynamic grand-potentials ${\varphi }_R^{(\nu )}=-{\beta }_{\nu }^{-1}\ln \mathrm{Tr}\left\{\exp e^{-{\beta }_{\nu }(H_R^{(\nu )}-{\mu }_{\nu }{N}_R^{(\nu )})}\right\}$ of reservoir $\nu =1,\cdots ,N$, N denoting the total number of reservoirs. This directly (or after suitable transformations) often implies that the Lindblad generator can be additively decomposed in the reservoir index ν. The master equation in the presence of slow driving and multiple reservoirs can thus be formally written as:

$$\begin{array}{c}\hfill {\dot{\rho }}_S=\mathcal{L}(t){\rho }_S(t)={\mathcal{L}}_0(t){\rho }_S(t)+\sum _{\nu }{\mathcal{L}}^{(\nu )}(t){\rho }_S(t),\end{array}$$

(14)

where ${\mathcal{L}}_0(t){\rho }_S(t)\widehat{=}-\mathrm{i}\left[H_S(t),{\rho }_S(t)\right]$ describes the action of the driven system Hamiltonian only. As discussed before for a single reservoir, the dissipator associated with reservoir ν will obey detailed balance relations leading to:

$$\begin{array}{c}\hfill {\mathcal{L}}^{(\nu )}(t){\rho }_{\mathrm{eq}}^{(\nu )}(t)=0,\end{array}$$

(15)

where we have introduced the time-dependent grand-canonical equilibrium state:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{eq}}^{(\nu )}(t)=e^{-{\beta }_{\nu }[H_S(t)-{\mu }_{\nu }{N}_S-{\varphi }_S^{(\nu )}(t)]}\end{array}$$

(16)

in terms of the system Hamiltonian $H_S(t)$, system particle number operator $N_S$ and the system grand-potentials ${\varphi }_S^{(\nu )}(t)=-{\beta }_{\nu }^{-1}\ln \mathrm{Tr}\left\{e^{-{\beta }_{\nu }[H_S(t)-{\mu }_{\nu }{N}_S]}\right\}$.

2.2. Average Thermodynamics

The change of the system energy under the quantum master equation dynamics can be decomposed as:

$$\begin{array}{ccc}\dot{E}\hfill & =\hfill & \frac{d}{dt}\mathrm{Tr}\left\{H_S(t)\rho (t)\right\}\hfill \\ & =\hfill & \mathrm{Tr}\left\{{\dot{H}}_S\rho \right\}+{\sum }_{\nu }{\mu }_{\nu }\mathrm{Tr}\left\{N_S{\mathcal{L}}^{(\nu )}\rho \right\}+{\sum }_{\nu }\mathrm{Tr}\left\{(H_S(t)-{\mu }_{\nu }{N}_S){\mathcal{L}}^{(\nu )}\rho \right\}\hfill \\ & =\hfill & \dot{W}+{\sum }_{\nu }{\dot{Q}}^{(\nu )},\hfill \end{array}$$

(17)

where we omit the system index S on the density matrix and the time dependence in the Liouvillians ${\mathcal{L}}^{(\nu )}$ for brevity. The work performed on the system contains a mechanical contribution (${\dot{W}}_m$) due to the external driving and a chemical one (${\dot{W}}_c$) due to the particle transfers with the reservoirs, $\dot{W}={\dot{W}}_m+{\dot{W}}_c$, where:

$$\begin{array}{c}\hfill {\dot{W}}_m=\mathrm{Tr}\left\{{\dot{H}}_S\rho \right\},\mathrm{and}{\dot{W}}_c={\sum }_{\nu }{\mu }_{\nu }\mathrm{Tr}\left\{N_S{\mathcal{L}}^{(\nu )}\rho \right\}.\end{array}$$

(18)

The heat current entering the system from reservoir ν is:

$$\begin{array}{c}\hfill {\dot{Q}}^{(\nu )}=\mathrm{Tr}\left\{(H_S(t)-{\mu }_{\nu }{N}_S){\mathcal{L}}^{(\nu )}\rho \right\}.\end{array}$$

(19)

After having established the first law, we now turn to the second law and introduce the von Neumann entropy, which represents the system entropy:

$$\begin{array}{c}\hfill S(t)=-\mathrm{Tr}\left\{\rho \ln \rho \right\}.\end{array}$$

(20)

Its time evolution is given by:

$$\begin{array}{c}\hfill \dot{S}=-\frac{d}{dt}\mathrm{Tr}\left\{\rho \ln \rho \right\}=-\mathrm{Tr}\left\{\dot{\rho }\ln \rho \right\},\end{array}$$

(21)

where we used $\mathrm{Tr}\left\{\rho \frac{d}{dt}\ln \rho \right\}=0$ (this can be shown using the fact that the density matrix can be diagonalized by unitary transformation). Using the Lindblad generator, one can directly see that the Hamiltonian driving does not directly contribute to the change of entropy, so that we have:

$$\begin{array}{c}\hfill \dot{S}=-\sum _{\nu }\mathrm{Tr}\left\{\left[{\mathcal{L}}^{(\nu )}\rho \right]\ln \rho \right\}.\end{array}$$

(22)

The entropy production is then given by the sum of the system entropy change plus the entropy change in the reservoirs (caused by the heat flows):

$$\begin{array}{c}\hfill {\dot{S}}_{\mathrm{i}}\equiv \dot{S}-\sum _{\nu }{\beta }_{\nu }{\dot{Q}}^{(\nu )}\ge 0.\end{array}$$

(23)

This expression can be proven to be positive by using Spohn’s inequality [46], but we also provide a direct proof in Appendix B.

We finish with a note on the Shannon entropy of the system, which by construction depends on the basis $S_{\mathrm{Sh}}=-{\sum }_i{\rho }_{ii}\ln {\rho }_{ii}$. For master equations in the rotating wave approximation, the basis chosen is the energy eigenbasis. This Shannon entropy does not depend on the eigenstate coherences, which anyway evolve independently of the populations. Furthermore, it is larger than or equal to the von Neumann entropy. Indeed, the relative entropy between a density matrix and its diagonal part ${\rho }_D$ reads:

$$\begin{array}{c}\hfill D(\rho ,{\rho }_D)=\mathrm{Tr}\left\{\rho \ln \rho -\rho \ln {\rho }_D\right\}=-S+S_{\mathrm{Sh}}-\mathrm{Tr}\left\{(\rho -{\rho }_D)\ln {\rho }_D\right\}.\end{array}$$

(24)

Since $(\rho -{\rho }_D)$ only contains off-diagonal matrix elements, whereas $\ln {\rho }_D$ has only entries on the diagonal, we have that:

$$\begin{array}{c}\hfill \mathrm{Tr}\left\{(\rho -{\rho }_D)\ln {\rho }_D\right\}=\sum _{ij}{(\rho -{\rho }_D)}_{ij}{\left(\ln {\rho }_D\right)}_{ji}=\sum _i{(\rho -{\rho }_D)}_{ii}{\left(\ln {\rho }_D\right)}_{ii}=0.\end{array}$$

(25)

Since the relative entropy is non-negative $D(\rho ,{\rho }_D)\ge 0$ under dynamics generated by a Lindblad master equation as we consider here, it follows that $S\le S_{\mathrm{Sh}}$. Note however that $\dot{S}$ and ${\dot{S}}_{\mathrm{Sh}}$ do not obey a general inequality. Similarly, the correct entropy production rate (23) and a Shannon-based entropy production rate ${\dot{S}}_{\mathrm{i}}^{\mathrm{Sh}}={\dot{S}}_{\mathrm{Sh}}-{\sum }_{\nu }{\beta }_{\nu }{\dot{Q}}^{(\nu )}$ are not generally related by an inequality.

2.3. Fluctuating Thermodynamics

2.3.1. Counting Statistics

Within the same approximations used to derive the quantum master equation, one can derive the full counting statistics for the energy and matter transfers using the dressed master equation formalism [9,38]. The measurement scheme corresponds to two point projective measurements of the energy $H_R^{(\nu )}$ and particle number $N_R^{(\nu )}$ in the reservoirs $\nu =1,\cdots ,N$. The energy and particle transfer generating function $G(\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)$ is then obtained by taking the trace of the dressed density matrix of the system $\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)$:

$$\begin{array}{c}\hfill G(\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)=\mathrm{Tr}\left\{\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\right\},\end{array}$$

(26)

where the counting field vectors $\{{\xi }_{\nu }\}=\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_N\}$ and $\{{\lambda }_{\nu }\}=\{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_N\}$ account for, respectively, the energy and matter currents out of the reservoirs. The dressed system density matrix satisfies the dressed quantum master equation:

$$\begin{array}{ccc}\dot{\rho }(\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\hfill & =\hfill & \mathcal{L}(\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\hfill \\ & \equiv \hfill & -\mathrm{i}\left[H_S(t)+{\sum }_{ab}{\sum }_{\tilde{\nu }}{\sigma }_{ab}^{(\tilde{\nu })}(t)L_{ab}(t),\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\right]+{\sum }_{ab,cd}{\sum }_{\tilde{\nu }}{\gamma }_{ab,cd}^{(\tilde{\nu })}(t)\times \hfill \\ & & \left[C_{ab,cd}({\xi }_{\tilde{\nu }},{\lambda }_{\tilde{\nu }},t)L_{ab}(t)\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)L_{cd}^{†}(t)-\frac{1}{2}\left\{L_{cd}^{†}(t)L_{ab}(t),\rho (\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)\right\}\right],\hfill \end{array}$$

(27)

whose dressed Liouvillian depends on the counting fields. The factors:

$$\begin{array}{c}\hfill C_{ab,cd}({\xi }_{\nu },{\lambda }_{\nu },t)=\exp \left\{\left[\mathrm{i}{\xi }_{\nu }(E_b(t)-E_a(t))+\mathrm{i}{\lambda }_{\nu }(N_b-N_a)\right]\right\},\end{array}$$

(28)

contain the counting fields keeping track of the energy and matter transfers with the reservoirs. The dressed quantum master Equation (3) reduces to the regular quantum master equation for the system reduced density matrix when the counting fields are set equal to zero, i.e., $\{{\xi }_{\nu }\}=\{{\lambda }_{\nu }\}=\{0\}$.

The joined distribution for the energy and matter currents out of the reservoirs $P(\{J_E^{(\nu )}\},\{J_M^{(\nu )}\},t)$ is obtained by using the Fourier transform:

$$\begin{array}{c}\hfill P(\{J_E^{(\nu )}\},\{J_M^{(\nu )}\},t)={\int }_{-\infty }^{\infty }\left[\prod _{\nu }t\frac{d{\xi }_{\nu }}{2\pi }\right]{\int }_0^{2\pi }\left[\prod _{\nu }t\frac{d{\lambda }_{\nu }}{2\pi }\right]e^{\mathrm{i}{\sum }_{\nu }\left({\xi }_{\nu }\Delta E_{\nu }+{\lambda }_{\nu }\Delta N_{\nu }\right)}G(\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t),\end{array}$$

(29)

where $\Delta E_{\nu }$ and $\Delta N_{\nu }$ are the energy and particle number changes in reservoir ν over a duration t and $J_E^{(\nu )}=\Delta E_{\nu }/t$ and $J_M^{(\nu )}=\Delta N_{\nu }/t$ denote the corresponding energy and matter currents, respectively.

To calculate the counting statistics of the mechanical work, a projective measurement in the system Hamiltonian $H_S(t)$ is required. The generating function for the associated counting statistics can be written as: [9]

$$\begin{array}{c}\hfill G(\alpha ,t)=\mathrm{Tr}\left\{e^{\mathrm{i}\alpha H_S(t)}\left(\mathcal{T}\exp {\int }_0^{t}d\tau \mathcal{L}(\tau )\right)\left(e^{-\mathrm{i}\alpha H_S(0)}\rho (0)\right)\right\},\end{array}$$

(30)

where the counting field α counts the energy changes in the system.

Since the mechanical work is the system energy change minus the total energy that has flown to the reservoirs, the generating function for mechanical power and energy and matter currents can be written as [47]:

$$\begin{array}{c}\hfill G(\alpha ,\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)=\mathrm{Tr}\left\{e^{\mathrm{i}\alpha H_S(t)}\left(\mathcal{T}\exp {\int }_0^{t}d\tau \mathcal{L}(\{{\xi }_{\nu }-\alpha \},\{{\lambda }_{\nu }\},\tau )\right)\left(e^{-\mathrm{i}\alpha H_S(0)}\rho (0)\right)\right\},\end{array}$$

(31)

where α is now the mechanical work counting field. Furthermore, $\mathcal{T}\exp \left\{·\right\}$ denotes the time-ordered exponential and $\rho (0)$ the initial density matrix of the system. By Fourier transform, we get the corresponding probability distribution:

$$\begin{array}{ccc}& & P(w,\{J_E^{(\nu )}\},\{J_M^{(\nu )}\},t)=\hfill \\ & & {\int }_{-\infty }^{\infty }\frac{d\alpha }{2\pi }{\int }_{-\infty }^{\infty }\left[\prod _{\nu }t\frac{d{\xi }_{\nu }}{2\pi }\right]{\int }_0^{2\pi }\left[\prod _{\nu }t\frac{d{\lambda }_{\nu }}{2\pi }\right]e^{\mathrm{i}\alpha w+\mathrm{i}{\sum }_{\nu }\left({\xi }_{\nu }\Delta E^{(\nu )}+{\lambda }_{\nu }\Delta N^{(\nu )}\right)}G(\alpha ,\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t),\hfill \end{array}$$

(32)

where w denotes the mechanical work performed on the system over time t.

2.3.2. Finite-Time Fluctuation Theorem

We now consider the generating function (31) when the system is driven by a time-dependent protocol, $H_S(\tau )$ for $\tau \in \left[0,t\right]$, and initially at equilibrium with reservoir $\nu =1$:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{eq}}^{(1)}(0)=e^{-{\beta }_1[H_S(0)-{\mu }_1{N}_S-{\varphi }_S^{(1)}(0)]}.\end{array}$$

(33)

We also consider the corresponding backward process where the system is driven by the time-reversed protocol, ${\tilde{H}}_S(\tau )=H_S(t-\tau )$ for $\tau \in \left[0,t\right]$ and initially at equilibrium with reservoir $\nu =1$ at the final time of the forward protocol:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{eq}}^{(1)}(t)=e^{-{\beta }_1[H_S(t)-{\mu }_1{N}_S-{\varphi }_S^{(1)}(t)]}.\end{array}$$

(34)

Since the Liouvillian depends parametrically on time through the system Hamiltonian, the generating function for the backward process is given by:

$$\begin{array}{c}\hfill \tilde{G}(\alpha ,\{{\xi }_{\nu }\},\{{\lambda }_{\nu }\},t)=\mathrm{Tr}\left\{e^{\mathrm{i}\alpha {\tilde{H}}_S(t)}\left(\mathcal{T}\exp {\int }_0^{t}d\tau \tilde{\mathcal{L}}(\{{\xi }_{\nu }-\alpha \},\{{\lambda }_{\nu }\},\tau )\right)\left(e^{-\mathrm{i}\alpha {\tilde{H}}_S(0)}{\rho }_{\mathrm{eq}}^{(1)}(t)\right)\right\},\end{array}$$

(35)

where $\tilde{\mathcal{L}}(\{{\xi }_{\nu }-\alpha \},\{{\lambda }_{\nu }\},\tau )=\mathcal{L}(\{{\xi }_{\nu }-\alpha \},\{{\lambda }_{\nu }\},t-\tau )$. In the following, we take reservoir $\nu =1$ as a reference for the energy and particle number counting. Accordingly, we set ${\xi }_1={\lambda }_1=0$ and introduce the new counting field vectors ${\{{\xi }_{\nu }\}}^{\prime }=\{{\xi }_2,{\xi }_3,\cdots ,{\xi }_N\}$ and ${\{{\lambda }_{\nu }\}}^{\prime }=\{{\lambda }_2,{\lambda }_3,\cdots ,{\lambda }_N\}$. Using the LDB relation (12), we find the symmetry relation:

$$\begin{array}{c}\hfill {\mathcal{L}}^{†}({\{{\xi }_{\nu }-\alpha \}}^{\prime },{\{{\lambda }_{\nu }\}}^{\prime },\tau )=e^{-{\beta }_1{\mu }_1{N}_S}\mathcal{L}({\{\mathrm{i}{A}_{\nu }^{ϵ}+{\xi }_{\nu }-(\mathrm{i}{\beta }_1+\alpha )\}}^{\prime },{\{-\mathrm{i}{A}_{\nu }^n+{\lambda }_{\nu }\}}^{\prime },\tau )e^{+{\beta }_1{\mu }_1{N}_S}\end{array}$$

(36)

expressed in terms of the thermodynamic affinities:

$$\begin{array}{c}\hfill A_{\nu }^E={\beta }_1-{\beta }_{\nu },A_{\nu }^N=-{\beta }_1{\mu }_1+{\beta }_{\nu }{\mu }_{\nu },\end{array}$$

(37)

and where ${\mathcal{L}}^{†}$ denotes the conjugate transpose in the system Liouville space, that is $\mathrm{Tr}\left\{|a_1\rangle \langle a_2|{\mathcal{L}}^{†}\left(|a_3\rangle \langle a_4|\right)\right\}=\mathrm{Tr}\left\{|a_3\rangle \langle a_4|\mathcal{L}\left(|a_1\rangle \langle a_2|\right)\right\}$, where $|a_i\rangle$ for $i=1,\cdots ,4$ are arbitrary quantum states in the system Hilbert space.

This symmetry (36) combined with the initial conditions (33) and (34) implies the finite-time fluctuation theorem:

$$\begin{array}{ccc}G(\alpha ,{\{{\xi }_{\nu }\}}^{\prime },{\{{\lambda }_{\nu }\}}^{\prime },t)\hfill & =\hfill & \mathrm{Tr}\left\{e^{\mathrm{i}\alpha H_S(t)}\left(\mathcal{T}{e}^{{\int }_0^t\tilde{\mathcal{L}}({\{{\xi }_{\nu }-\alpha \}}^{\prime },{\{{\lambda }_{\nu }\}}^{\prime },t-\tau )d\tau }\right)e^{-\mathrm{i}\alpha H_S(0)}{\rho }_{\mathrm{eq}}^{(1)}(0)\right\}\hfill \\ & =\hfill & \mathrm{Tr}{\left\{e^{\mathrm{i}\alpha H_S(0)}{\rho }_{\mathrm{eq}}^{(1)}(0)\left(\mathcal{T}{e}^{{\int }_0^t{\tilde{\mathcal{L}}}^{†}({\{{\xi }_{\nu }-\alpha \}}^{\prime },{\{{\lambda }_{\nu }\}}^{\prime },\tau )d\tau }\right)e^{-\mathrm{i}\alpha H_S(t)}\right\}}^{*}\hfill \\ & =\hfill & \mathrm{Tr}\left\{e^{\mathrm{i}(\alpha +\mathrm{i}{\beta }_1)H_S(0)}\left(\mathcal{T}{e}^{{\int }_0^t{\mathcal{L}}^{†}({\{\mathrm{i}{A}_{\nu }^E+{\xi }_{\nu }-(\mathrm{i}{\beta }_1+\alpha )\}}^{\prime },{\{-\mathrm{i}{A}_{\nu }^N+{\lambda }_{\nu }\}}^{\prime },\tau )d\tau }\right)\right.\times \hfill \\ & & {\left.e^{-\mathrm{i}(\alpha +\mathrm{i}{\beta }_1)H_S(t)}{\rho }_{\mathrm{eq}}^{(1)}(t)\right\}}^{*}{e}^{-{\beta }_1\Delta {\varphi }_1}\hfill \\ & =\hfill & \tilde{G}{(\mathrm{i}{\beta }_1+\alpha ,{\{-\mathrm{i}{A}_{\nu }^E+{\xi }_{\nu }\}}^{\prime },{\{-\mathrm{i}{A}_{\nu }^N+{\lambda }_{\nu }\}}^{\prime },t)}^{*}{e}^{-{\beta }_1\Delta {\varphi }_1}\hfill \\ & =\hfill & \tilde{G}(\mathrm{i}{\beta }_1-\alpha ,{\{-\mathrm{i}{A}_{\nu }^E-{\xi }_{\nu }\}}^{\prime },{\{-\mathrm{i}{A}_{\nu }^N-{\lambda }_{\nu }\}}^{\prime },t)e^{-{\beta }_1\Delta {\varphi }_1},\hfill \end{array}$$

(38)

where $\Delta {\varphi }_S^{(1)}={\varphi }_S^{(1)}(t)-{\varphi }_S^{(1)}(0)$. At the probability level, the finite time fluctuation theorem is given by:

$$\begin{array}{c}\hfill \ln \frac{P(+w,{\{+J_E^{(\nu )}\}}^{\prime },{\{+J_M^{(\nu )}\}}^{\prime },t)}{\tilde{P}(-w,{\{-J_E^{(\nu )}\}}^{\prime },{\{-J_M^{(\nu )}\}}^{\prime },t)}={\beta }_1(w-\Delta {\varphi }_S^{(1)})+\sum _{\nu =2}^N\left(A_{\nu }^E{J}_E^{(\nu )}+A_{\nu }^N{J}_M^{(\nu )}\right)t,\end{array}$$

(39)

where $P(+w,{\{+J_E^{(\nu )}\}}^{\prime },{\{+J_M^{(\nu )}\}}^{\prime },t)=\int d{J}_E^{(1)}\int d{J}_N^{(1)}P(w,\{J_E^{(\nu )}\},\{J_M^{(\nu )}\},t)$ and $\tilde{P}(-w,{\{-J_E^{(\nu )}\}}^{\prime },{\{-J_M^{(\nu )}\}}^{\prime },t)$ denotes the corresponding probability distribution along the backward process. This fluctuation theorem (39) holds for any given time t and is exclusively expressed in terms of the mechanical power and the energy and matter currents. It is the quantum analogue of the classical result derived in [34].

3. Degenerate Single Quantum Dot Circuit

We now illustrate our formalism by considering a specific model consisting of two degenerate quantum dots connected to two electron leads; see Figure 1. After defining the model, we first study its average thermodynamics. We then compare its counting statistics with and without eigenstate coherences and show that both satisfy the finite time fluctuation theorem derived above.

3.1. Model

We consider a double quantum dot with no direct tunneling between the dots, but exactly degenerate on-site energies. In general, it is well known that exact degeneracies may give rise to rich dynamics [36,48]. For our particular model, it is from a transport perspective also well known that negative differential conductance may arise from the Coulomb interaction due to coherences [35,37,38,39]. The effect has been observed experimentally [40] and is also present beyond the sequential tunneling regime [41]. A distinctive feature of this system is that the attached fermionic contacts allow for electron jumps into superposition states. The system, interaction and reservoir Hamiltonians read:

$$\begin{array}{ccc}{H}_S\hfill & =\hfill & ϵ\left(d_t^{†}{d}_t+d_b^{†}{d}_b\right)+U{d}_t^{†}{d}_t{d}_b^{†}{d}_b,\hfill \\ H_I\hfill & =\hfill & {\sum }_{\nu \in \{L,R\}}{\sum }_{i\in \{t,b\}}{\sum }_k\left[t_{k\nu i}{d}_i{c}_{k\nu }^{†}+t_{k\nu i}^{*}{c}_{k\nu }{d}_i^{†}\right],\hfill \\ H_R\hfill & =\hfill & {\sum }_{k\nu }{ϵ}_{k\nu }{c}_{k\nu }^{†}{c}_{k\nu }.\hfill \end{array}$$

(40)

Here, the on-site energies ϵ of the top (t) and bottom (b) dot are degenerate, and U denotes their Coulomb interaction. The $t_{k\nu ,i}$ denote the tunneling amplitudes into mode k of lead ν with energy ${ϵ}_{k\nu }$ from dot i (top or bottom). It is visible that both leads may trigger electronic jumps into both dots.

First, we remark that for charged states, not all superposition states are allowed. In particular, we cannot form superpositions of differently-charged states, such that coherences between, e.g., the empty and doubly-occupied states can be neglected from the beginning. Formally, they will evolve in a decoupled (and damped) fashion, but in reality, they cannot be created in a system-local state and will therefore vanish throughout. Denoting the diagonal matrix elements of the empty, the top occupied, the bottom occupied and the doubly-occupied state by ${\rho }_0$, ${\rho }_t$, ${\rho }_b$, and ${\rho }_2$, respectively, and the admissible coherences between the singly-charged states by ${\rho }_{tb}$ and ${\rho }_{bt}={\rho }_{tb}^{*}$, the BMS Lindblad master Equation (3) becomes (see Appendix C for more details on the derivation):

$$\begin{array}{ccc}{\dot{\rho }}_0\hfill & =\hfill & -\left[({\mathrm{\Gamma }}_{Lt}+{\mathrm{\Gamma }}_{Lb})f_L+({\mathrm{\Gamma }}_{Rt}+{\mathrm{\Gamma }}_{Rb})f_R\right]{\rho }_0\hfill \\ & & +\left[{\mathrm{\Gamma }}_{Lt}(1-f_L)+{\mathrm{\Gamma }}_{Rt}(1-f_R)\right]{\rho }_t+\left[{\mathrm{\Gamma }}_{Lb}(1-f_L)+{\mathrm{\Gamma }}_{Rb}(1-f_R)\right]{\rho }_b\hfill \\ & & +\left[{\gamma }_L(1-f_L)+{\gamma }_R(1-f_R)\right]{\rho }_{tb}+\left[{\gamma }_L^{*}(1-f_L)+{\gamma }_R^{*}(1-f_R)\right]{\rho }_{bt},\hfill \\ {\dot{\rho }}_t\hfill & =\hfill & -\left[{\mathrm{\Gamma }}_{Lt}(1-f_L)+{\mathrm{\Gamma }}_{Rt}(1-f_R)+{\mathrm{\Gamma }}_{Lb}{f}_L^U+{\mathrm{\Gamma }}_{Rb}{f}_R^U\right]{\rho }_t\hfill \\ & & +\left[{\mathrm{\Gamma }}_{Lt}{f}_L+{\mathrm{\Gamma }}_{Rt}{f}_R\right]{\rho }_0+\left[{\mathrm{\Gamma }}_{Lb}(1-f_L^U)+{\mathrm{\Gamma }}_{Rb}(1-f_R^U)\right]{\rho }_2\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L(f_L^U-(1-f_L))+{\gamma }_R(f_R^U-(1-f_R))-\mathrm{i}{\gamma }_L{\Sigma }_L-\mathrm{i}{\gamma }_R{\Sigma }_R\right]{\rho }_{tb}\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L^{*}(f_L^U-(1-f_L))+{\gamma }_R^{*}(f_R^U-(1-f_R))+\mathrm{i}{\gamma }_L^{*}{\Sigma }_L+\mathrm{i}{\gamma }_R^{*}{\Sigma }_R\right]{\rho }_{bt},\hfill \\ {\dot{\rho }}_b\hfill & =\hfill & -\left[{\mathrm{\Gamma }}_{Lb}(1-f_L)+{\mathrm{\Gamma }}_{Rb}(1-f_R)+{\mathrm{\Gamma }}_{Lt}{f}_L^U+{\mathrm{\Gamma }}_{Rt}{f}_R^U\right]{\rho }_b\hfill \\ & & +\left[{\mathrm{\Gamma }}_{Lb}{f}_L+{\mathrm{\Gamma }}_{Rb}{f}_R\right]{\rho }_0+\left[{\mathrm{\Gamma }}_{Lt}(1-f_L^U)+{\mathrm{\Gamma }}_{Rt}(1-f_R^U)\right]{\rho }_2\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L(f_L^U-(1-f_L))+{\gamma }_R(f_R^U-(1-f_R))+\mathrm{i}{\gamma }_L{\Sigma }_L+\mathrm{i}{\gamma }_R{\Sigma }_R\right]{\rho }_{tb}\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L^{*}(f_L^U-(1-f_L))+{\gamma }_R^{*}(f_R^U-(1-f_R))-\mathrm{i}{\gamma }_L^{*}{\Sigma }_L-\mathrm{i}{\gamma }_R^{*}{\Sigma }_R\right]{\rho }_{bt},\hfill \\ {\dot{\rho }}_2\hfill & =\hfill & -\left[({\mathrm{\Gamma }}_{Lt}+{\mathrm{\Gamma }}_{Lb})(1-f_L^U)+({\mathrm{\Gamma }}_{Rt}+{\mathrm{\Gamma }}_{Rb})(1-f_R^U)\right]{\rho }_2\hfill \\ & & +\left[{\mathrm{\Gamma }}_{Lb}{f}_L^U+{\mathrm{\Gamma }}_{Rb}{f}_R^U\right]{\rho }_t+\left[{\mathrm{\Gamma }}_{Lt}{f}_L^U+{\mathrm{\Gamma }}_{Rt}{f}_R^U\right]{\rho }_b\hfill \\ & & -\left[{\gamma }_L{f}_L^U+{\gamma }_R{f}_R^U\right]{\rho }_{tb}-\left[{\gamma }_L^{*}{f}_L^U+{\gamma }_R^{*}{f}_R^U\right]{\rho }_{bt},\hfill \\ {\dot{\rho }}_{tb}\hfill & =\hfill & -\frac{1}{2}[({\mathrm{\Gamma }}_{Lt}+{\mathrm{\Gamma }}_{Lb})(f_L^U+(1-f_L))+({\mathrm{\Gamma }}_{Rt}+{\mathrm{\Gamma }}_{Rb})(f_R^U+(1-f_R))\hfill \\ & & -\mathrm{i}{\mathrm{\Gamma }}_{Lt}{\Sigma }_L+\mathrm{i}{\mathrm{\Gamma }}_{Lb}{\Sigma }_L-\mathrm{i}{\mathrm{\Gamma }}_{Rt}{\Sigma }_R+\mathrm{i}{\mathrm{\Gamma }}_{Rb}{\Sigma }_R]{\rho }_{tb}\hfill \\ & & +\left[{\gamma }_L^{*}{f}_L+{\gamma }_R^{*}{f}_R\right]{\rho }_0-\left[{\gamma }_L^{*}(1-f_L^U)+{\gamma }_R^{*}(1-f_R^U)\right]{\rho }_2\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L^{*}(f_L^U-(1-f_L))+{\gamma }_R^{*}(f_R^U-(1-f_R))-\mathrm{i}{\gamma }_L^{*}{\Sigma }_L-\mathrm{i}{\gamma }_R^{*}{\Sigma }_R\right]{\rho }_t\hfill \\ & & +\frac{1}{2}\left[{\gamma }_L^{*}(f_L^U-(1-f_L))+{\gamma }_R^{*}(f_R^U-(1-f_R))+\mathrm{i}{\gamma }_L^{*}{\Sigma }_L+\mathrm{i}{\gamma }_R^{*}{\Sigma }_R\right]{\rho }_b,\hfill \\ {\dot{\rho }}_{bt}\hfill & =\hfill & {\dot{\rho }}_{tb}^{*},\hfill \end{array}$$

(41)

where we used the wide-band limit for the tunneling rates:

$$\begin{array}{ccc}{\mathrm{\Gamma }}_{\nu i}(\omega )\hfill & =\hfill & 2\pi {\sum }_k{\left|t_{k\nu ,i}\right|}^2\delta (\omega -{ϵ}_{k\nu })\to {\mathrm{\Gamma }}_{\nu i},\hfill \\ {\gamma }_{\nu }(\omega )\hfill & =\hfill & 2\pi {\sum }_k{t}_{k\nu ,t}{t}_{k\nu ,b}^{*}\delta (\omega -{ϵ}_{k\nu })\to {\gamma }_{\nu }.\hfill \end{array}$$

(42)

Whereas the tunneling rates ${\mathrm{\Gamma }}_{\nu i}$ are rates in the traditional sense ${\mathrm{\Gamma }}_{\nu i}\ge 0$ and describe tunneling processes into top- and bottom-localized electronic states, respectively, this is different for the unconventional complex-valued rates ${\gamma }_{\nu }$. Formally, we see that the ${\gamma }_{\nu }$ mediate the coupling between coherences and populations and, thus, allow the system to jump, e.g., from the empty state into a superposition of the singly-charged states. Depending on the microscopic details of the coupling, the phases of the tunneling amplitudes in Equation (42) may interfere destructively (such that ${\gamma }_{\nu }\to 0$, which is equivalent to taking the RWA limit) or constructively (when all tunneling amplitudes are equal, we have ${\left|{\gamma }_{\nu }\right|}^2={\mathrm{\Gamma }}_{\nu t}{\mathrm{\Gamma }}_{\nu b}$). This last limit of constructive interference ${\gamma }_{\nu }\to \sqrt{{\mathrm{\Gamma }}_{\nu t}{\mathrm{\Gamma }}_{\nu b}}$ will be used here as the wide-band limit of the secular approximation.

The thermal reservoir properties are contained in the Fermi functions and Lamb-shift terms:

$$\begin{array}{ccc}{f}_{\nu }\hfill & =\hfill & \frac{1}{e^{{\beta }_{\nu }(ϵ-{\mu }_{\nu })}+1},f_{\nu }^U=\frac{1}{e^{{\beta }_{\nu }(ϵ+U-{\mu }_{\nu })}+1},\hfill \\ {\Sigma }_{\nu }\hfill & =\hfill & \frac{1}{\pi }\mathsf{\Re }\left[\mathrm{\Psi }\left(\frac{1}{2}+\mathrm{i}\frac{{\beta }_{\nu }(ϵ+U-{\mu }_{\nu })}{2\pi }\right)-\mathrm{\Psi }\left(\frac{1}{2}+\mathrm{i}\frac{{\beta }_{\nu }(ϵ-{\mu }_{\nu })}{2\pi }\right)\right],\hfill \end{array}$$

(43)

where $\mathrm{\Psi }(x)$ denotes the digamma function.

We stress a few things before proceeding. First, as the master equation is of the Lindblad form by construction, the density matrix properties will be preserved. Second, we see that the dissipator is additive in the reservoirs $\mathcal{L}={\mathcal{L}}_L+{\mathcal{L}}_R$. Each dissipator annihilates its associated Gibbs state, cf. Equation (15). Consequently, at global equilibrium (${\beta }_L={\beta }_R$ and ${\mu }_L={\mu }_R$), the thermal Gibbs state (with vanishing coherences) is the stationary state. Finite coherences in the steady state can however arise in nonequilibrium setups, as will be discussed below. Finally, we mention that the total Liouvillian becomes bistable when ${\mathrm{\Gamma }}_{Lt}={\mathrm{\Gamma }}_{Lb}={\mathrm{\Gamma }}_L$ and ${\mathrm{\Gamma }}_{Rt}={\mathrm{\Gamma }}_{Rb}={\mathrm{\Gamma }}_R$, and we will in the following avoid this situation. The graph of the master equation is depicted in Figure 1b.

3.2. Model Thermodynamics

In what follows, we will consider mainly for simplicity the limit ${\mathrm{\Gamma }}_{Lt}={\mathrm{\Gamma }}_{Rb}={\mathrm{\Gamma }}_A$ and ${\mathrm{\Gamma }}_{Lb}={\mathrm{\Gamma }}_{Rt}={\mathrm{\Gamma }}_B$ and $\gamma =\sqrt{{\mathrm{\Gamma }}_A{\mathrm{\Gamma }}_B}$ (or, for the RWA limit, $\gamma =0$); see Figure 1a. We assume ${\mathrm{\Gamma }}_A\ne {\mathrm{\Gamma }}_B$, so that we will not consider the bistable situation in the present paper [38].

We can extract the time-dependent energy and matter currents into and from both reservoirs, e.g., from the full counting statistics methods, as discussed in Section 2.3.1. In Appendix C.2, we provide the required counting fields exemplarily for transitions triggered by the left junction. Alternatively, we may also use definitions analogous to the heat current (19) to calculate energy and matter currents. In Figure 2, we plot the time-dependent matter and energy currents for our model versus the potential difference.

Previous investigations of this particular model [37,38] have already revealed a significant suppression of the steady state matter current due to coherences. The suppression of the currents is linked to a pure nonequilibrium steady state arising at low temperatures ${\beta }_{L}U={\beta }_{R}U\gg 1$ when $\Delta \mu ={\mu }_L-{\mu }_R=\pm (2ϵ+U)$, which for our particular parameters can be understood analytically; see Appendix C.4. Here, we complete this picture by the time-dependent evolution and the time-dependent energy current. Most important, we note the striking difference between the steady-state currents of energy and matter currents of the BMS (solid and dashed black) and the RWA (dotted black) versions. Interestingly, the time-dependent average current obtained in the current suppression valley within the BMS approximation first approaches the steady-state value obtained within the RWA approximation, before ultimately converging to the steady-state value given within the BMS master equation (see bottom right inset of Figure 2).

Furthermore, from the time-dependent solution of the master Equation (3), we can compute the Shannon entropy (in the original energy eigenbasis $\{\left|0\right.〉,\left|t\right.〉,\left|b\right.〉,\left|2\right.〉\}$) and the von Neumann entropy (basis independent). Since the first neglects the coherences, these will obviously differ in regions where coherences are present; see Figure 3.

In particular, we can see that the steady-state von Neumann entropy vanishes when $\Delta \mu =\pm (2ϵ+U)$, whereas the Shannon entropy does not, which nicely illustrates that the system reaches a stationary pure state at this nonequilibrium configuration; cf. Appendix C.4. The time-dependent von Neumann entropy of the system as calculated within the BMS and RWA relates in a similar fashion as described above for the time-dependent average current out of the left reservoir (bottom left inset of Figure 3). For a bias voltage out of the current suppression valley, entropies computed within the BMS and RWA remain approximately equal at all times. Within the current suppression valley, the former mirrors the later at short times before converging to its own steady-state value (vanishing system entropy at the current suppression point $\Delta \mu =\pm (2ϵ+U)$).

From the difference between the change of the system entropy and the heat currents, we can obtain the entropy production rate, Equation (23), which we plot in Figure 4.

Beyond the evident sanity check that it is positive, we see that even at steady state, coherences between the degenerate states may survive in a nonequilibrium setup, which goes along with a suppression of the steady-state entropy production rate.

We now briefly consider how our model can operate as a thermoelectric device. We consider the situation in which a thermal gradient is applied between the reservoirs (${\beta }_L>{\beta }_R$) to drive a current against a chemical potential bias ($\Delta \mu ={\mu }_L-{\mu }_R>0$). Denoting the electronic and energy current entering the system from the left reservoir by $J_M$ and $J_E$ (dropping the reservoir index L), the thermoelectric efficiency of this process is defined as the ratio between the generated power $P=-J_M\Delta \mu$ and the heat extracted from the hot right reservoir $-(J_E-{\mu }_R{J}_M)$:

$$\begin{array}{c}\hfill \eta =\frac{J_M\Delta \mu }{J_E-{\mu }_R{J}_M}*\Theta (-J_M\Delta \mu ).\end{array}$$

(44)

The Heaviside function is introduced to indicate that this efficiency is only meaningful in regions of positive power. Positivity of the steady-state entropy production rate implies, as usual, that this efficiency is upper-bounded by the Carnot efficiency, $\eta \le 1-{\beta }_R/{\beta }_L$. A strong thermoelectric effect requires a large temperature gradient, which in our model reduces the impact of the coherences. In Figure 5, we observe numerically that the region of positive power is outside the region where quantum coherences suppress the current. To obtain a non-negligible power output, we have to consider parameter ranges where the coherences do not significantly modify the energetics. Consequently, the BMS and RWA results are qualitatively the same. In particular, the quantum efficiencies (${\gamma }_{\nu }=\sqrt{{\mathrm{\Gamma }}_{\nu t}{\mathrm{\Gamma }}_{\nu b}}$, solid black) and classical efficiencies (${\gamma }_{\nu }=0$, dashed black) are rather close, although the quantum efficiency is always smaller than the classical one. We have numerically observed this inequality also for other parameters.

3.3. Statistics and Fluctuation Theorem

The generating function of work and currents (31) can be evaluated numerically by solving the dressed quantum master Equation (27) for the specific model (40), namely Equation (C9). The corresponding joined probability distribution is then obtained by a Fourier transform.

As an illustration, we now consider the system introduced in Section 3.1 in the isothermal regime $\beta ={\beta }_L={\beta }_R$. The initial condition of the system is taken as the grand canonical equilibrium with respect to the right reservoir:

$$\begin{array}{c}\hfill \rho (0)=\exp \left[-\beta (H_S-{\mu }_R{N}_S-{\varphi }_S^R)\right],\end{array}$$

(45)

where the equilibrium grand-potential is ${\varphi }_S^R=-{\beta }^{-1}\ln \mathrm{Tr}\left\{\exp \left[-\beta (H_S-{\mu }_R{N}_S)\right]\right\}$ and where $N_S=d_t^{†}{d}_t+d_b^{†}{d}_b$ is the particle number operator in the system.

The distribution $P(\Delta n,\tau )$ of the particle changes in the left reservoir $\Delta n=J_M\tau$ during time τ is numerically evaluated for three different values of the measurement time. The results from the master equation in the secular approximation are compared to those obtained from the RWA master Equation (9), in which one neglects the influence of quantum coherences on the dynamics and current statistics. In Figure 6, we see that the distributions obtained in the former case (i.e., BMS) exhibit a bimodal behavior in the transient regime, which approaches a long-tail distribution for large times. This was observed in a wide range of parameters close to the current suppression point $\Delta {\mu }^{*}=\pm (2ϵ+U)$. Qualitatively, this can be well understood from the fact that the system is close to a bistable configuration, associated with a near-block form of the Liouvillian: whereas one block supports a finite steady-state current, the current associated with the other subspace (with a dark state) vanishes, and telegraph-type averaging over the two distributions yields the visible long-tail distribution [49,50]. The diagonal initial state (45) then also explains why the long-term distribution starts at $\Delta n=1$: since the dark state is a superposition of the two singly-charged states, at least a single jump event is required to create it. This effect is totally absent in the latter case (i.e., RWA), where the distribution has the usual bell-shape whose drift gives the finite average current at steady state. The BMS drift instead is, as expected, very small, close to the current suppression point when coherences are taken into account (see Section 2.2). The BMS distribution is thus non-trivial and converges to a distribution with a large tail. This example shows that not only average currents are affected by sustained coherences, but also their statistics.

For our choice of initial condition (45), the statistics of the current flowing out of the left reservoir must satisfy the fluctuation symmetries (38) and (39). In the present case, the fluctuation relation (39) reduces to the finite-time fluctuation theorem for the net number of particles transferred to the left reservoir ($\Delta n=J_M\tau$):

$$\ln \frac{P(+\Delta n,\tau )}{P(-\Delta n,\tau )}=\beta \Delta \mu \Delta n,$$

(46)

where $\Delta \mu ={\mu }_L-{\mu }_R$. Since we first numerically evaluate the current generating function $G(\lambda ,\tau )={\tau }^{-1}{\sum }_{\Delta n}P(\Delta n,\tau )e^{-i\lambda \Delta n}$ and due to the highly oscillating integrals involved in obtaining the distribution $P(\Delta n,\tau )$ at large particle number changes $\Delta n$, it is however simpler to test the equivalent fluctuation theorem symmetry (38), which here reduces to:

$$G(\lambda ,\tau )=G(\mathrm{i}\beta \Delta \mu -\lambda ,\tau ).$$

(47)

This symmetry is indeed verified by the generating functions of the distributions shown in Figure 6. We note that the fluctuation theorem is also satisfied by the generating function obtained within the RWA (Figure 7) even though the two statistics significantly differ. The fact that the statistics obtained within the RWA also satisfies a finite time (FT) directly results from the fact that transition rates of the stochastic master Equation (9) satisfy the LDB relation (12).

4. Summary

In the present paper, we established the nonequilibrium thermodynamics of open quantum systems exhibiting degeneracies and described by quantum master Equation (3). We established the first and and second law, as well as a finite-time fluctuation theorem solely expressed in terms of the mechanical work and the energy and particle counting statistics. Using a simple model with two degenerate quantum dots, we showed that eigenbasis coherences at steady state can generate non-trivial counting statistics, such as bi-modality and diverging second and higher cumulants. These findings will help to elucidate the role of coherences in stochastic thermodynamics. A remaining open issue is to be able to treat close-to-degenerate eigenstates within the quantum master equation formalism. This is particularly important to treat drivings, which can induce crossings between the system eigenenergies.