*2.2. Self-Consistent Equations*

From the effective Hamiltonian Equation (7), one can derive four equations of the motion of slave-boson operators, which serve as the basic equations together with the three constraints. Then, we further apply the mean-field approximation in the statistical expectations of these equations, where all the boson operators are replaced by their respective expectation values. After a lengthy and tedious calculation employing the Langreth technique (please see the Appendix B for the details of derivation), we can obtain the self-consistent equations as follows [49–52]:

$$
\Gamma\_s \frac{\partial (z\_1 z\_2)}{\partial \varepsilon} (R + R^\*) + \lambda^1 \varepsilon + \sum\_{\sigma} \frac{\partial z\_{\sigma}}{\partial \varepsilon} (Q\_{\sigma} + Q\_{\sigma}^\*) = 0,\tag{9}
$$

$$
\Gamma\_s \frac{\partial (z\_1 z\_2)}{\partial p\_1} (R + R^\*) + (\lambda^1 - \lambda\_1^2) p\_1 + \frac{\partial z\_1}{\partial p\_1} (Q\_1 + Q\_1^\*) + \frac{\partial z\_2}{\partial p\_1} (Q\_2 + Q\_2^\*) = 0,\tag{10}
$$

$$\Gamma\_s \frac{\partial (z\_1 z\_2)}{\partial p\_2} (R + R^\*) + (\lambda^1 - \lambda\_2^2) p\_2 + \frac{\partial z\_1}{\partial p\_2} (Q\_1 + Q\_1^\*) + \frac{\partial z\_2}{\partial p\_2} (Q\_2 + Q\_2^\*) = 0,\tag{11}$$

$$
\Gamma\_s \frac{\partial (z\_1 z\_2)}{\partial d} (R + R^\*) + (\mathcal{U} + \lambda^1 - \sum\_{\sigma} \lambda\_{\sigma}^2) d + \sum\_{\sigma} \frac{\partial z\_{\sigma}}{\partial d} (Q\_{\sigma} + Q\_{\sigma}^\*) = 0,\tag{12}
$$

$$\sum\_{\sigma} |p\_{\sigma}|^2 + |e|^2 + |d|^2 - 1 = 0,\tag{13}$$

$$K\_{\sigma} - |p\_{\sigma}|^2 - |d|^2 = 0,\tag{14}$$

where

$$\mathcal{K}\_1 = \left< \mathbf{c}\_{d1}^\dagger \mathbf{c}\_{d1} \right> = \int \frac{d\omega}{2\pi i} \left< \left< \mathbf{c}\_{d1}; \mathbf{c}\_{d1}^\dagger \right> \right>^\leqslant (\omega) = \frac{1}{2\pi i} \int d\omega \, \mathbf{G}\_{d11}^{<}(\omega),\tag{15}$$

$$\mathcal{K}\_2 = \left< \mathbf{c}\_{d2}^\dagger \mathbf{c}\_{d2} \right> = \int \frac{-d\omega}{2\pi i} \left< \left< \mathbf{c}\_{d2}^\dagger \mathbf{c}\_{d2} \right> \right>^\flat (\omega) = \frac{-1}{2\pi i} \int d\omega \mathbf{G}\_{d22}^\triangleright(\omega),\tag{16}$$

$$R = \left\langle \mathbf{c}\_{d1}^{\dagger} \mathbf{c}\_{d2}^{\dagger} \right\rangle = \int \frac{d\omega}{2\pi i} \left\langle \left\langle \mathbf{c}\_{d2}^{\dagger} \mathbf{c}\_{d1}^{\dagger} \right\rangle \right\rangle^{\less} (\omega) = \frac{1}{2\pi i} \int d\omega \,\mathbf{G}\_{d21}^{<}(\omega),\tag{17}$$

$$\begin{split} \mathcal{Q}\_{1\eta} &= \mathbf{z}\_{1} \Gamma\_{\eta 1} \int \frac{d\omega}{2\pi} \left\{ -\frac{i}{2} \left[ \widetilde{\Gamma}\_{L1} f\_{L}(\omega) + \widetilde{\Gamma}\_{R1} f\_{R}(\omega) \right] |\mathcal{G}\_{d11}^{R}(\omega)|^{2} \\ &- \frac{i}{2} \left[ \widetilde{\Gamma}\_{L2} (1 - f\_{L}(-\omega)) + \widetilde{\Gamma}\_{R1} (1 - f\_{R}(-\omega)) \right] |\mathcal{G}\_{d21}^{R}(\omega)|^{2} + f\_{\eta}(\omega) \mathcal{G}\_{d11}^{A}(\omega) \right\}, \end{split} \tag{18}$$

