**2. Theoretical Model**

### *2.1. Description of the System and Underlying Framework*

We consider the case of two conductors that are mutually connected via the Coulomb interaction. Each conductor consists of a quantum dot with a single level active for transport. We omit spin indices due to spin degeneracy. Besides, we consider a large on-site Coulomb interaction that prevents the double occupancy in each dot. Each quantum dot is tunnel-coupled to two electronic reservoirs

that can be biased with electrostatic and thermal gradients. Each tunneling barrier is modeled by capacitors denoted by *Ci* with *i* = 1, 2, 3, 4. As mentioned above, the two quantum dots interact electrostatically through a capacitor *C*. A sketch for this system is depicted in Figure 1b. Under these circumstances, we describe the system using four possible charge states |0 = |<sup>0</sup>*u*0*d*, |*u* = |<sup>1</sup>*u*0*d*, |*d* = |<sup>0</sup>*u*1*d*, and |2 = |<sup>1</sup>*u*1*d*, where *nund* denotes the charge state with *nu* electrons in the upper dot and *nd* electrons in the lower dot. For simplicity, we consider an isothermal configuration in which all reservoirs are held at a common temperature *T*. We also keep different bias voltages *Vi* applied to the four terminals.

We are interested in the charge and heat transport in the sequential tunneling regime, in which the tunneling rate (denoted by Γ) satisfies *h*¯ Γ *kBT*. In this regime, transport of electrons along each quantum dot occurs in a sequence of one electron transfer event at each time. Electrons can hop into a quantum dot, and then relax before they jump again. This restriction eliminates the transitions |0(2)→|2(0) and |*u*(*d*)→|*d*(*u*). Additionally, we consider that there is no particle transfer from one dot to another by tunneling. The only interaction between the dots is then due to their mutual influence caused by the electrostatic interactions.

**Figure 1.** (**a**) Double quantum dot capacitively coupled to four terminals held at potentials *Vi* and temperatures *Ti*, for *i* = 1, 2, 3, 4. The transition rates Γ±*i* and *γ*±*i* for each barrier are described in the main text. (**b**) Electrostatic sketch showing the capacitors and voltages involved in the description of the energy levels of the quantum dots.

The theoretical framework employed to describe the quantum transport in our system is called *stochastic thermodynamics* [14,20–22]. Quite generally, we can consider a setup with an arbitrary number of states *n* ∈ {1, 2, ..., *N*} and picture each state as a node in a connected network. We draw edges **e** connecting states between which a transition may occur, and require these to be possible in both directions. However, transitions along ±**e** are not required to happen at the same rate or with the same probability. Note that two nodes may be connected with several edges if there are various physical mechanisms through which the system can transition between the associated states. The evolution of the system is modeled as a Markov jump process, i.e., the probability that the system jumps from one state to another is independent of its previous history. This evolution can also be visualized as a random walk on the network. A physical model is defined by prescribing the forward and backward transition rates *w*±**e**, which evidently may be functions of the physical parameters involved. The fluctuating current along an edge **e**, *j***e**(*t*) = ∑*k δ*(*t* − *tk*)(*<sup>δ</sup>*+**e**,**e***<sup>k</sup>* − *<sup>δ</sup>*−**e**,**e***<sup>k</sup>* ), is a stochastic variable that peaks if the system transitions along the directed edge **e***k* at time *tk*. Physical currents, i.e., currents associated to the transport of physical quantities such as charge or heat, are weighted currents *Jα* = ∑**e** *dα***e** *j***e**, where *<sup>d</sup><sup>α</sup>*+**e** = −*d<sup>α</sup>*−**e** specifies the amount of a physical variable *α* exchanged with an external reservoir along a transition edge **e**.

When applying the previous theoretical treatment to our particular system, we consider that the tunneling rates depend on the energy of the system. Specifically, we consider the value of Γ*i* for the tunneling of electrons between a reservoir *i* and a quantum dot whenever the other dot is empty, and *γi* when the other dot is occupied. Then, the transition rates (previously called *<sup>w</sup>*±**e**) are thus dependent on the dot charge states. The transition rates are defined according to Fermi's golden rule as

