**Appendix B**

It is easily noticed that there are seven averages in total to be determined within the SBMF approach, and the constraints provide three conditions: Equations (13) and (14). Hence, four more conditions are necessary. Originally, Kotliar and Ruckenstein [46] derived them by the saddle-point approximation for the equilibrium free energy, but this approach cannot be applied in this paper, since we are dealing with the NESS. Instead, herein, we can derive those conditions from the equations of motion of the four slave boson fields according to Equation (A3). For instance, the empty boson field *e* obeys the equation of motion:

$$i\hbar\frac{d\mathbf{e}}{dt} = \Gamma\_{\delta}\left( [\mathbf{e}, z\_1^{\dagger} z\_2^{\dagger}] \mathbf{c}\_{d1}^{\dagger} \mathbf{c}\_{d2}^{\dagger} + [\mathbf{e}, z\_1 z z\_2] \mathbf{c}\_{d2} \mathbf{c}\_{d1} \right) \\ + \sum\_{\eta kr} V\_{\eta k} \left( \mathbf{c}\_{\eta kr}^{\dagger} \mathbf{c}\_{d\sigma} [\mathbf{e}, z\_\sigma] + \mathbf{c}\_{d\sigma}^{\dagger} \mathbf{c}\_{\eta kr} [\mathbf{e}, z\_\sigma^{\dagger}] \right) \\ + \lambda^1 \mathbf{e}. \tag{A4}$$

By replacing the slave-boson operators and Lagrange multipliers to their mean values and evaluating NESS averages of the fermionic operators, one obtains, in the condition of steady state,

$$i\hbar \left\langle \frac{d\boldsymbol{\nu}}{dt} \right\rangle = \Gamma\_s \left( \frac{\partial \boldsymbol{\varepsilon}\_1^\* \boldsymbol{\varepsilon}\_2^\*}{\mathrm{d}t^\*} \langle \boldsymbol{\varepsilon}\_{d1}^\dagger \boldsymbol{\varepsilon}\_{d2}^\dagger \rangle + \frac{\partial \boldsymbol{\varepsilon}\_1 \boldsymbol{\varepsilon}\_2}{\mathrm{d}t^\*} \langle \boldsymbol{\varepsilon}\_{d2} \boldsymbol{\varepsilon}\_{d1} \rangle \right) + \sum\_{\eta \neq \nu} V\_{\eta k} \left( \frac{\partial \boldsymbol{\varepsilon}\_\sigma}{\mathrm{d}t^\*} \langle \boldsymbol{\varepsilon}\_{\eta k \sigma}^\dagger \boldsymbol{\varepsilon}\_{d\sigma} \rangle + \frac{\partial \boldsymbol{\varepsilon}\_r^\*}{\mathrm{d}t^\*} \langle \boldsymbol{\varepsilon}\_{d\sigma}^\dagger \boldsymbol{\varepsilon}\_{\eta k \sigma} \rangle \right) + \lambda^1 \boldsymbol{\varepsilon} = 0. \tag{A5}$$

