Coherence enhanced quantum-dot heat engine

We show that quantum coherence can enhance the performance of a continuous quantum heat engine in the Lindblad description. We investigate the steady-state solutions of the particle-exchanging quantum heat engine, composed of degenerate double quantum dots coupled to two heat baths in parallel, where quantum coherence may be induced due to interference between relaxation channels. We find that the engine power can be enhanced by the coherence in the nonlinear response regime, when the symmetry of coupling configurations between dots and two baths is broken. In the symmetric case, the coherence cannot be maintained in the steady state, except for the maximum interference degenerate case, where initial-condition-dependent multiple steady states appear with a dark state.

Introduction -Quantum thermodynamics is an emerging field in view of significant progress of technology which allows to scale down heat-energy converting devices to nanoscale where quantum effects cannot be neglected [1]. Examples of such quantum heat engines (QHE) include lasers, solar cells, photosynthetic organisms, etc where along with few-level quantum structure [2][3][4] a phenomenon of quantum coherence plays an important role [5][6][7]. In particular, coherence in system-bath interactions that originates from the interference may enhance the power [8,9] and efficiency at maximum power [10] of the laser and solar cell and is responsible for highly efficient energy transfer in photosynthetic systems [11]. These effects have been confirmed in the experimental studies of polymer solar cells [12]. The noiseinduced coherence is a different kind from internal coherence usually involved with a system Hamiltonian [13] which was recently demonstrated in the nitrogen vacancy-based microscopic QHE in diamond [14], and manifests as an improved efficiency in spectroscopic pump-probe measurements [15].
In much of the literature, the quantum coherence effect has been studied on continuously-working bosonic devices [8][9][10][11]. In this Letter, we focus on the fermionic QHE autonomously working without external work agents like a driving laser, made up of repulsively interacting double quantum dots with the same energy levels, coupled to fermionic baths in parallel, depicted in Fig. 1. In contrast to previous studies for such a system [16][17][18], we introduce a parameter for the strength of interference between relaxation channels, which plays a crucial role. We derive the condition for maintaining quantum coherence in the steady state and investigate the engine performance, not only in terms of tunneling coefficients between dots and baths, but also interference strengths.
We find that the power enhancement of the QHE can be achieved in the nonlinear response regime. When coupling configurations assigned to each bath are symmetric to each other, a quantum coherence initially induced by interference between relaxation channels would eventually disappear * slung@postech.ac.kr † dorfmank@lps.ecnu.edu.cn ‡ hgpark@kias.re.kr FIG. 1. A schematic illustration of the QHE, composed of two heat baths and a two-dot system. The dot energies, E 1 and E 2 (in this work, E 1 = E 2 ) are higher than chemical potentials, µ h and µ c . w a d± represents transfer rate of a particle between dot d and bath a, and φ a w a 1± w a 2± denotes the interference amplitude. Inset: A circuit analogy of resistors in parallel.
in the long-time (steady-state) limit. The exceptional case emerges for the degenerate energy level configuration at the maximum interference strength, where the dynamics is found to be localized, manifested as a mathematical singularity in the evolution operator evoking the so-called dark state [19], characterized by multiple steady states with finite quantum coherence depending on a given initial state. This singularity should be also found for more general settings with coherent dynamics originated from the energy-level degeneracy and parallel couplings, including a single bath case. Note that a spurious quantum coherence can be observed for a very long time (quasi-stationary state regime) near the maximum interference.
When the coupling configuration symmetry is broken in terms of either tunnelling coefficients or interference strengths, a genuine new steady state emerges with nonvanishing quantum coherence in general, producing an extra quantum current between two baths through dots in addition to the conventional classical current. This quantum current yields an extra contribution to the engine power, which can be positive in a specific parameter regime.
Model -We first derive the quantum master equation (QME) [20] for the density operatorρ S (t) of the fermionic QHE in the limit of weak coupling to hot (h) and cold (c) baths, where a temperature difference T h − T c > 0 and a potential bias µ c − µ h > 0 are applied. For simplicity, we assume a single energy level for each quantum dot with the degenerate energy levels E 1 = E 2 = E and infinitely large repulsion between particles in dots. The system then can be described using three two-particle eigenstates: |0 denotes empty dots, |1 and |2 stand for occupation of dot 1 and 2, respectively by a single particle. In addition, coherent hopping between dots is also forbidden and the only source of coherence is due to coupling to thermal baths. The interaction between system and bath a(= h , c) is given byĤ a SB = d,k g a db a † k |0 d| + h.c., whereb a † k is the operator creating a single particle with momentum k in bath a and g a d is the tunneling coefficient between dot d(= 1, 2) and bath a. Following the standard procedure of tracing out bath degrees of freedom, we obtain the QME which reads (1) where the system HamiltonianĤ S = E (|1 1| + |2 2|) and the Lindblad operators areL 1 = |1 0|,L 2 = |2 0|,L 3 =L † 1 , and L 4 =L † 2 . Note that we neglected the Lamb shift term [21]. The dissipation matrix Γ a is given by (2) where w a d± represents transfer rate of a particle between dot d and bath a; the subscript + (−) denotes the inflow (outflow) with respect to the dot, respectively. These rates are given by w a d+ = 2π|g a d | 2 N a and w a d− = 2π|g a d | 2 N a , where N a = N a (E) is the Fermi-Dirac distribution in bath a and N a = 1−N a (see the derivation in Sec. ?? of the Supplemental Material(SM) [21]).
The off-diagonal terms in Eq. (2) represent interference between particle transfer associated with different dots. The interference effect will be manifested in nonzero off-diagonal terms ofρ S , e.g., 1|ρ S |2 0. The coherence, however, can be destroyed by other environmental noises not captured in the interaction Hamiltonian, like the gate voltage noise [22] which may dephase the system. Instead of the secular approximation removing such coherence entirely in the dissipation [23], we introduce a phenomenological parameter φ a for the interference term in Eq. (2), assigned to each bath; |φ a | = 1 stands for permitting the full interference of relaxations with bath a, while φ a = 0 corresponds to no quantum effect of system-bath interactions. In earlier bosonic QHE models, this interference parameter is governed by the angle between dipole moments corresponding to two dots which ensures that |φ a | ≤ 1 [8]. For convenience, φ a is treated as a real number. Note that the second term in Eq. (1) is a standard form of quantum dynamical semigroup [20], which guarantees the positive and trace preserving dynamics since Γ a in Eq. (2) is the positivesemidefinite matrix for |φ a | ≤ 1.
To solve the QME, it is convenient to map the density operator to a vector: P = (ρ 00 , ρ 11 , ρ 22 , ρ 12 , ρ 21 , ρ 01 , ρ 02 , ρ 10 , ρ 20 ) T , where ρ i j = i|ρ S | j . The last four components vanish in the long-time limit because there is no dynamics producing the coherence between the empty and occupied states so that only dephasing is allowed, as seen in Sec. ?? of SM [21]. Thus, we write the corresponding Liouville equation as where L is a 5 × 5 matrix with the reduced vector P = (ρ 00 , ρ 11 , ρ 22 , ρ 12 , ρ 21 ) T . Introducing W d = a w a d+ , W d = a w a d− , Φ = a φ a w a 1+ w a 2+ , and Φ = a φ a w a 1− w a 2− , the L matrix then reads Steady-state solutions -From the steady-state condition, LP(∞) = 0, we find the relations as with the population conservation (ρ 00 + ρ 11 + ρ 22 = 1) and Note that the classical solution is recovered from Eq. (5), when the coherence term vanishes (ρ 12 (∞) = 0). This classical incoherent condition is determined by Eq (6) as which is obviously satisfied for the trivial case with Φ = Φ = 0 (or equivalently φ a = 0). Note that the equilibrium case (T h = T c and µ h = µ c ) also satisfies this incoherent condition In general, Eqs. (5) and (6) leads to a 2 × 2 matrix equation for ρ 11 and ρ 22 as with Unless the determinant |L ss | vanishes, the steady state solution is uniquely defined, which is given explicitly in Eq. (??) of SM [21].
We next consider a special case, where the coupling strength ratio of two dots with a bath is the same for both baths, i.e. g h 2 /g h 1 = g c 2 /g c 1 ≡ r, leading to w a 2± /w a 1± = r 2 . This may be a natural situation in real experiments and will be called the r-symmetric configuration [17]. We take r > 0 for simplicity. Assuming an additional symmetry for the coherence parameter as φ h = φ c ≡ φ, one can show W d /W d = Φ/Φ even in nonequilibrium (N h N c ), satisfying the incoherence condition in Eq. (7). However, at the maximum interference (|φ| = 1), the matrix L ss becomes singular with |L ss | = 0 and multiple steady-state solutions emerge, which will be discussed later. With the broken symmetry (φ h φ c ), the quantum coherence survives with a non-classical solution (ρ 12 (∞) 0). In a more general case with g h 2 /g h 1 g c 2 /g c 1 , the classical solution is still possible by adjusting φ h and φ c appropriately to satisfy the incoherent condition, but L ss cannot be singular.
Steady-state currents -A particle current J a d representing a particle net flow to dot d from bath a can be obtained from the QME in Eq. (1) or Eq. (3) as which represents time increment of the particle density of dot d due to bath a. In the steady state, we expect that the particle density increment should be balanced by two reservoirs such and the total current is given by J = d J d (∞). Transferring an electron from bath h to bath c, the electron gains an energy governed by the difference between underlying chemical potentials µ c − µ h , thus the output power of the QHE is simply given by P = (µ c − µ h )J. As the heat flux from the high-temperature reservoir is given aṡ Q h = (E − µ h )J, the engine efficiency does not vary with the particle current as η = P/Q h = (µ c − µ h )/(E − µ h ).
The particle current can be further separated into the classical and quantum part as where the classical current is defined by setting φ a = 0 in Eqs. (5) and (10) as . For a proper classical engine to generate positive power (J cl d > 0), we consider ∆N > 0 only.
The second term represents the quantum current J q d ≡ Ψ d ρ 12 (∞), induced by the coherence, and the explicit expressions for the quantum speed Ψ d and ρ 12 (∞) are given in Sec. ?? of SM [21]. Note that the quantum current for each dot can be both positive and negative, depending on the parameter values, as well as the total quantum current J q = d J q d (see Fig. ?? of of SM [21]).
As ρ 12 (∞) is also proportional to bias ∆N, the QHE can be viewed as an analogue of an electronic circuit with resistors Relative quantum conductances of (a) dot 1 and (b) dot 2, denoted as σ q 1 /σ cl 1 and σ q 2 /σ cl 2 , respectively in the (φ c , φ h ) plane. Here, we used N h = 0.2 and N c = 0.1, and the r-symmetric configuration with |g a 1 | 2 = 8π/(1 + r 2 ) and |g a 2 | 2 = 8πr 2 /(1 + r 2 ) at r = 4. Along the line of symmetry (purple), ρ 12 (∞) = 0, while Ψ d = 0 defines the black line. The quantum conductances vanish along both lines. Note that a back flow where the negative quantum current overmatches the positive classical current. R 1 and R 2 in parallel under the external potential bias (see the inset of Fig. 1 The conductance is also divided into the classical and quantum part as . The classical part σ cl d is always positive, while the quantum part can be either positive or negative. In Fig. 2, we plot the relative quantum conductance σ q d /σ cl d in the (φ c , φ h ) plane in the r-symmetric configuration. Near but off the symmetric line of φ h = φ c , we find the total quantum conductance σ q = d σ q d > 0, which means that the performance of the QHE can be enhanced beyond the classical limit in this parameter regime.
For small ∆N, we expand the relative quantum conductance as for the r-symmetric configuration (see Sec. ?? of SM [21] for detailed calculations). Interestingly, σ q d is always non-positive in the linear response regime (S 0 d ≤ 0), but may become positive due to S 1 d in the nonlinear regime as ∆N increases for is relatively stronger for dot 1 which has a weaker coupling with baths, as also seen in Fig. 2, which might be applicable to a filtering circuit.
Although ρ 12 (∞) becomes finite off the symmetric line (φ h φ c ), the quantum current may vanish again when Ψ d = 0 in Eq. (11), which is denoted by black lines in Fig. 2. This can happen by balancing the quantum contributions from the stochastic part and the interference part, which are represented by first two terms and the third term in the right-hand-side of Eq. (10), respectively. The quantum enhancement occurs only between two lines of Ψ d = 0 and ρ 12 (∞) = 0. For general cases outside of the r-symmetric configuration, these two lines are simply tilted (see Fig. ?? in SM [21]), but the general features of the QHE are essentially unchanged.
Coupling-configuration symmetric case -We focus on the most symmetric case with φ h = φ c = φ in the r-symmetric configuration, where we find the simple relations that W 2 = r 2 W 1 , W 2 = r 2 W 1 , Φ = rφW 1 , and Φ = rφW 1 , yielding W d /W d = Φ/Φ. Then, the QME in Eq. (1) can be reduced to the QME with a single effective bath, defined by a single coherence parameter φ and a single rate W 1 . As is well known for the QME with a single bath, the system should reach a classical equilibrium state in the long-time limit. However, with degenerate energy levels, the off-diagonal (coherent) terms in the dissipation matrix Γ in Eq. (2) cannot be ignored even under the rotational wave approximation. It turns out that these coherent terms slow down the quantum dynamics significantly in general (|φ| < 1), approaching the classical steady state via a long lived quasi-stationary state with nonzero coherence.
First, we calculate the eigenvectors v i and the corresponding eigenvalues λ j of the Liouville matrix L. Details are given in Sec. ?? of SM [21]. We find the steady-state eigenvector , which corresponds to the classical fixed point. Other eigenvalues are negative except for |φ| = 1, thus the classical fixed point represents the unique steady state. At |φ| = 1, however, another eigenvector v 4 also has the zero eigenvalue, allowing multiple fixed points spanned by v 1 and v 4 . Note that |L ss | = r 2 (1 − φ 2 )(2W 1 + W 1 )W 1 from Eq. (9), which vanishes at these singular points of |φ| = 1.
with the probability conservation, the formal solution for P(t) reads where χ i depends on the initial condition P(0) in general. At |φ| = 1, λ 1 = λ 4 = 0, so the steady state P(∞) still depends on P(0). In Fig. 3, we display typical dynamic trajectories in the (ρ 12 , ρ 11 ) space with r = 1, starting from the empty initial condition of ρ i j (0) = 0 except for ρ 00 (0) = 1. As expected, all trajectories end up in the single (classical) fixed point in the long-time limit except for |φ| = 1, where the new coherent fixed point appears with ρ 12 (∞) 0. Note that the dynamics for φ 1 detours around the coherent fixed point for a significantly long time (known as a quasi-stationary state before), approaching the classical fixed point, thus it may be observed experimentally as a quasi-stationary state even in the presence of small decoherence. The additional zero eigenvalue (λ 4 = 0) at the singular points (|φ| = 1) implies another conservation law beside the probability conservation. Specifically, we find r 2ρ 11 +ρ 22 − rρ 12 − rρ 21 = 0 for φ = 1 from Eq. (4), or r 2 ρ 11 (t) + ρ 22 (t) − rρ 12 (t) − rρ 21 (t) = I 0 for all time t, where I 0 is a constant determined by the initial condition. We obtain the steady state solutions using Eq. (5) and the conservation law, written as ρ 11 (∞) = α − [rᾱ − 1−r 2 r α] ρ 12 (∞) and ρ 22 (∞) = α − [ᾱ r + 1−r 2 r α] ρ 12 (∞) with which depends on the initial state. In Fig. 3, we set r = 1 and I 0 = 0, so the coherent fixed point is determined by the intersection of two lines, ρ 11 = ρ 12 and ρ 11 = α −ᾱρ 12 . For I 0 0, the coherent fixed point is shifted along the curve of ρ 11 = α −ᾱρ 12 . In case of φ = −1, we get the same results except for changing the signs of ρ 12 and ρ 21 (see Eq. (??) of SM [21]). Note that the coherence can be finite and initialstate dependent even for ∆N = 0 (equilibrium). This may raise a doubt that the quantum current J The phenomena of multiple fixed points are observed not only in fermionic systems [17], but also in bosonic systems [19,24] which result from the existence of a dark state. In our case, the system state can be rewritten in a rotated orthonormal bases as |0 , |+ = (|1 + r|2 )/N r , and |− = (r|1 − |2 )/N r with N r = √ 1 + r 2 . Then, the system Hamiltonian is given asĤ S = E(|+ +| + |− −|) and the interaction Hamiltonian becomesĤ a SB = N r k g a 1b a † k |0 +| + h.c. at the singular points. Note that the state |− remains unchanged under the evolution operator, which corresponds to the dark state at φ = 1, i.e. any initial population in the dark state remains intact or −|ρ S |− = (r 2 ρ 11 + ρ 22 − rρ 12 − rρ 21 )/N 2 r should be conserved. We can easily extend our result to the degenerate multiple dots with multiple occupancy allowed. As the dark state decouples with baths, it may be useful to protect quantum information from decoherence [25].
We remark that the Lindblad description of degenerate quantum dots coupled to a single bath also yields multiple steady states with coherence at the maximum interference. This might be against the conventional wisdom that a system coupled to a single bath should reach the incoherent thermal equilibrium, regardless of its initial state. In this sense, the phenomenological introduction of φ is natural to guarantee the thermal steady state for |φ| < 1. Near the singular points, one may observe a long-living quasi-stationary state with the information of initial-state dependent coherent solutions.
Conclusion -We have investigated all possible steady-state solutions in the continuous quantum-dot heat engine coupled to terminals in parallel for various tunneling coefficients and interference strengths. The interference strength used in this work plays a similar role of the alignment of dipoles [26] in the bosonic system and acts as a source for decoherence from environments. We found that, unless the interference is completely negated, the steady states possess the coherence in general, which generates an extra quantum current, thus the engine performance can be enhanced in a specific region of the parameter space. Recently, the single quantum-dot (fermion) heat engine was realized experimentally [27]. The parallel-double-dot engine is also expected to be synthesized to confirm the enhancement of the QHE performance by thermal noises.

A. Model
We start with the Hamiltonian for the system (quantum dots) interacting with heat bathes, which are given byĤ whereĤ S ,Ĥ B , andĤ SB denote the Hamiltonians for the quantum-dot system, heat baths and interactions between the system and baths, respectively. The Hamiltonian of double quantum dots is given byĤ where E 1 and E 2 denote energies for dot 1 and dot 2, respectively, and E 12 is the Coulomb repulsion between electrons at dots. In the case of the degenerated dots, is the fermionic operator annihilating (creating) a single particle at dot d. We assume that only a single spinless fermion is allowed for each dot. Note that coherent hoppings between dots are not allowed.
The bath Hamiltonian is the sum of each bath HamiltonianĤ a B of bath a , which can be also written in terms of fermionic operators aŝ whereb a k (b a † k ) denotes the operator annihilating (creating) a particle with momentum k and energy ω a k in bath a (for simplicity, we assume here that the momentum is a scalar variable). The interaction Hamiltonian H SB is also the simple sum of the interaction HamiltonianĤ a SB for each bath a, which is also expressed with the fermionic operators aŝ which describes an electron hopping between dot d and bath a with a coupling strength g a d , usually depending on the momentum or energy.
In the limit of E 12 → ∞ (infinite repulsion), the simultaneous occupation at both dots is prohibited, thus the system state can be described with the three orthonormal bases of |0 (empty), |1 (single occupation in dot 1), and |2 (single occupation in dot 2). Then, the operatord d at dot d can be replaced by a jump operator |0 d|. Using these bases, we rewritê where E 0 means the empty-state energy (here, we set E 0 = 0).

B. Derivation of the QME
We derive the QME with an assumption that dots and baths are weakly coupled [1]. Instead of exploiting fermionic operators of dots, used in previous works [2][3][4], we use the jump operators. Starting from the von Neumann equation of the total system, ∂ tρ = −i Ĥ ,ρ , whereρ(t) is the density operator in the Schrödinger picture, the system dynamics expressed by the reduced density operator,ρ S = tr Bρ , is obtained by tracing out the bath degrees of freedom in the total system dynamic equation. In the weak coupling limit where the interaction Hamiltonian is small in comparison to the system and bath Hamiltonian, it is convenient to take the interaction picture, where the initial condition satisfies tr B [Ĥ SB ,ρ (0)] = 0. Substituting τ = t − s, we obtain Now we take the so-called Born-Markov approximation, where it is assumed thatρ (t) ≈ρ S (t)⊗ ρ B with the canonical heat bath density operator,

S3
with the temperature T a , the chemical potential µ a , and the number operatorn a = kb a † kb a k for each bath a, and the partition function Z = tr B e − a(Ĥ a B −µ ana) /T a (the Boltzmann constant is set as k B = 1). As the total density operator is given in the product form, this assumption implies that each bath is always in its equilibrium, regardless of the system evolution. This happens when the bath time scale τ a B is much smaller than the system time scale, thus the time scale separation between the system and baths is taken for granted, leading to the approximate replacement ofρ S (t − s) →ρ S (t).
Since the correlation tr B Ĥ SB (t), Ĥ SB (t − s),ρ (t) in Eq. (S9) may vanish for s τ a B , the integral upper bound can be extended to ∞, yielding a simpler approximate dynamic equation as (S11) Inserting Eq. (S7) into Eq.(S11), one can write each term in the commutation in Eq. (S11) as and the remainders are the Hermitian conjugates of Eqs. (S12) and (S13). Note that each bath contributes additively to Eq. (S11). and the correlators for bath a are defined by With the Fock-state description of bath particles in Eq. (S10), we find where N a is the Fermi-Dirac distribution in bath a, given as . (S16) Since the integral over time s in Eq. (S11) yields the delta function, i.e., a single mode for each bath satisfying ω a k = E survives in Eq. (S11) for the degenerate case with E 1 = E 2 = E. Note that we have omitted the Lamb shift correction, which is the order of E −1 , negligible in the high energy limit.
Changing k → N dk with a proper normalization N and integrating over k, we calculate the transition rates. First, consider the incoherent terms such as |d 0|ρ S |0 d| and |0 d|ρ S |d 0|.
For transitions between |0 and |d due to bath a, the transitions rates are obtained as where N a = 1 − N a and g a d = g a d (E). The + sign in Eq. (S18) stands for the transition from |0 to |d and the − sign stands for the opposite direction. Now we consider interference terms such as |0 2|ρ S |1 0| or |1 0|ρ S |0 2|. Due to the phase factor exp[±i(E 1 − E 2 )t], the interference terms vanish in long-time limit unless E 1 = E 2 (rotational wave approximation). In this work with E 1 = E 2 = E, we find the nonvanishing interference terms as where θ a is the difference of phase angles between g a 1 and g a 2 , defined as g a * 1 g a 2 = |g a 1 ||g a 2 |e iθ a . Defining the Lindblad operators aŝ the dynamic equation for the density operator in Eq. (S11) is rewritten as where { , } denotes the anticommutator. We introduce a phenomenological prefactor φ a for the interference terms in Eq. (S21) by replacing e iθ a → φ a with |φ a | ≤ 1 to take into account decoherence effects by other unknown environmental noises [5]. Note that the phase difference θ a is absorbed into φ a . In the Schrödinger picture withρ S (t) = e −iĤ S tρ S (t)e iĤ S t , Eq. (S21) is rewritten as a matrix form in Eqs. (??) and (??)of the main text.

S2. EIGENVECTORS AND EIGENVALUES OF THE LIOUVILLE OPERATOR
The density operator can be written in a form of vector: P = (ρ 00 , ρ 11 , ρ 22 , ρ 12 , ρ 21 , ρ 01 , ρ 02 , ρ 10 , ρ 20 ) T , with ρ i j = i|ρ S | j . Then, the equation of motion is given by ∂ t P = L tot P, where the Liouville operator L tot has a form of where the upper 5 × 5 block is given by L = a L a and the first term of the lower 4 × 4 block L irr = a L a irr . Each term in the summation is given as and the second term of the lower block E is It is easy to see that each 2 × 2 subblock of L a irr has negative eigenvalues only for |φ a | ≤ 1, thus ρ 01 , ρ 02 , ρ 10 , and ρ 20 , associated with L irr will vanish in long-time limit as the pure imaginary E contributes to a modulation only.
From now on, we focus on the 5×5 matrix L with the reduced vector P = (ρ 00 , ρ 11 , ρ 22 , ρ 12 , ρ 21 ) T , satisfying the dynamic equation ∂ t P = L P. For convenience, we take φ a as a real number. We introduce collective parameters for the sake of brevity as and then Eq. (??) of the main text is obtained.
From Eq. (??), we find the steady-state solution by inverting the 2 × 2 matrix L ss when its determinant |L ss | 0 as with The eigenvalues of the remaining two eigenvectors are the two roots of the characteristic equation of λ 2 + λ W 1 + W 1 + W 2 + W 2 + |L ss | = 0. Thus, we find the eigenvalues λ 4 and λ 5 as Note that |L ss | = 0 yields the additional zero eigenvalue, λ 4 = 0 and we expect multiple steady-state solutions given by a linear combination of v 1 and v 4 .

S3. STEADY-STATE CURRENTS
In this section, we calculate steady-state currents explicitly. The net particle currents J a d from bath a to dot d can be calculated from the Liouville equation in Eqs. (??) and (??) of the main text by sorting out the contributions to the time increment of the particle density of each dot d (ρ dd ) from each bath a. Then, we can easily identify where the first term in the right-hand-side represents particle transfer from bath a to the empty dot d, the second term represents particle transfer form the occupied dot d to bath a, and finally the third term represents the interference between relaxation channels to both dots. In the steady state, the particle density at dots is stationary, so the currents from both baths should be balanced in such . We can divide the particle current into the classical and quantum part as where the quantum current J q d is given as a product of the quantum speed Ψ d and the coherence ρ 12 (∞). The classical current is easily obtained by simply setting φ a = 0 in Eq. (S33) and Eq. (??) of the main text, and using the rates w a d+ = 2π|g a d | 2 N a and w a d− = 2π|g a d | 2 N a , yielding where |L 0 | = |L ss | φ a =0 = W 1 W 2 + W 1 W 2 + W 1 W 2 and the external bias ∆N = N h − N c > 0. Note that J cl d is always positive. The quantum part, J q d = Ψ d ρ 12 (∞), is also obtained from Eq. (S33) and Eq. (??) of the main text, yielding and the coherence term ρ 12 (∞) = ρ 21 (∞) is obtained from Eq. (S27), after some algebra, as which is valid except for the singular case of |L ss | = 0 (see Sec. S4 for the singular case). In contrast to the classical current, the quantum current J q d can be both positive and negative, which can vanish either by the zero quantum speed (Ψ d = 0) or by the zero coherence (ρ 12 = 0). In Note that the quantum conductance in Eq. (S40) cannot be positive from , implying the current enhancement is not possible in the linear response regime.
Next, we will investigate the nonlinear regime in the r-symmetric configuration (g a 2 = rg a 1 ), where the algebra becomes simplified. From Eqs. (S35) and (S36), we can easily see that the classical currents for two dots are simply related as J cl 2 = r 2 J cl 1 . For convenience, we set the tunnelling coefficients as |g a 1 | 2 = k a 1 + r 2 , |g a 2 | 2 = r 2 k a 1 + r 2 , which satisfies the r-symmetric condition. Then, we find where N a = 1 + N a . After some algebra, we also find the relative quantum conductance σ q d /σ cl d as where L reads In the expansion of σ q d /σ cl d = S 0 d + S 1 d ∆N + O(∆N 2 ), the leading terms are given by Setting φ c = φ − ∆φ and φ h = φ with a finite φ, we obtain and where O and O are higher order terms and the expansion of L eq yields L eq ≈ 1 − N 2 k h + k c 2 1 − φ 2 + 2k c k h + k c φ∆φ .
We find that both linear coefficients are negative and S 0 1 = r 2 S 0 2 . For r > 1, the negative contribution from dot 1 (weaker coupling) is stronger. As the second-order coefficients are positive for φ∆φ > 0 and stronger than the linear coefficients for very small ∆φ, the nonlinear contribution may overcome the linear response to make the quantum conductance positive.
Approaching to φ h = φ c = ±1, the leading order of L eq becomes linear in ∆φ, thus S 1 d remains finite, while S 0 d goes to zero. Therefore, a strong enhancement of the current is expected.

S4. FULLY SYMMETRIC CASE
We consider the most symmetric case with φ h = φ c = φ in the r-symmetric configuration, where we find simple relations as W 2 = r 2 W 1 , W 2 = r 2 W 1 , Φ = rφW 1 , and Φ = rφW 1 . Then, the