$$
\Gamma\_i^- = \Gamma\_i f(\mu\_{\ell,0} - qV\_i) \tag{1}
$$

$$
\Gamma\_i^+ = \Gamma\_i \left( 1 - f(\mu\_{\ell,0} - qV\_i) \right) \tag{2}
$$

$$
\gamma\_i^- = \gamma\_i f(\mu\_{\ell,1} - qV\_i) \tag{3}
$$

$$
\gamma\_i^+ = \gamma\_i \left( 1 - f(\mu\_{\ell,1} - qV\_i) \right) \tag{4}
$$

where *f*(*x*) = 1 + *e<sup>x</sup>*/*kBT*−<sup>1</sup> is the Fermi–Dirac distribution function, the − superscript stands for the tunneling from the lead to the dot, and + for the reverse process. The transition rates are schematically represented in a network diagram in Figure 2. In our arrangement, we take the *up* dot (- = *u*) connected to left and right reservoirs with *i* = {1, <sup>2</sup>}, and the *down* dot (- = *d*) to reservoirs with *i* = {3, <sup>4</sup>}. Note that the numerical subindex in the previous transition rates thus indicates the reservoir involved in the transition, as shown in Figure 1a. The chemical potential for the dot -, i.e., *μ*-,0 (*μ*-,1), corresponds to the situation in which the *other* dot is empty (occupied).

**Figure 2.** Scheme for the transition rates between the two dot states.

To determine the effective chemical potentials of the dots, we must develop a model that takes into account how their energy levels are influenced by electrostatic interactions. When interactions are properly included as in our description, all currents are gauge invariant, as they depend only on voltage differences. Hereafter, we shorten the notation and define *Vij* ≡ *Vi* − *Vj*. Under these considerations, the dot levels become

$$
\varepsilon\_{\rm u,n} \rightarrow \mu\_{\rm u,n} = \varepsilon\_{\rm u} + \mathcal{U}\left(1,0\right) - \mathcal{U}\left(0,0\right) + E\_{\mathbb{C}}\delta\_{1n} \tag{5}
$$

$$
\varepsilon\_{d,n} \rightarrow \mu\_{d,n} = \varepsilon\_d + \mathcal{U}\left(0,1\right) - \mathcal{U}\left(0,0\right) + E\_{\mathbb{C}}\delta\_{1n} \tag{6}
$$

where *εu* and *εd* are the bare energy levels, and *n* = 0 (1) corresponds to the case where *other* dot is empty (occupied). Here, *EC* = <sup>2</sup>*q*<sup>2</sup>*C*/ *<sup>C</sup>*Σ*uC*Σ*<sup>d</sup>* − *C*<sup>2</sup> is the charging energy with *<sup>C</sup>*Σ*d*(*u*) = *<sup>C</sup>*1(3) + *<sup>C</sup>*2(4) + *C*. The chemical potential *μu*(*d*),*<sup>n</sup>* is defined as the change in the electrostatic energy when the charge number *Nu*(*d*) changes by one when the dot *d*(*u*) is either empty (*n* = 0) or occupied by one electron (*n* = 1). The electrostatic energy is computed from *<sup>U</sup>*(*Nu*, *Nd*) = ∑*i qNi* 0 *dQi φi*(*Qi*) where *φi* is the internal potential in each quantum dot obtained by means of elementary electrostatic relations. Then, the arguments of the Fermi functions appearing in the tunneling rates read [16]:

$$
\mu\_{\boldsymbol{\mu},\boldsymbol{n}} - qV\_1 = \varepsilon\_{\boldsymbol{n}} + \frac{1}{\mathbb{C}\_{\Sigma\boldsymbol{n}}\mathbb{C}\_{\Sigma\boldsymbol{d}} - \mathbb{C}^2} \left[ \frac{q^2}{2} \mathbb{C}\_{\Sigma\boldsymbol{d}} + q \left( \mathbb{C}\_{\Sigma\boldsymbol{d}} \mathbb{C}\_2 V\_{21} + \mathbb{C} \mathbb{C}\_3 V\_{31} + \mathbb{C} \mathbb{C}\_4 V\_{41} \right) \right] + E\_{\mathbb{C}} \delta\_{1\boldsymbol{n}} \tag{7}
$$

