**Appendix A. Bisection Method**

To find the roots of the drag current at a certain point of parameter space, we make extensive use of the bisection algorithm. Firstly, we fix the value of the biases *V*12 and *V*13. We are now interested in finding *V*34 such that *I*1 (*<sup>V</sup>*12, *V*13, *<sup>V</sup>*34) = 0. Let *f*(*x*) = *I*1 (*<sup>V</sup>*12, *V*13, *<sup>x</sup>*). We iteratively seek two points *a*

and *b* such that *f*(*a*)*f*(*b*) < 0. Since *I*1 is continuous, which implies there is a root lying in the interval (*a*, *b*). To find this root with a tolerance , once we have localised such points, where  represents the largest value for the width of the interval centered at the point that we ultimately accept as a root, we proceed as follows:


Despite the fact that the width of the interval decreases only linearly with this method, thus making it rather slow, it is always ensured to converge if there is a root lying inside the original interval. Furthermore, for a highly nonlinear function, which may present a behavior that is difficult to picture, such as the drag current, this method is far superior to non-fixed interval algorithms such as the Newton method.

### **Appendix B. Full Counting Statistics and Computation of Cumulants**

Despite this method being extensively discussed in References [27,28], a good understanding of it has been crucial for this work. However, several different approaches are found in the literature, and therefore we consider it adequate to clarify the precise path we have taken. The derivation that follows closely follows the one found in the additional material of Reference [26].

The central quantity we are involved with is *P* ({*<sup>n</sup>*1, *n*2, *n*3, *<sup>n</sup>*4}; *t*) ≡ *P* ({*n*}; *t*), the probability that *nj* electrons have been transferred through the terminal *j*. The associated cumulant generating function (CGF) F ({*χ*}; *t*) follows from

$$\exp\left[\mathcal{F}(\{\chi\};t)\right] = \sum\_{\{n\}} P\left(\{n\};t\right) \exp\left(i\sum\_{j} \chi\_{j} n\_{j}\right) \tag{A1}$$

From the CGF, we can obtain the desired cumulants by taking partial derivatives with respect to the counting fields *χj* at *χj* = 0:

$$\mathcal{C}\_{pqrs} = \left. \partial\_{i\chi\_1}^p \partial\_{i\chi\_2}^q \partial\_{i\chi\_3}^r \partial\_{i\chi\_4}^s \mathcal{F}(\{\chi\};t) \right|\_{\chi=0} \tag{A2}$$

Then, the current cumulants in the long-time limit are simply given by

$$\left\langle I\_1^p I\_2^q I\_3^r I\_4^s \right\rangle = (-q)^{p+q+r+s} \frac{\text{dC}\_{pqrs}}{\text{d}t} \bigg|\_{t \to \infty} \tag{A3}$$

However, in general, the expression for the CGF is difficult to obtain and handle, and must be computed recursively. We first introduce the operator L(*χ*) as the Fourier transform of the entries of our matrix L(0) ≡ L satisfying *ρ*˙ = L*ρ*. When only sequential tunneling is considered, it is obtained by adding counting fields to the off-diagonal entries of L, with a plus (minus) sign when the transition corresponds to an electron entering (leaving) the corresponding lead (we skip the details). For our matrix (Equation (11)), it assumes the form

$$
\mathcal{L}(\chi) = \begin{pmatrix}
\Gamma^{-}\_{1}\epsilon^{-i\chi\_{1}} + \Gamma^{-}\_{2}\epsilon^{-i\chi\_{2}} & -\Gamma^{+}\_{u} - \gamma^{-}\_{d} & 0 & \gamma\_{3}^{+}\epsilon^{j\chi\_{3}} + \gamma\_{4}^{+}\epsilon^{j\chi\_{4}} \\
\Gamma^{-}\_{3}\epsilon^{-i\chi\_{3}} + \Gamma^{-}\_{4}\epsilon^{-i\chi\_{4}} & 0 & -\gamma\_{u}^{-} - \Gamma^{+}\_{d} & \gamma\_{1}^{+}\epsilon^{i\chi\_{1}} + \gamma\_{2}^{+}\epsilon^{i\chi\_{2}} \\
0 & \gamma\_{3}^{-}\epsilon^{-i\chi\_{3}} + \gamma\_{4}^{-}\epsilon^{-i\chi\_{4}} & \gamma\_{1}^{-}\epsilon^{-i\chi\_{1}} + \gamma\_{2}^{-}\epsilon^{-i\chi\_{2}} & -\gamma\_{u}^{+} - \gamma\_{d}^{+} \\
\end{pmatrix} \tag{A4}$$

