The proposed model uses transmission spectrum (TS) to calculate the characteristics of single electron transistors. In our previous works, we have implemented earlier versions of the model to molecular devices \cite{zahir2016ee} and silicon nanocrystals \cite{mehdy2016efficient}. The current model is implemented on SET and the difference from \cite{zahir2016ee} is the implementation of two spectra for different spin as explained in \cite{mehdy2016efficient}. In this work, we extend the model presented in \cite{mehdy2016efficient} to implement the gate voltage in order to use the silicon nanocrystal as SET. Further explanation on the modeling of gate voltage can be found in the following description and in the results section. The model takes as input the transmission spectrum at equilibrium from the simulations. From the transmission spectrum, current is calculated using Landauers formula. The evolution in transmission spectrum with applied voltage is calculated by the model using self-consistent field technique. By evaluating the TS at all voltages ($V_{ds},V_g$) and using the Landauers formula, IV characteristics of the device can be calculated \cite{datta2005quantum}.

\begin{figure}[t]
\centering
{\includegraphics[width=6cm]{./figure/flow.eps}
}
\caption{The flow chart of the model: at equilibrium, the transmission spectrum is constructed using parameters from simulations, while at non-equilibrium conditions, shift is calculated by the self-consistent loop.}
\label{fig:flow}
\vspace{0.0em}
\end{figure}


We calculate the TS at equilibrium using Atlas TCAD. The TS parameters i.e. position, amplitude and width of peaks are given as input to the model. The model uses these parameters to construct the TS at equilibrium. The shape of the peaks is defined by the Lorentzian function.
\begin{equation}
D_0=k_0 \times Norm(\frac {\gamma_{1} \gamma_{2}} {(E- \epsilon)^2+ \gamma^2})
\end{equation}
Here \(D_0\) is the transmission peak defined over the range of energies \(E\), \(\epsilon\) is the position of the corresponding peak and \(\gamma\) is the sum of drain and source rate constants \(\gamma_1\),\(\gamma_2\). The width (or broadening) of the transmission peaks is defined by \(\gamma\). After normalizing the Lorentzian function, it is multiplied by the amplitude of the corresponding peak \(k_0\). Other peaks are calculated similarly and then TS is constructed by adding all the transmission peaks. The same procedure is repeated for both spins and two spectra are achieved. Although at this point both the spectra are similar, the difference arises when the shift in spectra is calculated as explained below.

When the voltage is applied, transmission spectrum shifts on the energy axis. The shift is caused by the external potential and because of the change in number of electrons. At equilibrium, the Fermi level of island is aligned with the drain and source chemical potential. When the drain to source voltage (\(V_{ds}\)) is applied the drain and source chemical potential split and the separation is so called bias window. If there is a transmission peak inside or at the edge of the bias window then the number of electrons (\(N\)) at the island change. If the transmission peak corresponds to the spin up electrons the change in number of electrons are calculated using following formula. 
\begin{equation}
N_{UP/DN}=\frac {\gamma_1 f_{1UP/DN}(\epsilon) + \gamma_2 f_{2UP/DN}(\epsilon)} {\gamma_1 + \gamma_2}
\label{equ:ne}
\end{equation}
Here \(f_{1UP/DN}(\epsilon)\) and \(f_{2UP/DN}(\epsilon)\) are the drain and source Fermi functions with respect to the corresponding spin up/down transmission peaks \((\epsilon)\) If the number of electrons increase in the island the energy of the system increases so the transmission spectrum shifts towards higher energies and vice versa. The overall shift is the result of both i.e. change in number of electrons and the applied drain to source bias. The positive gate voltage lowers the energy of the island so it shifts the transmission spectrum towards low energies. The overall shift (\(U\)) can be calculated using:
\begin{equation}
U_{UP/DN}=U_L+ U_0 \Delta N_{DN/UP}
\label{equ:u}
\end{equation}
where \(U_0\) is the single electron charging energy, \(\Delta N_{DN/UP}\) is the change in number of spin down/up electrons and \(U_L=\alpha_d V_{ds}+\alpha_g V_g\) is Laplace potential. The effect of the external voltages is modeled by the coupling coefficients \(\alpha_d,\alpha_g\). Equation \ref{equ:ne} and \ref{equ:u} are solved self-consistently to calculate the shift in transmission spectrum with respect to the applied voltage. The shift in both the spectra for spin up and spin down is calculated. The shift in the spin up transmission spectrum is calculated from the change in number of spin down electrons as spin up electron feels the potential due to spin down electron and vice versa \cite{datta2005quantum}. 

\begin{figure}[b!]
\centering
{\includegraphics[width=8cm]{./figure/ts_1.eps}
}
\caption{Transmission spectra at different \(V_g\) and \(V_{ds}\) achieved from simulations and model. The bias window is shown as the dashed lines. a) \(V_g=0.0V\), \(V_{ds}=0.0V\) b) \(V_g=1.0V\), \(V_{ds}=0.0V\) c) \(V_g=1.0V\), \(V_{ds}=1.0V\) d) \(V_g=1.5V\), \(V_{ds}=1.0V\) e) \(V_g=2.5V\), \(V_{ds}=1.0V\).}
\label{fig:ts_1}
\end{figure}

When a positive gate voltage is applied and a lowest unoccupied molecular level (LUMO) goes below the drain and source chemical potentials, the number of electrons increase on the island. This increase in number of electrons causes the rate constant (\(\gamma\)) to decrease so the width of the transmission peaks also reduces. In the current version of the model, the change in width is calculated using fitting parameters. Another effect of the increase in number of electrons on the island is to reduce the gate coupling (\(\alpha_g\)). With the increase in number of electrons on the island the amount of shift in transmission spectra due to increase in gate voltage reduces. We implement this effect in the model by linear fitting of \(\alpha_g\) with respect to the change in number of electrons \((N)\). After evaluating the transmission spectra evolved with the applied voltage we can find the current using Landauers formula.
\begin{equation}
I_{UP/DN}= \frac{q^2}{h}\int_{-\infty} ^{+\infty} T_{UP/DN}(E)[f_{1}(E)-f_{2}(E)]dE
\label{equ:i}
\end{equation}
Here \(q\) is the electron charge, \(h\) is plank's constant and \(T(E)\) is the transmission spectrum. The difference between the source and drain Fermi functions models the effect of bias window so that only the part of transmission spectrum which is inside the bias window is integrated. The total current is calculated by the sum of spin up and spin down currents.

