Derivation

We assume a Hamiltonian \(H\) of the form

\[H=H_{S}+\sum_\alpha H_{\alpha,SL}+H_{\alpha,L}\]

where

\[ H_s = \sum_n \epsilon_n |n\rangle\langle n| \]

is the (diagonalizable) system Hamiltonian,

\[ H_{\alpha,L}=\sum_{k}\epsilon_{\alpha,k}a_{\alpha,k}^\dagger a_{\alpha,k} \]

and the coupling, \(H_{\alpha,SL}\) is given as

\[ H_{\alpha,SL}=\sum_{r\in\alpha}(t_r c_{r}^\dagger a_{\alpha,k}+t^\dagger_r a_{\alpha,k}^\dagger c_{r}). \]

We assume the system to be in an initial state described by the density matrix \(\rho_0 = \sum_n p_n |n\rangle\langle n|\bigotimes_\alpha \rho_\alpha^{(eq)}\), where \(\rho_i^{(eq)}=Tr[e^{-\beta H_{i,L}}]/Z\). Using the argument presented in C. Timm's work, we assume that off-diagonal elements of the density matrix vanish, Markovianity, and we can restrict ourself to the rate equation of the type

\[ \partial_t p_n = \sum_{m}\Gamma_{m\rightarrow n}p_m - \sum_m \Gamma_{n\rightarrow m}p_n, \]

subject to the normalization \(\sum_n p_n=1\). In the tunneling limit and second order in the coupling we can obtain the rates via Fermi's golden rule

\[ \Gamma_{n' n}=\Gamma_{n\rightarrow n'}=\frac{2\pi}{\hbar}\sum_{k_f, k_i}|\langle f|\sum_\alpha H_{\alpha,SL}|i \rangle|^2 \delta[\epsilon_f-\epsilon_i]. \]

We are in the end only interested in the system probabilities which is why we sum over initial and final lead states, \(k_i, k_f\). For the matrix element, we consider specified initial and final lead states which we only sum over after.

For the matrix element we use \(|i\rangle = |n,\vec{k}\rangle\) and \(|n',\vec{k}'\rangle\), where \(\vec{k}\) indicates the fact that there are \(n_L\) different momenta, and calculate

$$ \Gamma_{n',n}=\frac{2\pi}{\hbar}\sum_{k,k'}\langle n,\vec{k}|\sum_{\alpha,k''}(t_\alpha c_{\alpha}^\dagger a_{\alpha, k''}+t_\alpha^\dagger a_{\alpha,k''}^\dagger c_{\alpha})|n',\vec{k}'\rangle\langle n',\vec{k}'| \sum_{\beta,k'''}(t_\beta c_{\beta}^\dagger a_{\beta,k'''}+t_\beta^\dagger a^\dagger_{\beta,k'''}c_{\beta}) |n,\vec{k}\rangle \delta \left[\epsilon_{n}'-\epsilon_n+\epsilon_{\vec{k}}-\epsilon_{\vec{k}'}\right]\ $$ We use \(T_{r,n'n}^+ = t_r \langle n'|c_{r}^\dagger |n \rangle\), and \(T_{r,n'n}^-=t^\dagger_r \langle n'|c_{r} |n \rangle\) and expand

\[ \Gamma_{n'n}=\frac{2\pi}{\hbar}\sum_{k,k',k'',k'''}\sum_{\alpha,\beta}\bigg(T^+_{\alpha,n,n'}T^-_{\beta,n'n}\langle \vec{k} | a_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a^\dagger_{\beta,k'''} |\vec{k}\rangle \\ T^+_{\alpha,n,n'}T^+_{\beta,n'n}\langle \vec{k} | a_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a_{\beta,k'''} |\vec{k}\rangle \\ T^-_{\alpha,n,n'}T^-_{\beta,n'n}\langle \vec{k} | a^\dagger_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a^\dagger_{\beta,k'''} |\vec{k}\rangle \\ T^-_{\alpha,n,n'}T^+_{\beta,n'n}\langle \vec{k} | a^\dagger_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a_{\beta,k'''} |\vec{k}\rangle \bigg) \delta \left[\epsilon_{n'}-\epsilon_n+\epsilon_{\vec{k}'}-\epsilon_{\vec{k}}\right] $$ We see that we can discard the second and third terms. For the first and fourth we calculate the matrix element of lead states. $$ \sum_{k'}\langle \vec{k} | a_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a^\dagger_{\beta,k'''} |\vec{k}\rangle \delta \left[\epsilon_{n'}-\epsilon_n+\epsilon_{\vec{k}'}-\epsilon_{\vec{k}}\right] = \delta_{\alpha,\beta}\delta_{k'',k'''}\langle 1-n(\epsilon_{\beta k''}) \rangle \delta\left[\epsilon_{n'}-\epsilon_n+\epsilon_{k''}\right] \\ \sum_{k}\langle \vec{k} | a^\dagger_{\alpha,k''} | \vec{k}' \rangle\langle \vec{k}' | a_{\beta,k'''} |\vec{k}\rangle\delta \left[\epsilon_{n'}-\epsilon_n+\epsilon_{\vec{k}'}-\epsilon_{\vec{k}}\right]'= \delta_{\alpha,\beta}\delta_{k'',k'''}\langle n(\epsilon_{k''})\rangle \delta\left[\epsilon_{n'}-\epsilon_n-\epsilon_{k''}\right] $$ We now use $T^+_{a,n'n}=(T^-_{a,nn'})^*$, and perform the thermal average over $k$ which yields the lead chemical potential shift, and transforming $\sum_k f(\epsilon_k) \rightarrow \int f(\epsilon)\rho(\epsilon)d\epsilon$ we find $$ \Gamma_{\alpha, n'n}=\frac{2\pi}{\hbar}\int \nu_\alpha |T_{\alpha,n'n}^-|^2[1-n_F(\epsilon - \mu_\alpha)]\delta[\epsilon_{n'}-\epsilon_n+\epsilon]d\epsilon \\ + \frac{2\pi}{\hbar}\int \nu_\alpha |T_{\alpha,n'n}^+|^2n_F(\epsilon - \mu_\alpha)\delta[\epsilon_{n'}-\epsilon_n-\epsilon]d\epsilon \\ =\frac{2\pi}{\hbar}\nu_\alpha |T^-_{\alpha,n'n}|^2[1-n_F(\epsilon_n-\epsilon_{n'}-\mu_\alpha)] \\ +\frac{2\pi}{\hbar}\nu_\alpha |T^+_{\alpha,n'n}|^2n_F(\epsilon_{n'}-\epsilon_n-\mu_\alpha) \]

We find that, to second order, the rate equations do not contain direct co-tunneling terms. This needs to be corrected in the code. Furthermore we see that all Fermi functions are dependent on only a single Fermi distribution. This obsolves the currently present structure containing differences of lead potentials. These would occur in fourth order when different virtual states are present. What this demonstrate is that only tensors of the structure [n_l,n_e,n_e] are sensible.