In the long-time limit, the CGF can be written as

$$\mathcal{F}\left(\{\chi\};t\right) = \lambda\_0(\chi)t \tag{A5}$$

where *<sup>λ</sup>*0(*χ*) is the minimum eigenvalue of L(*χ*) [27]. If we are able to compute it, then the current correlations easily follow from Equation (A3) as

$$\left\langle I\_1^p I\_2^q I\_3^r I\_4^s \right\rangle = (-q)^{p+q+r+s} \left. \partial\_{i\chi\_1}^p \partial\_{i\chi\_2}^q \partial\_{i\chi\_3}^r \partial\_{i\chi\_4}^s \right\vert\_{\chi\_4} \lambda\_0(\chi) \bigg\vert\_{\chi=0} \tag{A6}$$

To calculate *<sup>λ</sup>*0(*χ*), we have employed Flidnt's method, which is discussed now. First, we write L(*χ*) as

$$
\mathcal{L}(\chi) = \mathcal{L} + \mathcal{L}(\chi) \tag{A7}
$$

where L = L(0) and <sup>L</sup>˜(*χ*) = L(*χ*) − L. Next, we define operators P = |00| and Q = 1 − P, where |0 = (*p*0, *pu*, *pd*, *<sup>p</sup>*2)*<sup>T</sup>* and 0| = (1, 1, 1, 1) are the left and right null eigenvectors of L. Clearly, PL = LP = 0 and QL = LQ = L. To determine the CGF from Equation (A5), we must solve

$$\mathcal{L}(\chi)|0(\chi)\rangle = \left[\mathcal{L} + \mathcal{L}(\chi)\right]|0(\chi)\rangle = \lambda\_0(\chi)|0(\chi)\rangle \tag{A8}$$

By choosing 0| <sup>0</sup>(*χ*) = 1, it follows that

$$<\langle 0|\lambda\_0(\chi) - \mathcal{L}|0(\chi)\rangle = \lambda\_0(\chi) = \langle 0|\mathcal{L}(\chi)|0(\chi)\rangle \tag{A9}$$

and using Q on |0(*χ*) we also find

$$|0(\chi)\rangle = |0\rangle + \mathcal{Q}|0(\chi)\rangle\tag{A10}$$

From Equation (A8), using that L and Q commute and Q<sup>2</sup> = Q, we obtain

$$\mathcal{Q}|0(\chi)\rangle = \mathcal{Q}\left[\lambda\_0(\chi) - \mathcal{L}\right]^{-1} \mathcal{Q}\mathcal{L}(\chi)|0(\chi)\rangle \tag{A11}$$

We now define

$$\mathcal{R}\left[\lambda\_0(\chi)\right] = \mathcal{Q}\left[\mathcal{L} - \lambda\_0(\chi)\right]^{-1}\mathcal{Q} \tag{A12}$$

and substitute Equation (A11) into Equation (A10) to find

$$|0(\chi)\rangle = |0\rangle - \mathcal{R}\left[\lambda\_0(\chi)\right]\mathcal{L}(\chi)|0(\chi)\rangle\tag{A13}$$

so that finally

$$|0(\chi)\rangle = \left\{1 + \mathcal{R}\left[\lambda\_0(\chi)\right] \mathcal{L}(\chi)\right\}^{-1} |0\rangle \tag{A14}$$

and therefore using Equation (A9) we arrive at

$$\lambda\_0(\chi) = \langle 0 \vert \mathcal{L}(\chi) \left\{ 1 + \mathcal{R} \left[ \lambda\_0(\chi) \right] \mathcal{L}(\chi) \right\}^{-1} \vert 0 \rangle \tag{A15}$$

We now Taylor expand the previous expression. In our case of four counting fields, L ˜(*χ*) is expanded as

$$\begin{split} \mathcal{L}(\boldsymbol{\chi}) &= \mathcal{L}^{(1,0,0,0)}(i\chi\_{1}) + \mathcal{L}^{(0,1,0,0)}(i\chi\_{2}) + \mathcal{L}^{(0,0,1,0)}(i\chi\_{3}) + \mathcal{L}^{(0,0,0,1)}(i\chi\_{4}) + \\ &+ \frac{1}{2!} \left[ \tilde{\mathcal{L}}^{(2,0,0,0)}(i\chi\_{1})^{2} + \tilde{\mathcal{L}}^{(0,2,0,0)}(i\chi\_{2})^{2} + \tilde{\mathcal{L}}^{(0,0,2,0)}(i\chi\_{3})^{2} + \tilde{\mathcal{L}}^{(0,0,0,2)}(i\chi\_{4})^{2} + \\ &+ 2\mathcal{L}^{(1,1,0,0)}(i\chi\_{1})(i\chi\_{2}) + 2\mathcal{L}^{(1,0,1,0)}(i\chi\_{1})(i\chi\_{3}) + 2\mathcal{L}^{(1,0,0,1)}(i\chi\_{1})(i\chi\_{4}) + \\ &+ 2\mathcal{L}^{(0,1,1,0)}(i\chi\_{2})(i\chi\_{3}) + 2\mathcal{L}^{(0,1,0,1)}(i\chi\_{2})(i\chi\_{4}) + 2\mathcal{L}^{(0,0,1,1)}(i\chi\_{3})(i\chi\_{4}) \right] + \mathcal{O}\left(\chi^{3}\right) \end{split} \tag{A16}$$