$$
\mu\_{\rm u,n} - qV\_2 = \varepsilon\_{\rm u} + \frac{1}{\mathbb{C}\_{\Sigma \rm d} \mathbb{C}\_{\Sigma d} - \mathbb{C}^2} \left[ \frac{q^2}{2} \mathbb{C}\_{\Sigma d} + q \left( \mathbb{C}\_{\Sigma d} \mathbb{C}\_1 V\_{12} + \mathbb{C} \mathbb{C}\_3 V\_{32} + \mathbb{C} \mathbb{C}\_4 V\_{42} \right) \right] + E\_{\mathbb{C}} \delta\_{1n} \tag{8}
$$

$$p\mu\_{d,\mathfrak{n}} - qV\_3 = \varepsilon\_d + \frac{1}{\mathbb{C}\_{\Sigma\mathfrak{u}}\mathbb{C}\_{\Sigma d} - \mathbb{C}^2} \left[ \frac{q^2}{2} \mathbb{C}\_{\Sigma\mathfrak{u}} + q \left( \mathbb{C}\_{\Sigma\mathfrak{u}}\mathbb{C}\_4 V\_{43} + \mathbb{C}\mathbb{C}\_1 V\_{13} + \mathbb{C}\mathbb{C}\_2 V\_{23} \right) \right] + E\_\mathbb{C} \delta\_{1\mathfrak{n}} \tag{9}$$

$$\mu\_{d,n} - qV\_4 = \varepsilon\_d + \frac{1}{\mathbb{C}\_{\Sigma u}\mathbb{C}\_{\Sigma d} - \mathbb{C}^2} \left[ \frac{q^2}{2}\mathbb{C}\_{\Sigma u} + q \left( \mathbb{C}\_{\Sigma u}\mathbb{C}\_3 V\_{34} + \mathbb{C}\mathbb{C}\_1 V\_{14} + \mathbb{C}\mathbb{C}\_2 V\_{24} \right) \right] + E\_\mathbb{C} \delta\_{1n} \tag{10}$$

that now depend only on voltage differences. The four *μ<sup>u</sup>*,*<sup>n</sup>* − *<sup>V</sup>*1(2), and *μd*,*<sup>n</sup>* − *<sup>V</sup>*3(4) are the electrochemical potentials. We take *V*12, *V*13 and *V*34 as the only independent biases, since the rest of voltage differences can be expressed as linear combinations of these values.

As discussed above, we apply the Markov approximation in order to determine the dynamics of the probabilities of finding the system in one of the four states. Specifically, we employ the master equation formalism, where the time evolution of the system is governed by a master equation that gives the probability distribution of the considered stochastic variables in terms of the transition rates between the different states. Defining Γ<sup>±</sup>*u*(*d*) = <sup>Γ</sup>±1(3) + <sup>Γ</sup>±2(4), the following relations are found:

$$
\begin{pmatrix}
\dot{p}\_0 \\ \dot{p}\_u \\ \dot{p}\_d \\ \dot{p}\_2
\end{pmatrix} = \begin{pmatrix}
\Gamma^-\_u & -\Gamma^+\_u - \gamma\_d^- & 0 & \gamma\_d^+ \\
\Gamma^-\_d & 0 & -\gamma\_u^- - \Gamma^+\_d & \gamma\_u^+ \\
0 & \gamma\_d^- & \gamma\_u^- & -\gamma\_u^+ - \gamma\_d^+
\end{pmatrix} \begin{pmatrix} p\_0 \\ p\_u \\ p\_d \\ p\_2 \end{pmatrix} \tag{11}
$$

As we are interested in the steady state, we set all *p*˙*i* = 0. Considering the normalization condition ∑*i pi* = 1, we obtain:

$$\mathcal{p}\_0 = \frac{1}{a} \left[ \Gamma\_d^+ \gamma\_\mu^+ \left( \Gamma\_\mu^+ + \gamma\_d^- \right) + \Gamma\_\mu^+ \gamma\_d^+ \left( \Gamma\_d^+ + \gamma\_\mu^- \right) \right] \tag{12}$$