Here the commutators of the boson fields are evaluated as [*e*, *<sup>z</sup>σ*] = *∂<sup>z</sup>σ*/*∂e*† and [*e*, *<sup>z</sup>*1*z*2] = *<sup>∂</sup><sup>z</sup>*1*<sup>z</sup>*2/*∂e*†, according to the Dirac's quantum algebraic scheme.

Now, we turn to discuss the averages of fermionic operators which are evaluated at a NESS with respect to the fermionic part of the effective Hamiltonian Equation (7). Since we discuss quantum transport through a nanodevice in this paper, a NESS can be defined as being when all electrodes are

at equilibrium with their own temperatures 1/*βη* and chemical potentials *μη*, and thus the statistical operator can be introduced as

$$\varrho\_0 = \frac{1}{\Xi} e^{-\beta\_L \left(H\_L - \mu\_L N\_L\right)} e^{-\beta\_R \left(H\_R - \mu\_R N\_R\right)}.\tag{A6}$$

In the above, Ξ is the normalization constant, *<sup>H</sup>η* and *<sup>N</sup>η* are the Hamiltonian and the total number of particles in the lead *η*, respectively. According to the nonequilibrium statistical theory, the average of an observable *A* ˆ with respect to such a NESS can then be evaluated as *A*<sup>ˆ</sup> = *Tr*(*A*<sup>ˆ</sup>0). Thereby, we can evaluate the averages of fermionic operators with the help of the NGF technique, for example, the operators of occupation number on the dot, *K*1 and *K*2 (Equations (15) and (16)), and the order parameter *R* in Equation (17). Moreover, we can define

$$Q\_{1\eta} = \sum\_{k} V\_{\eta k} \left< c\_{d1}^{\dagger} c\_{\eta k1} \right> = \sum\_{k} V\_{\eta k} \int \frac{d\omega}{2\pi i} \left< \left< c\_{\eta k1}; c\_{d1}^{\dagger} \right> \right>^{<} (\omega) = \sum\_{k} V\_{\eta k} \int \frac{d\omega}{2\pi i} G\_{\eta k, d; 11}^{<}(\omega), \tag{A7}$$

and

$$Q\_{2\eta} = \sum\_{k} V\_{\eta k} \left< \mathfrak{c}\_{\eta k2}^{\dagger} \mathfrak{c}\_{d2} \right> = \sum\_{k} V\_{\eta k} \int \frac{d\omega}{2\pi i} \left< \left< \mathfrak{c}\_{\eta k2}^{\dagger} ; \mathfrak{c}\_{d2} \right> \right>^{\varsigma} (\omega) = \sum\_{k} V\_{\eta k} \int \frac{d\omega}{2\pi i} \mathcal{G}\_{\eta k, d22}^{\varsigma}(\omega). \tag{A8}$$

Here, the hybrid NGFs, *<sup>G</sup>*<sup>&</sup>lt;*ηk*,*<sup>d</sup>*(*ω*) are the matrix elements of the 2 × 2 hybrid contour-order GF matrix *<sup>G</sup>γη<sup>k</sup>*,*<sup>d</sup>*(*ω*) = L*ηk*; *φ*†*γ* (*γ* = *R*, *A*, <, >) with the Nambu representation in the lead <sup>L</sup>*η<sup>k</sup>* = (*<sup>c</sup>ηk*1, *<sup>c</sup>*†*ηk*2)*<sup>T</sup>*. Applying the Langreth theorem, we obtain

$$\mathcal{G}^{<'}\_{\eta k,d}(\omega) = \mathcal{g}^{\mathbb{R}}\_{\eta k}(\omega)\hat{\mathcal{V}}\_{\eta k}\mathcal{G}^{<'}\_{d}(\omega) + \mathcal{g}^{<'}\_{\eta k}(\omega)\hat{\mathcal{V}}\_{\eta k}\mathcal{G}^{A}\_{d}(\omega),\tag{A9}$$

where *gγηk* is the decoupled NGF of the lead *η* defined as *<sup>g</sup>γηk*(*ω*) = L*ηk*;L†*η<sup>k</sup><sup>γ</sup>* and

$$
\hat{V}\_{\eta k} = \begin{pmatrix} z\_1 V\_{\eta k} & 0 \\ 0 & -z\_2 V\_{\eta k}^\* \end{pmatrix} . \tag{A10}
$$

Subsequently, we have

$$\begin{split} \sum\_{k} \hat{\mathcal{V}}\_{\eta k} \mathcal{G}\_{\eta k, d}^{<}(\omega) &= \sum\_{k} \hat{\mathcal{V}}\_{\eta k} \mathcal{S}\_{\eta k}^{R}(\omega) \hat{\mathcal{V}}\_{\eta k} \mathcal{G}\_{d}^{<}(\omega) + \hat{\mathcal{V}}\_{\eta k} \mathcal{g}\_{\eta k}^{<}(\omega) \hat{\mathcal{V}}\_{\eta k} \mathcal{G}\_{d}^{A}(\omega), \\ &= \Sigma\_{\eta}^{R}(\omega) \mathcal{G}\_{d}^{<}(\omega) + \Sigma\_{\eta}^{<}(\omega) \mathcal{G}\_{d}^{A}(\omega), \end{split} \tag{A11}$$

with

$$\Sigma^{\gamma}\_{\eta}(\omega) = \sum\_{k} \hat{\mathsf{V}}\_{\eta k}^{\*} \underline{\mathbf{g}}^{\gamma}\_{\eta k}(\omega) \boldsymbol{\mathcal{V}}\_{\eta k} = \begin{pmatrix} \sum\_{k} |V\_{\eta k}|^{2} \mathbf{g}^{\gamma}\_{\eta k; \mathfrak{U}}(\omega) & \mathbf{0} \\ \mathbf{0} & \sum\_{k} |V\_{\eta k}|^{2} \mathbf{g}^{\gamma}\_{\eta k; \mathfrak{U} \mathbf{2}}(\omega) \end{pmatrix}. \tag{A12}$$

Noticing the retarded and advanced self energies in the wide band limit,

$$
\Sigma^{\mathbb{R}(A)}\_{\eta} = \pm \frac{i}{2} \begin{pmatrix} \widetilde{\Gamma}^{\mathbb{s}}\_{\eta 1} & 0 \\ 0 & \widetilde{\Gamma}\_{\eta 2} \end{pmatrix} . \tag{A13}
$$

and Equations (21)–(24), we can further obtain an explicit expression of Equation (A11) and insert this expression into Equations (A7) and (A8) to yield Equations (18) and (19). Finally, we obtain the self-consistent equation Equation (9).

Similarly, from the equations of motion of *pσ* and *d*,

$$\begin{split} i\hbar \frac{dp\_{\sigma}}{dt} = \Gamma\_{s} \left( [p\_{\sigma}, z\_{1}^{\dagger} z\_{2}^{\dagger}] c\_{d1}^{\dagger} c\_{d2}^{\dagger} + [p\_{\sigma}, z\_{1} z\_{2}] c\_{d2} c\_{d1} \right) + \sum\_{\eta k\sigma'} V\_{\eta k} \left( c\_{\eta k\sigma'}^{\dagger} c\_{d\sigma'} [p\_{\sigma}, z\_{\sigma}^{\dagger}] + c\_{d\sigma'}^{\dagger} c\_{\eta k\sigma'} [p\_{\sigma}, z\_{\sigma}^{\dagger}] \right) \\ + (\lambda^{1} - \lambda\_{\sigma}^{2}) p\_{\sigma\prime} \\ i\hbar \frac{d d}{dt} = \Gamma\_{s} \left( [d, z\_{1}^{\dagger} z\_{2}^{\dagger}] c\_{d1}^{\dagger} c\_{d2}^{\dagger} + [d, z\_{1} z\_{2}] c\_{d2} c\_{d1} \right) + \sum\_{\eta k\sigma} V\_{\eta k} \left( c\_{\eta k\sigma}^{\dagger} c\_{d\sigma} [d, z\_{\sigma}] + c\_{d\sigma}^{\dagger} c\_{\eta k\sigma} [d, z\_{\sigma}^{\dagger}] \right) \\ + (\lambda I + \lambda^{1} - \Sigma\_{\sigma} \lambda\_{\sigma}^{2}) d\_{\tau} \end{split} \tag{A14}$$

we can ge<sup>t</sup> self-consistent equations, Equations (10)–(12).