where we use L ˜(0) = 0 as follows from the definition. Similarly we obtain for R [*<sup>λ</sup>*0(*χ*)] (with R(0) ≡ R)

$$\mathcal{R}\left[\lambda\_0(\chi)\right] = \mathcal{R} + \mathcal{R}^{(1,0,0,0)}(i\chi\_1) + \mathcal{R}^{(0,1,0,0)}(i\chi\_2) + \mathcal{R}^{(0,0,1,0)}(i\chi\_3) + \mathcal{R}^{(0,0,0,1)}(i\chi\_4) + \mathcal{O}\left(\chi^2\right) \tag{A17}$$

where

$$\tilde{\mathcal{L}}^{(p,q,r,s)} = \left. \partial\_{i\chi\_1}^p \partial\_{i\chi\_2}^q \partial\_{i\chi\_3}^r \partial\_{i\chi\_4}^s \tilde{\mathcal{L}}(\chi) \right|\_{\chi=0}, \qquad \mathcal{R}^{(p,q,r,s)} = \left. \partial\_{i\chi\_1}^p \partial\_{i\chi\_2}^q \partial\_{i\chi\_3}^r \partial\_{i\chi\_4}^s \mathcal{R} \left[ \lambda o(\chi) \right] \right|\_{\chi=0} \tag{A18}$$

From the definition in Equation (A1), it follows that F({0}; *t*) = 0, so that *<sup>λ</sup>*0(0) = 0, and therefore R = QL−<sup>1</sup>Q. This operator satisfies LR = RL, RLR = *R*, and LRL = *L* and is therefore called the *Drazin pseudoinverse*. It can be shown that for a matrix of rank 4 it can be computed as

$$\mathcal{R} = (a\_3)^{-2} \mathcal{L} B\_2^2 \tag{A19}$$

with the rules

$$\mathcal{L}\_0 = 1 \qquad \qquad a\_0 = 1 \qquad \qquad \qquad B\_0 = I\_4 \tag{A20}$$

$$\mathcal{L}\_1 = \mathcal{L}B\_0 \qquad a\_1 = -\text{Tr}(\mathcal{L}\_1)/1 \qquad B\_1 = \mathcal{L}\_1 + a\_1 I\_4 \tag{A21}$$

$$\mathcal{L}\_2 = \mathcal{L}B\_1 \qquad a\_2 = -\text{Tr}(\mathcal{L}\_2)/2 \qquad B\_1 = \mathcal{L}\_2 + a\_2 I\_4 \tag{A22}$$

$$
\mathcal{L}\_3 = \mathcal{L}B\_2 \qquad a\_3 = -\text{Tr}(\mathcal{L}\_3)/3 \tag{A23}
$$

Furthermore, it is possible to prove that the first derivatives of R, which are required for the computation of the third-order cumulants, satisfy

$$\mathcal{R}^{(1,0,0,0)} = \mathcal{R}^2 \left< 0 \right| \mathcal{L}^{(1,0,0,0)} \left| 0 \right>\tag{A24}$$

$$\mathcal{R}^{(0,1,0,0)} = \mathcal{R}^2 \left< 0 \right| \mathcal{L}^{(0,1,0,0)} |0\rangle \tag{A25}$$

$$\mathcal{R}^{(0,0,1,0)} = \mathcal{R}^2 \left< 0 \right| \mathcal{L}^{(0,0,1,0)} \left| 0 \right> \tag{A26}$$

$$\mathcal{R}^{(0,0,0,1)} = \mathcal{R}^2 \left< 0 \right| \mathcal{L}^{(0,0,0,1)} \left| 0 \right>\tag{A27}$$

Finally, we find *<sup>λ</sup>*0(*χ*) in the form of a power series,