$$p\_{\mu} = \frac{1}{a} \left[ \Gamma\_{\mu}^{-} \Gamma\_{d}^{+} \left( \gamma\_{u}^{+} + \gamma\_{d}^{+} \right) + \gamma\_{u}^{-} \gamma\_{d}^{+} \left( \Gamma\_{u}^{-} + \Gamma\_{d}^{-} \right) \right] \tag{13}$$

$$\mathcal{V}\_d = \frac{1}{a\_\cdot} \left[ \Gamma\_u^+ \Gamma\_d^- \left( \gamma\_u^+ + \gamma\_d^+ \right) + \gamma\_u^+ \gamma\_d^- \left( \Gamma\_u^- + \Gamma\_d^- \right) \right] \tag{14}$$

$$p\_2 = \frac{1}{a} \left[ \gamma\_u^- \gamma\_d^- \left( \Gamma\_u^- + \Gamma\_d^- \right) + \Gamma\_u^- \Gamma\_d^+ \gamma\_d^- + \Gamma\_u^+ \Gamma\_d^- \gamma\_u^- \right] \tag{15}$$

with

$$\begin{split} \boldsymbol{u} &= \Gamma\_{\boldsymbol{u}}^{-} \left[ \Gamma\_{\boldsymbol{d}}^{+} \left( \gamma\_{\boldsymbol{u}}^{+} + \gamma\_{\boldsymbol{d}}^{+} \right) + \gamma\_{\boldsymbol{d}}^{-} \left( \gamma\_{\boldsymbol{u}} + \Gamma\_{\boldsymbol{d}}^{+} \right) \right] + \Gamma\_{\boldsymbol{u}}^{+} \Gamma\_{\boldsymbol{d}}^{+} \left( \gamma\_{\boldsymbol{u}}^{+} + \gamma\_{\boldsymbol{d}}^{+} \right) + \Gamma\_{\boldsymbol{d}} \gamma\_{\boldsymbol{d}}^{-} \gamma\_{\boldsymbol{u}}^{+} + \Gamma\_{\boldsymbol{u}} \gamma\_{\boldsymbol{d}}^{+} \gamma\_{\boldsymbol{u}}^{-} + \\ &+ \gamma\_{\boldsymbol{d}}^{-} \left[ \Gamma\_{\boldsymbol{u}}^{+} \left( \gamma\_{\boldsymbol{u}} + \gamma\_{\boldsymbol{d}}^{+} \right) + \gamma\_{\boldsymbol{u}}^{-} \gamma\_{\boldsymbol{d}} \right] \end{split} \tag{16}$$

and <sup>Γ</sup>*u*(*d*) = Γ<sup>+</sup>*u*(*d*) + Γ−*u*(*d*) (similar for *<sup>γ</sup>u*(*d*)).

We now compute the electrical current *I*1 that flows between the first lead and the upper dot, which, we from now on call *drag current* for historical reasons (note that since generally *V*12 = 0 it is not a current arising solely from the drag effect). This current is obtained by weighting the transition probabilities with the electron charge *q*. The result is

$$I\_1 = q \left( \Gamma\_1^- p\_0 - \Gamma\_1^+ p\_u + \gamma\_1^- p\_d - \gamma\_1^+ p\_2 \right) \tag{17}$$

Because of electric charge conservation, we immediately know *I*2 = −*I*<sup>1</sup> = *I* for the current between the second terminal and the *up* dot (we assign a + sign whenever the current flows from a lead into a dot, and a − sign otherwise). We can also compute the heat current by weighting the transitions with the amount of transferred effective energy (the electrochemical potential),

$$J\_1 = \tilde{\mu}\_{\rm u,0} \left( \Gamma\_1^- p\_0 - \Gamma\_1^+ p\_u \right) + \tilde{\mu}\_{\rm u,1} \left( \gamma\_1^- p\_d - \gamma\_1^+ p\_2 \right) \tag{18}$$

where *μ*˜*u*,*<sup>n</sup>* = *μ<sup>u</sup>*,*<sup>n</sup>* − *qV*1. Similar expressions are obtained for the rest of the *Ji*. Energy conservation leads to *J*1 + *J*2 + *J*3 + *J*4 = *I*1*V*21 + *I*3*V*43. These currents were investigated by Sánchez et al. [16] when the *up* dot is at equilibrium with *V*1 = *V*2; a nonzero drag current *I*1 then appears when Γ1*γ*2 = *γ*1Γ2. This means that the current in the lower terminals (drive system) drives the upper dot towards a