$$\begin{split} Q\_{2\eta} &= z\_2 \Gamma\_{\eta 2} \int \frac{d\omega}{2\pi} \left\{ \frac{i}{2} \left[ \widetilde{\Gamma}\_{L1} (1 - f\_L(\omega)) + \widetilde{\Gamma}\_{R1} (1 - f\_R(\omega)) \right] |G\_{d21}^R(\omega)|^2 \\ &+ \frac{i}{2} \left[ \widetilde{\Gamma}\_{L2} f\_L(-\omega) + \widetilde{\Gamma}\_{R1} f\_R(-\omega) \right] |G\_{d22}^R(\omega)|^2 - f\_{\eta}(-\omega) G\_{d22}^A(\omega) \right\}, \end{split} \tag{19}$$

and

$$Q\_{\sigma} = \sum\_{\eta} Q\_{\sigma \eta}.\tag{20}$$

Here, the QD Keldysh NGFs, *<sup>G</sup><sup>R</sup>*(*<sup>A</sup>*,<,<sup>&</sup>gt;) *dσσ* (*ω*), are the matrix elements of the 2 × 2 retarded (advanced and correlation) GF matrix *<sup>G</sup><sup>R</sup>*(*<sup>A</sup>*,<,<sup>&</sup>gt;) *d* (*ω*) = *φ*; *φ*†*<sup>R</sup>*(*<sup>A</sup>*,<,<sup>&</sup>gt;) defined in the Nambu presentation, in which the mixture Fermion operator, *φ* = (*cd*1, *<sup>c</sup>*†*d*2)*<sup>T</sup>*, has to be introduced to describe the electronic

dynamics due to the superconducting proximity effect. For the effective noninteracting Hamiltonian, the retarded and advanced GFs *G<sup>R</sup>*(*A*) *d*can be easily written in the frequency domain as

$$\left(G\_{d}^{\mathbb{R}(A)}(\omega)\right)^{-1} = \begin{bmatrix} \omega - \mathfrak{e}\_{d} - \lambda\_{1}^{2} \pm \frac{i}{2} (\widetilde{\Gamma}\_{L1} + \widetilde{\Gamma}\_{\mathbb{R}1}) & -\Gamma\_{s} z\_{1} z\_{2} \\ -\Gamma\_{s} z\_{1}^{\*} z\_{2}^{\*} & \omega + \mathfrak{e}\_{d} + \lambda\_{2}^{2} \pm \frac{i}{2} (\widetilde{\Gamma}\_{L2} + \widetilde{\Gamma}\_{\mathbb{R}2}) \end{bmatrix},\tag{21}$$

with the renormalized parameters, Γ ; *ησ* = |*<sup>z</sup>σ*|<sup>2</sup>Γ*ησ*. In addition, the correlation GFs *G*<sup>&</sup>lt;(>) *d* (*ω*) can be obtained with the help of the following Keldysh relation typical for a noninteracting system:

$$\mathcal{G}\_d^{<\langle>\rangle}(\omega) = \mathcal{G}\_d^R(\omega) \left[ \Sigma\_L^{<\langle>\rangle}(\omega) + \Sigma\_R^{<\langle>\rangle}(\omega) \right] \mathcal{G}\_d^A(\omega),\tag{22}$$

with the self-energies

$$\Sigma\_{\eta}^{<}(\omega) = i \begin{bmatrix} \check{\Gamma}\_{\eta 1}^{<} f\_{\eta}(\omega) & 0\\ 0 & \check{\Gamma}\_{\eta 2}^{<} [1 - f\_{\eta}(-\omega)] \end{bmatrix} \,' \tag{23}$$

and

$$\Sigma\_{\eta}^{>}(\omega) = -i \begin{bmatrix} \widetilde{\Gamma}\_{\eta 1} [1 - f\_{\eta}(\omega)] & 0 \\ 0 & \widetilde{\Gamma}\_{\eta 2} f\_{\eta}(-\omega) \end{bmatrix} \, , \tag{24}$$