$$\lambda\_0(\chi) = \langle 0 \vert \, \bar{\mathcal{L}}(\chi) \left[ 1 - \mathcal{R}\bar{\mathcal{L}}(\chi) + (\mathcal{R}\bar{\mathcal{L}}(\chi))^2 - \dots \right] \vert 0 \rangle \tag{A28}$$

or, explicitly,

*<sup>λ</sup>*0(*χ*) = 0| #L˜(1,0,0,0)(*<sup>i</sup>χ*1) + <sup>L</sup>˜(0,1,0,0)(*<sup>i</sup>χ*2) + <sup>L</sup>˜(0,0,1,0)(*<sup>i</sup>χ*3) + <sup>L</sup>˜(0,0,0,1)(*<sup>i</sup>χ*4)+ + 1 2! L˜(2,0,0,0)(*<sup>i</sup>χ*1)<sup>2</sup> + <sup>L</sup>˜(0,2,0,0)(*<sup>i</sup>χ*2)<sup>2</sup> + <sup>L</sup>˜(0,0,2,0)(*<sup>i</sup>χ*3)<sup>2</sup> + <sup>L</sup>˜(0,0,0,2)(*<sup>i</sup>χ*4)<sup>2</sup><sup>+</sup> + 2L ˜(1,1,0,0)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*2) + <sup>2</sup>L˜(1,0,1,0)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*3) + <sup>2</sup>L˜(1,0,0,1)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*4)+ + 2L ˜(0,1,1,0)(*<sup>i</sup>χ*2)(*<sup>i</sup>χ*3) + <sup>2</sup>L˜(0,1,0,1)(*<sup>i</sup>χ*2)(*<sup>i</sup>χ*4) + <sup>2</sup>L˜(0,0,1,1)(*<sup>i</sup>χ*3)(*<sup>i</sup>χ*4)− − 2L ˜(1,0,0,0)RL˜(1,0,0,0)(*<sup>i</sup>χ*1)<sup>2</sup> − <sup>2</sup>L˜(1,0,0,0)RL˜(0,1,0,0)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*2)− − 2L ˜(1,0,0,0)RL˜(0,0,1,0)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*3) − <sup>2</sup>L˜(1,0,0,0)RL˜(0,0,0,1)(*<sup>i</sup>χ*1)(*<sup>i</sup>χ*4)− − 2L ˜(0,1,0,0)RL˜(1,0,0,0)(*<sup>i</sup>χ*2)(*<sup>i</sup>χ*1) − <sup>2</sup>L˜(0,1,0,0)RL˜(0,1,0,0)(*<sup>i</sup>χ*2)<sup>2</sup>− − 2L ˜(0,1,0,0)RL˜(0,0,1,0)(*<sup>i</sup>χ*2)(*<sup>i</sup>χ*3) − <sup>2</sup>L˜(0,1,0,0)RL˜(0,0,0,1)(*<sup>i</sup>χ*2)(*<sup>i</sup>χ*4)− − 2L ˜(0,0,1,0)RL˜(1,0,0,0)(*<sup>i</sup>χ*3)(*<sup>i</sup>χ*1) − <sup>2</sup>L˜(0,0,1,0)RL˜(0,1,0,0)(*<sup>i</sup>χ*3)(*<sup>i</sup>χ*2)− − 2L ˜(0,0,1,0)RL˜(0,0,1,0)(*<sup>i</sup>χ*3)<sup>2</sup> − <sup>2</sup>L˜(0,0,1,0)RL˜(0,0,0,1)(*<sup>i</sup>χ*3)(*<sup>i</sup>χ*4)−

$$\begin{aligned} &-2\mathcal{L}^{(0,0,0,1)}\mathcal{R}\mathcal{L}^{(1,0,0,0)}(i\chi\_4)(i\chi\_1) - 2\mathcal{L}^{(0,0,0,1)}\mathcal{R}\mathcal{L}^{(0,1,0,0)}(i\chi\_4)(i\chi\_2) - \\ &-2\mathcal{L}^{(0,0,0,1)}\mathcal{R}\mathcal{L}^{(0,0,1,0)}(i\chi\_4)(i\chi\_3) - 2\mathcal{L}^{(0,0,0,1)}\mathcal{R}\mathcal{L}^{(0,0,0,1)}(i\chi\_4)^2\Big] + \mathcal{O}\left(\chi^3\right)\Big]|0\rangle \end{aligned} \tag{A29}$$

We can now apply Equation (A2) to compute all cumulants recursively. To first order, we ge<sup>t</sup> (except for a factor of *t*)

$$C\_{1000} = \langle 0 \vert \mathcal{L}^{(1,0,0,0)} \vert 0 \rangle \tag{A30}$$