non-equilibrium situation by the appearance of a drag current. The drag phenomenon can be clearly understood from the following Joule relation found for this setup

$$\begin{array}{ccccc}J\_1 + J\_2 + J\_c & = & -I\_1V\_{12} \\ J\_3 + J\_4 - J\_c & = & -I\_3V\_{34} \end{array} \tag{19}$$

where

$$J\_{\mathbb{C}} = \frac{E\_{\mathbb{C}}}{\mathfrak{a}} (\gamma\_{\mu}^{+} \gamma\_{d}^{-} \Gamma\_{\mu}^{-} \Gamma\_{d}^{+} - \gamma\_{\mu}^{-} \gamma\_{d}^{+} \Gamma\_{\mu}^{+} \Gamma\_{d}^{-}) \tag{20}$$

is the heat flow between the drag and the drive system. This expression generalizes the relation found in Reference [23] for a three-terminal double quantum dot. In such a system, the drag conductor is connected to two reservoirs, whereas the drive dot is coupled to a single contact. Therefore, the drive subsystem does not support any charge current. Under these considerations, the drive dot carries a heat flow *J* (with a similar form to Equation (20), which is proportional to the drag charge current when Γ2 = *γ*1 = 0, i.e. in the so-called strong coupling regime. Here, Equation (19) demonstrates the existence of a heat flow *Jc* between the drag and drive subsystems. This energy flow appears in addition to the heat flows *J*1, *J*2 through the drag conductor and the heat currents *J*3, *J*4 in the drive subsystem, even when they are held at common temperature and no particle transfer exists between them.

Returning to our purpose, which is to find a route to an effective equilibrium state, we address the issue of whether the opposite phenomenon to the drag is possible, i.e., if we can achieve a non-equilibrium configuration with *V*1 = *V*2 for which the drag effect causes the stalling of the upper currents. Under this novel situation, we check whether our system reaches an "effective linear response regime" by testing the microreversibility property through the Onsager relations and the fulfilment of the FDT. To this end, we focus on the *up* dot and consider three stalling configurations: (i) when both the electrical and the heat flow vanish, i.e., *I*1 = 0 and *J*1 = 0, which we call the globally stalled scenario; (ii) when the charge current is nullified, *I*1 = 0, but there is a finite heat flow *J*1 = 0; and (iii) when there is a finite electrical current *I*1 = 0 but no heat flow, *J*1 = 0. These situations correspond to the locally charge-stalled and heat-stalled cases, respectively. The simplest manner to achieve the globally stalled case is tuning the system to the strong coupling configuration by setting *γ*1 = Γ2 = 0. Under this situation, electrons can only tunnel in and out of the top-left reservoir if the lower dot is empty, and of the top-right reservoir if the lower dot is occupied.

### *2.2. Detailed Balance and Behavior at Equilibrium*

Before presenting our results, we carefully revise the behavior of systems near thermodynamic equilibrium. In this situation, all existing currents in a system tend to zero on average. This behavior is called *global detailed balance*. According to statistical mechanics, systems subject to these conditions exhibit the property that the correlations of the spontaneous fluctuations and the dissipative response to an external perturbation obey the same rules, which is primarily known as Onsager's regression hypothesis [14]. This important statement is the heart of the fluctuation–dissipation theorem (FDT). If we consider an arbitrary physical current *Jα* (such as a heat or charge current) and its affinity or conjugate force *hα* (which in these cases would correspond to gradients in temperature or electrical potential, respectively), the theorem can be expressed as

$$D\_{h\_{\mathfrak{a}}}I\_{\mathfrak{a}}\left(\mathfrak{x}^{\mathfrak{e}\mathfrak{q}}\right) = D\_{\mathfrak{a},\mathfrak{a}}\left(\mathfrak{x}^{\mathfrak{e}\mathfrak{q}}\right) \tag{21}$$

where *D<sup>α</sup>*,*<sup>α</sup>* is a generalized diffusion constant proportional to *Jα Jα*. The vector *x* contains all the parameters the current may depend on, and satisfies *Jα* (*xeq*) = 0 for all currents in the system; their conjugate forces are evidently also required to vanish. The previous equation can be generalized in such a way that it expresses the FDT for the combination of two currents and their conjugate

forces by changing one index *α* for a different one and symmetrizing both sides of the expression (see complementary material of Reference [14]).

Another major result in thermodynamics close to equilibrium is found in Onsager's reciprocal relations (RRs), which actually follow from the FDT if the system enjoys the property of being time-reversible [15]. In the following, we restrict ourselves to relations between heat and charge currents, following Onsager's original article [1]. For a system where transport of these quantities exists, the mechanisms are usually not independent, but interfere with each other leading to the well known thermoelectric effects. If we consider a system at equilibrium, small fluctuations or external perturbations may allow for the transport of small quantities of charge and heat while the system is returning to its original state. Onsager established that, in these situations, the responses of a current due to a variation of the other current's conjugate force are equal, i.e. the heat current responds in the same way to a variation of the electrical potential as the charge current to a temperature fluctuation. This result is best visualized by writing the currents in matrix form. For a simple system with a single heat and charge current, we have:

$$
\begin{pmatrix} J\_{\text{charge}} \\ J\_{\text{heat}} \end{pmatrix} = \begin{pmatrix} L\_{11} & L\_{12} \\ L\_{21} & L\_{22} \end{pmatrix} \begin{pmatrix} \delta \left( \Delta V \right) \\ \delta \left( \Delta T / T \right) \end{pmatrix} \tag{22}
$$

where *L*11 and *L*22 are the electrical and thermal conductances, and *L*12 = *∂*Δ*T*/*T Jcharge* and *L*21 = *∂*Δ*V Jheat* represent the electrothermal and thermoelectrical coefficients that arise from the interference of the two transport mechanisms. Onsager's statement is then equivalent to the requirement that the conductance matrix be symmetric, *L*12 = *L*12. In addition to these relations, the scattering theory formalism ensures that both the thermal and the electrical conductances are semipositive.

Despite these theorems being major cornerstones in our understanding of the behavior of systems obeying global detailed balance, most complex systems live out of equilibrium. Accordingly, similar relations have been sought for systems where detailed balance is explicitly broken, since their finding would allow us to characterize and study out-of-equilibrium systems in a similar manner as when detailed balance is satisfied.

### *2.3. Local Detailed Balance and Equilibrium-Like Relations*

A central assumption in stochastic thermodynamics far from equilibrium, when global detailed balance is not satisfied, is local detailed balance (LDB). It relates the forward and backward transition rates *w* into and out of a state *A* by means of a mechanism *ν* and reads [24]

$$\frac{w\_{A \to B}^{\nu}}{w\_{B \to A}^{\nu}} = \varepsilon^{-\beta\_{\nu}\Delta x} \tag{23}$$

where *βν* is the inverse temperature of the reservoir involved in the transition and Δ*ε* is the difference between the energies of states *A* and *B*. It can be easily checked that the rates in Equations (1)–(4) indeed satisfy the LDB condition.

In Reference [14], it was reported that, if LDB is satisfied in a system driven arbitrarily far from equilibrium, its response to a perturbation or a spontaneous fluctuation may obey a relation similar to the equilibrium FDT if certain additional conditions are fulfilled. More precisely, it has been established that a current *Jα* in such a system obeys Equation (21) with *xeq* replaced by *<sup>x</sup>st*, where *xst* corresponds to a configuration of the parameters of the current such that *Jα xst* = 0, i.e., the considered current stalls. This is valid if the force *hα* couples exclusively to those transitions that contribute to the conjugate current *J<sup>α</sup>*. It is important to notice the difference between this statement and the first FDT valid only near equilibrium, since we now only require a given current to stall internally. This may be a consequence of the appropriate tuning of the rest of the currents in the system, which are no longer required to vanish, and can in fact assume arbitrary magnitudes.

Similarly, Onsager's reciprocal relations have also been extended to non-equilibrium situations, under the condition of a marginal time-reversibility [25]. Again, it is required that the currents stall in order for the RRs to hold far from equilibrium.