where *fη*(*ω*) = 1/(*eβ*(*<sup>ω</sup>*−*μη* ) + 1) is the Fermi distribution function of the lead *η* with the chemical potential *μη* and temperature 1/*β*.

### *2.3. The Current and Linear Conductance*

The electric current flowing from the lead *η* into the QD can be obtained from the rate of change of the electron number operator of the left lead:

$$I\_{\eta} = \sum\_{\sigma} I\_{\eta \sigma} = -e \sum\_{\sigma} \left\langle \frac{d}{dt} \sum\_{k} c\_{\eta k \sigma}^{\dagger} c\_{\eta k \sigma} \right\rangle. \tag{25}$$

After standard calculation, the current for the left lead can be written as [44,45]

$$I\_L = I\_L^{ET} + I\_L^{DAR} + I\_L^{CAR},\tag{26}$$

with

$$I\_L^{\rm ET} = \frac{\varepsilon}{\hbar} \int d\omega \left\{ \check{\Gamma}\_{L1} \check{\Gamma}\_{R1} \left[ f\_L(\omega) - f\_R(\omega) \right] |\mathcal{G}\_{d11}^R(\omega)|^2 + \check{\Gamma}\_{L2} \check{\Gamma}\_{R2} \left[ f\_L(-\omega) - f\_R(-\omega) \right] |\mathcal{G}\_{d22}^R(\omega)|^2 \right\}, \tag{27}$$

$$I\_{L}^{DAR} = \frac{2\varepsilon}{\hbar} \int d\omega \,\widetilde{\Gamma}\_{L1} \widetilde{\Gamma}\_{L2} \left[ f\_L(\omega) + f\_L(-\omega) - 1 \right] \times |G\_{d12}^{R}(\omega)|^2,\tag{28}$$

$$d\_L^{\rm CAR} = \frac{\varepsilon}{h} \int d\omega \left\{ \tilde{\Gamma}\_{L1} \tilde{\Gamma}\_{R2} \left[ f\_L(\omega) + f\_R(-\omega) - 1 \right] + \tilde{\Gamma}\_{L2} \tilde{\Gamma}\_{R1} \left[ f\_L(-\omega) + f\_R(\omega) - 1 \right] \right\} |G\_{d12}^{\rm R}(\omega)|^2. \tag{29}$$

The corresponding currents for the right lead can be readily obtained by simply exchanging the subscripts L and R in Equations (27)–(29). It is found that the current can be divided into three parts: *IETL* describes the single-particle tunneling current caused by the normal electron transfer (ET) processes from the left lead directly to the right lead; *IDAR L* denotes the local Andreev current caused by the direct AR (DAR) processes in which an electron injecting from the left lead forms a Cooper pair in the superconducting lead, and at the same time, is reflected as a hole back into the left lead; and *ICAR L* is the nonlocal Andreev current caused by the crossed AR (CAR) processes, which is similar to DAR except that the hole is reflected into another lead, i.e., here, the right lead.

Since we are interested in the interplay between the Andreev bound state and the Kondo effect in the nonlocal subgap tunneling, we choose the bias voltage configuration in this hybrid three-terminal

nanodevice as follows: The left lead is biased with the chemical potential *V*, while the right lead and the superconducting electrode are both in contact with the ground. Therefore, one can define two different linear conductances: The usual local conductance *GL* = *∂IL*/*∂V*|*<sup>V</sup>*=<sup>0</sup> and the unusual nonlocal (cross) conductance *GC* = *∂IR*/*∂V*|*<sup>V</sup>*=0, which is related to the nonlocal current response of the hybrid three-terminal nanodevice to external driving field, i.e., current flowing in the right lead caused by the bias voltage applied to the left lead. From Equations (27)–(29), the local conductance reads

$$G\_L = \left. \frac{\partial I\_L}{\partial V} \right|\_{V=0} = G^{ET} + G^{DAR} + G^{CAR}, \tag{30}$$

and the cross conductance is

$$\mathbf{G}\_{\mathbf{C}} = \left. \frac{\partial I\_{\mathbf{R}}}{\partial V} \right|\_{V=0} = \mathbf{G}^{ET} - \mathbf{G}^{CAR},\tag{31}$$

where

$$\mathbf{G}^{ET} = \frac{e^2}{h} \left( \widetilde{\Gamma}\_{L1} \widetilde{\Gamma}\_{R1} |\mathbf{G}\_{d11}^R(0)|^2 + \widetilde{\Gamma}\_{L2} \widetilde{\Gamma}\_{R2} |\mathbf{G}\_{d22}^R(0)|^2 \right),\tag{32}$$

$$G^{DAR} = \frac{4e^2}{h} \tilde{\Gamma}\_{L1}^{\epsilon} \tilde{\Gamma}\_{L2}^{\epsilon} |G^R\_{d12}(0)|^2,\tag{33}$$

$$\mathcal{G}^{\mathcal{C}AR} = \frac{e^2}{h} \left( \widetilde{\Gamma}\_{L1} \widetilde{\Gamma}\_{R2} + \widetilde{\Gamma}\_{R1} \widetilde{\Gamma}\_{L2} \right) |G^R\_{d12}(0)|^2. \tag{34}$$

It is obvious that all of the three different tunneling processes contribute to the local conductance. Nevertheless, the DAR tunneling process, as expected, has no contribution to the cross conductance. More interestingly, the CAR tunneling process provides a contrary contribution, in comparison with the ET process, to the cross-conductance Equation (31), which is responsible for the negative value of the cross conductance in certain appropriate conditions, as shown in the following section. This opposite role of the CAR can be interpreted in an intuitive way: A hole entering the right lead is physically equivalent to an electron breaking into the QD from the right lead, thus resulting in an opposite current flowing in the right lead. It is important to point out that if the superconducting coupling is switched off (Γ*s* = 0), there are no DAR and CAR processes, and as a result, the cross conductance reduces to the local conductance.