$$\mathcal{L}\_{0100} = \langle 0 \vert \, \tilde{\mathcal{L}}^{(0,1,0,0)} \vert 0 \rangle \tag{A31}$$

$$C\_{0010} = \langle 0 \vert \mathcal{L}^{(0,0,1,0)} \vert 0 \rangle \tag{A32}$$

$$C\_{0001} = \langle 0 \vert \, \tilde{\mathcal{L}}^{(0,0,0,1)} \vert 0 \rangle \tag{A33}$$

To second order, we obtain

$$\mathcal{L}\_{2000} = \left< 0 \right| \mathcal{L}^{(2,0,0,0)} - 2\mathcal{L}^{(1,0,0,0)} \mathcal{R} \mathcal{L}^{(1,0,0,0)} \left| 0 \right> \tag{A34}$$

$$\mathcal{L}\_{0200} = \left< 0 \right| \tilde{\mathcal{L}}^{(0,2,0,0)} - 2 \tilde{\mathcal{L}}^{(0,1,0,0)} \mathcal{R} \tilde{\mathcal{L}}^{(0,1,0,0)} \left| 0 \right> \tag{A35}$$

$$\mathcal{L}\_{0020} = \left< 0 \right| \mathcal{L}^{(0,0,2,0)} - 2\mathcal{L}^{(0,0,1,0)} \mathcal{R} \mathcal{L}^{(0,0,1,0)} |0\rangle \tag{A36}$$

$$\mathcal{L}\_{0002} = \left< 0 \right| \tilde{\mathcal{L}}^{(0,0,0,2)} - 2 \tilde{\mathcal{L}}^{(0,0,0,1)} \mathcal{R} \tilde{\mathcal{L}}^{(0,0,0,1)} \left| 0 \right> \tag{A37}$$

$$\mathcal{L}\_{1100} = \langle 0 \vert \mathcal{L}^{(1,1,0,0)} - \mathcal{L}^{(1,0,0,0)} \mathcal{R} \mathcal{L}^{(0,1,0,0)} - \mathcal{L}^{(0,1,0,0)} \mathcal{R} \mathcal{L}^{(1,0,0,0)} \vert 0 \rangle \tag{A38}$$

$$\mathcal{L}\_{1010} = \langle 0 \vert \, \bar{\mathcal{L}}^{(1,0,1,0)} - \bar{\mathcal{L}}^{(1,0,0,0)} \mathcal{R} \bar{\mathcal{L}}^{(0,0,1,0)} - \bar{\mathcal{L}}^{(0,0,1,0)} \mathcal{R} \bar{\mathcal{L}}^{(1,0,0,0)} \vert 0 \rangle \tag{A39}$$

$$\mathcal{L}\_{1001} = \langle 0 \vert \mathcal{L}^{(1,0,0,1)} - \mathcal{L}^{(1,0,0,0)} \mathcal{R} \mathcal{L}^{(0,0,0,1)} - \mathcal{L}^{(0,0,0,1)} \mathcal{R} \mathcal{L}^{(1,0,0,0)} \vert 0 \rangle \tag{A40}$$

$$\mathcal{L}\_{0110} = \langle 0 \vert \,\bar{\mathcal{L}}^{(0,1,1,0)} - \bar{\mathcal{L}}^{(0,1,0,0)} \mathcal{R}\bar{\mathcal{L}}^{(0,0,1,0)} - \bar{\mathcal{L}}^{(0,0,1,0)} \mathcal{R}\bar{\mathcal{L}}^{(0,1,0,0)} \vert 0 \rangle \tag{A41}$$

$$\mathcal{L}\_{0101} = \langle 0 \vert \mathcal{L}^{(0,1,0,1)} - \mathcal{L}^{(0,1,0,0)} \mathcal{R} \mathcal{L}^{(0,0,0,1)} - \mathcal{L}^{(0,0,0,1)} \mathcal{R} \mathcal{L}^{(0,1,0,0)} \vert 0 \rangle \tag{A42}$$

$$\mathcal{L}\_{0011} = \langle 0 \vert \, \bar{\mathcal{L}}^{(0,0,1,1)} - \bar{\mathcal{L}}^{(0,0,1,0)} \mathcal{R} \bar{\mathcal{L}}^{(0,0,0,1)} - \bar{\mathcal{L}}^{(0,0,0,1)} \mathcal{R} \bar{\mathcal{L}}^{(0,0,1,0)} \vert 0 \rangle \tag{A43}$$

and the same procedure is applied to higher-order cumulants. The expressions become rather cumbersome, but the recipe is clear and easy to use with some algebra.
