To study the tunneling current through the double quantum dots, it is very convenient to characterize this dynamic using the equation of motion for the time evolution of density matrix \(\rho \left(t\right)\) of the quantum system.

Within the Born-Markov approximation, the dynamics equation of motion can be written as follows [8];

$${\dot{\rho }}_{mn}\left(t\right)=-i⟨m|\left[{H}_{dots}^{0},\rho \right]|n⟩+\sum _{k\ne n}\left({\varGamma }_{nk}{\rho }_{kk}-{\varGamma }_{nk}{\rho }_{kk}\right){\delta }_{mn}-{{\Lambda }}_{mn}{\rho }_{mn}\left(1-{\delta }_{mn}\right)$$

(3)

\({\varGamma }_{nk}\) describes the transition rate from the state \(⟨\left.n\right|\) to the state \(\left|k>\right.\) and the non-adiabatic parameter, \({{\Lambda }}_{mn}\), whose real part is responsible for the time decay of the off-diagonal matrix elements (coherence) is,

$${ᴧ}_{mn}=\frac{1}{2}\sum _{k}\left({\varGamma }_{km}+{\varGamma }_{kn}\right)$$

4

one can calculate the transition rates using the Fermi Golden rule approximation;

$${\varGamma }_{mn}=\sum _{\alpha =L,R}{\varGamma }_{\alpha }\left\{f\left({E}_{m}-{E}_{n}-{\mu }_{\alpha }\right){\delta }_{{N}_{m},{N}_{n}+1}+\left[1-f\left({E}_{n}-{E}_{m}-{\mu }_{\alpha }\right){\delta }_{{N}_{m},{N}_{n}-1}\right]\right\}$$

(5)

\({N}_{m}\) is the number of electrons in the system when it is in state \(⟨\left.m\right|\) and \(f\left({E}_{m}-{E}_{n}-{\mu }_{\alpha }\right)\)is the Fermi distribution function at \(\alpha =L, R\) lead, and \({\delta }_{{N}_{m},{N}_{n}}\) is the delta function.

The available transitions that take place between the stats \(\text{n}=\text{1,2},\text{3,4}\) and the states \(\text{m}=\text{1,2},\text{3,4}\) reduced the solution of Eq. (3) at the steady states \({\dot{\rho }}_{mn}\left(t\right)=0\) to;

$${\rho }_{11}=\frac{{\varGamma }_{12}{\rho }_{22}+{\varGamma }_{13}{\rho }_{33}+{\varGamma }_{14}{\rho }_{44}}{x}$$

6

$${\rho }_{22}=\frac{{\varGamma }_{21}{\rho }_{11}+{\varGamma }_{23}{\rho }_{33}+{\varGamma }_{24}{\rho }_{44}}{y}$$

7

\({\rho }_{33}=\frac{{\varGamma }_{31}{\rho }_{11}+{\varGamma }_{32}{\rho }_{22}+{\varGamma }_{34}{\rho }_{44}}{z}\) (8)

\({\rho }_{44}=\frac{{\varGamma }_{41}{\rho }_{11}+{\varGamma }_{42}{\rho }_{22}+{\varGamma }_{43}{\rho }_{33}}{D}\) (9)

Where,

$$x=\left({\varGamma }_{21}+{\varGamma }_{31}+{\varGamma }_{41}\right)$$

$$y=\left({\varGamma }_{12}+{\varGamma }_{32}+{\varGamma }_{42}\right)$$

$$z=\left({{\Gamma }}_{13}+{{\Gamma }}_{23}+{{\Gamma }}_{43}\right)$$

$$D=\left({\varGamma }_{14}+{\varGamma }_{24}+{\varGamma }_{34}\right)$$

The total current that tunnels through the system \({I}_{tot}\) is considered as the summation of the spin-up current \({I}_{\uparrow }\) and the spin-down current \({I}_{\downarrow }\) as [9],

$${I}_{tot} = {I}_{\uparrow } + {I}_{\downarrow } = \varGamma \left[ \rho \right|0, \uparrow ⟩ + \rho |0, \downarrow ⟩]$$

10

by representing that \(\rho \left|0, \uparrow \right.⟩= {\rho }_{33}\) and \(\rho \left|0, \downarrow \right.⟩={\rho }_{44}\), so the total current in Eq. (10) is taken the form;

$${I}_{tot} ={\Gamma }\left({\rho }_{33}+{\rho }_{44}\right)$$

11

To calculate \({I}_{tot}\) it is convenient to reduce \({\rho }_{33}\), and \({\rho }_{44}\) in terms of \({\rho }_{22}\) such that,

$${\rho }_{33}=\left[\frac{\left(x1x6-x3x4\right)}{x3x5+x2x6}\right]{\rho }_{22}$$

12

$${\rho }_{44}=\left[\frac{\left({\varGamma }_{31}D+{\varGamma }_{41}{\varGamma }_{34}\right)\left(z{\varGamma }_{41}+{\varGamma }_{31}{\varGamma }_{43}\right)\left(x1x6-x3x4\right)}{\left(x3x5+x2x6\right)}+\frac{\left({\varGamma }_{31}{\varGamma }_{42}-{\varGamma }_{41}{\varGamma }_{32}\right)}{\left({\varGamma }_{31}D+{\varGamma }_{41}{\varGamma }_{34}\right)}\right]{\rho }_{22}$$

13

Where;

x1=\(\left(xy-{{\varGamma }_{21}\varGamma }_{12}\right)\) x4=\(\left({\varGamma }_{31}{\varGamma }_{42}-{\varGamma }_{41}{\varGamma }_{32}\right)\)

x2=\(\left(x{\varGamma }_{23}+{{\varGamma }_{21}\varGamma }_{13}\right)\) x5=\(\left(z{\varGamma }_{41}+{\varGamma }_{31}{\varGamma }_{43}\right)\) (14)

x3=\(\left({{\varGamma }_{21}\varGamma }_{14}+{x\varGamma }_{24}\right)\) x6=\(\left({\varGamma }_{31}D+{\varGamma }_{41}{\varGamma }_{34}\right)\)

By subtitling the values of \({\rho }_{33}\) and \({\rho }_{44}\) defined in Eq. (12), (13) and get used to the definitions in Eq. (14), then an expression for \({I}_{tot}\) will be,

$${I}_{tot}={\Gamma }\left\{\left[\frac{\left(x1x6-x3x4\right){\rho }_{22}}{x3x5+x2x6}\right]+\left[\frac{\left({\varGamma }_{31}D+{\varGamma }_{41}{\varGamma }_{34}\right)\left(z{\varGamma }_{41}+{\varGamma }_{31}{\varGamma }_{43}\right)\left(x1x6-x3x4\right)}{\left(x3x5+x2x6\right)}+ \frac{\left({\varGamma }_{31}{\varGamma }_{42}-{\varGamma }_{41}{\varGamma }_{32}\right)}{\left({\varGamma }_{31}D+{\varGamma }_{41}{\varGamma }_{34}\right)}\right]{\rho }_{22}\right\}$$

(15)