Theoretical background

We follow the derivation of the rate equation as described in Phys. Rev. B 106, L201404 (2022).

We consider a system described by the Hamiltonian:

$$H = \sum_i E_i |E_i\rangle \langle E_i| + \sum_{r,k,\sigma} \epsilon_{rk\sigma} c_{rk\sigma}^\dagger c_{rk\sigma} + H_{SL},$$

where $E_i$ are the eigenenergies of the isolated system, $\epsilon_{rk\sigma}$ are the energies of the lead states, $c_{rk\sigma}^\dagger$ and $c_{rk\sigma}$ are the creation and annihilation operators for lead states, and $H_{SL}$ is the coupling between the system and the leads.

In this basis, the coupling between the system and the leads is given by:

$$H_{SL} = \sum_{i,rk\sigma} T_{ijr\sigma} c_{rk\sigma}^\dagger |E_i\rangle \langle E_j| + \text{h.c.},$$

where $T_{irk\sigma}$ are the coupling matrix elements between the system states and the lead states. We look in more detail at the couplings.

Coupling to Leads

The transition matrix element between dot state $|E_i\rangle$ and $|E_j\rangle$ coupled to lead $r$ is

$$ T_{ijr\sigma} = \sum_{j} \Gamma_{ijr\sigma} \langle E_i | d^\dagger_{j\sigma} | E_j \rangle, $$

where $\Gamma_{ijr\sigma}$ are the coupling strengths between the lead states and the system states, and $d^\dagger_{j\sigma}$ are the creation operators for the system states.

Note: The index $j$ runs only over the normal sites of the system.

Transition Rates

We asumme that $2\pi\rho_r = 1$ so that $T_{ijr\sigma}$ has the correct units. Then we build the transition rates as $$ W_{ij\sigma+} = \sum_r T_{ijr\sigma} f_{r}(E_i - E_j), $$ $$ W_{ij\sigma-} = \sum_r T_{ijr\sigma} (1 - f_{r}(E_i-E_j)), $$ where $f_{r\sigma}(\epsilon_{ij})$ is the Fermi-Dirac distribution for lead $r$.

Rate Equation

Finally, we can write the rate equation for the occupation of the system states $n_i$ as $$ \frac{d}{dt} P_i = \sum_{j,\sigma} \left(W_{ij\sigma+} + W_{ij\sigma-}\right) P_j - \sum_{j\sigma} \left(W_{ji\sigma+} + W_{ji\sigma-} \right) P_i, $$ where $P_i$ is the probability of finding the system in state $|E_i\rangle$. We solve for the steady state by setting $\frac{d}{dt} P_i = 0$. Furthermore, we impose the normalization condition $\sum_i P_i = 1$.