\documentclass{article}

\usepackage{graphicx}
\usepackage{setspace}
\usepackage{listings}
\usepackage{color}
\usepackage{amsmath}
\usepackage{float}
\usepackage[justification=centering]{caption}
\usepackage{subcaption}
\usepackage[margin=0.75in]{geometry}

\renewcommand{\thesubsection}{\indent(\alph{subsection})}

\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{mauve}{rgb}{0.58,0,0.82}


\definecolor{mygreen}{RGB}{28,172,0} % color values Red, Green, Blue
\definecolor{mylilas}{RGB}{170,55,241}


\lstset{language=Matlab,%
    %basicstyle=\color{red},
    breaklines=true,%
    morekeywords={matlab2tikz},
    keywordstyle=\color{blue},%
    morekeywords=[2]{1}, keywordstyle=[2]{\color{black}},
    identifierstyle=\color{black},%
    stringstyle=\color{mylilas},
    commentstyle=\color{mygreen},%
    showstringspaces=false,%without this there will be a symbol in the places where there is a space
    numbers=left,%
    numberstyle={\tiny \color{black}},% size of th_sume numbers
    numbersep=9pt, % this defines how far the numbers are from the text
    emph=[1]{for,end,break},emphstyle=[1]\color{red}, %some words to emphasise
    %emph=[2]{word1,word2}, emphstyle=[2]{style},    
    xleftmargin=5em
}
%\lstset{basicstyle=\small,
%        keywordstyle=\color{mauve},
%        identifierstyle=\color{dkgreen},
%        stringstyle=\color{gray},
%        numbers=left,
%        xleftmargin=5em
%        }

\title{ECE 456 - Problem Set 1}
\date{2021-02-06}
\author{David Lenfesty \\ lenfesty@ualberta.ca
    \and Phillip Kirwin \\ pkirwin@ualberta.ca}

\setcounter{tocdepth}{2} % Show subsections

\begin{document}

    \doublespacing
    \pagenumbering{gobble}
    \maketitle
    \newpage

    \singlespacing
    \pagenumbering{arabic}

    \section*{Question 1}
        \subsection*{(a)}
        Beginning with the following two equations:
        \begin{equation}
            \label{eq:N_old}
            N = \int_{-\infty}^{\infty}\frac{\gamma_1 f_1(E) + \gamma_2 f_2(E)}{\gamma_1 + \gamma_2} D(E-U) dE,
        \end{equation}

        \begin{equation}
            \label{eq:I_old}
            I = \frac{q}{\hbar}\frac{\gamma_1 \gamma_2}{\gamma_1 + \gamma_2} \int_{-\infty}^{\infty}[f_1(E) - f_2(E)] D(E-U) dE,
        \end{equation}
        
        and changing the variable of integration to $E' = E - U$:
        
        \begin{equation*}
            N = \int_{-\infty}^{\infty}\frac{\gamma_1 f_1(E' + U) + \gamma_2 f_2(E' + U)}{\gamma_1 + \gamma_2} D(E') dE',
        \end{equation*}

        \begin{equation*}
            I = \frac{q}{\hbar}\frac{\gamma_1 \gamma_2}{\gamma_1 + \gamma_2} \int_{-\infty}^{\infty}[f_1(E' + U) - f_2(E' + U)] D(E') dE'.
        \end{equation*}
        
        Replacing $E' \rightarrow E$, we obtain equations \ref{eq:N_new} and \ref{eq:I_new}.

        \begin{equation}
            \label{eq:N_new}
            N = \int_{-\infty}^{\infty}\frac{\gamma_1 f_1(E + U) + \gamma_2 f_2(E + U)}{\gamma_1 + \gamma_2} DE dE,
        \end{equation}

        \begin{equation}
            \label{eq:I_new}
            I = \frac{q}{\hbar}\frac{\gamma_1 \gamma_2}{\gamma_1 + \gamma_2} \int_{-\infty}^{\infty}[f_1(E + U) - f_2(E + U)] DE dE.
        \end{equation}

        \subsection*{(b)}

        For the provided constants, the plots of the number of channel electrons and the channel current follow:


        \begin{figure}[H]
            \centering
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1b_electrons.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1b_electrons}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}%
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1b_current.png}
                \caption{Plot of channel current vs. drain voltage.}
                \label{fig:q1b_current}
            \end{subfigure}

            \caption{Number of electrons and current versus drain voltage.}
        \end{figure}

            Below is our code. Note that some variable names are different from those in the example code.
            \begin{lstlisting}[language=Matlab]
            
clear all;

%% Constants

% Physical constants
hbar = 1.052e-34;

% Single-charge coupling energy (eV)
U_0 = 0.25;
% (eV)
kBT = 0.025;
% Contact coupling coefficients (eV)
gamma_1 = 0.005;
gamma_2 = gamma_1;
gamma_sum = gamma_1 + gamma_2;
% Capacitive gate coefficient
a_G = 0.5;
% Capacitive drain coefficient
a_D = 0.5;
a_S = 1 - a_G - a_D;

% Central energy level
mu = 0;

% Energy grid, from -1eV to 1eV
NE = 501;
E = linspace(-1, 1, NE);
dE = E(2) - E(1);
% TODO name this better
cal_E = 0.2;

% Lorentzian density of states, normalized so the integral is 1
D = (gamma_sum / (2*pi)) ./ ( (E-cal_E).^2 + (gamma_sum/2).^2 );
D = D ./ (dE*sum(D));

% Reference no. of electrons in channel
N_0 = 0;

voltages = linspace(0, 1, 101);

% Terminal Voltages
V_G = 0;
V_S = 0;

for n = 1:length(voltages)
    % Set varying drain voltage
    V_D = voltages(n);

    % Shifted energy levels of the contacts
    mu_1 = mu - V_S;
    mu_2 = mu - V_D;

    % Laplace potential, does not change as solution is found (eV)
    % q is factored out here, we are working in eV
    U_L = - (a_G*V_G) - (a_D*V_D) - (a_S*V_S);

    % Poisson potential must change, assume 0 initially (eV)
    U_P = 0;

    % Assume large rate of change
    dU_P = 1;

    % Run until we get close enough to the answer
    while dU_P > 1e-6
        % source Fermi function
        f_1 = 1 ./ (1 + exp((E + U_L + U_P - mu_1) ./ kBT));
        % drain Fermi function
        f_2 = 1 ./ (1 + exp((E + U_L + U_P - mu_2) ./ kBT));

        % Update channel electrons against potential
        N(n) = dE * sum( ((gamma_1/gamma_sum) .* f_1 + (gamma_2/gamma_sum) .* f_2) .* D);

        % Re-update Poisson portion of potential
        tmpU_P = U_0 * ( N(n) - N_0);
        dU_P = abs(U_P - tmpU_P);

        % Unsure why U_P is updated incrementally, perhaps to avoid oscillations?
        %U_P = tmpU_P;
        %U_P = U_P + 0.1 * (tmpU_P - U_P);
    end

    % Calculate current based on solved potential.
    % Note: f1 is dependent on changes in U but has been updated prior in the loop
    I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f_1-f_2).*D);

end


%%Plotting commands

figure(1); 
h = plot(voltages, N,'k'); 
grid on;
set(h,'linewidth',[2.0]); 
set(gca,'Fontsize',[18]);
xlabel('Drain voltage [V]'); 
ylabel('Number of electrons');

figure(2); 
h = plot(voltages, I,'k'); 
grid on;
set(h,'linewidth',[2.0]); 
set(gca,'Fontsize',[18]);
xlabel('Drain voltage [V]'); 
ylabel('Current [A]');
            \end{lstlisting}

        \subsection*{(c)}

        Code for this section can be found in appendix A.

        \subsubsection*{(i)}
        \begin{figure}[H]
            \centering
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1c_1.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1c_1}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}%
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1c_2.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1c_2}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1c_3.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1c_3}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}%
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1c_4.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1c_4}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}
            \begin{subfigure}{0.5\textwidth}
                % Mess around with widths later
                \centering
                \includegraphics[width=\textwidth]{q1c_5.png}
                \caption{Plot of channel electrons vs. drain voltage.}
                \label{fig:q1c_5}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{subfigure}%

            \caption{Visual representation of the Fermi functions of the contacts and channel.}
        \end{figure}
        
        Table \ref{table:q1ci} shows the variation in the difference between $f_1$ and $f_2$ at $E = \varepsilon$, and $I$ at different drain voltages. 
        We compare the differences between values at $V_D = 1.0\:V$ and $V_D = 0.8\:V$:

        \begin{equation*}
            \frac{([f_1(E + U) - f_2(E + U)]|_{E = \varepsilon})|_{V_D = 1.0\:V}}{([f_1(E + U) - f_2(E + U)]|_{E = \varepsilon})|_{V_D = 0.8\:V}} = 1.0363,
        \end{equation*}
        \begin{equation*}
        \frac{I|_{V_D = 1.0\:V}}{I|_{V_D = 0.8\:V}} = 1.0504,
        \end{equation*}
        %
        and between values at $V_D = 0.8\:V$ and $V_D = 0.5\:V$:
        %
        \begin{equation*}
            \frac{([f_1(E + U) - f_2(E + U)]|_{E = \varepsilon})|_{V_D = 0.8\:V}}{([f_1(E + U) - f_2(E + U)]|_{E = \varepsilon})|_{V_D = 0.5\:V}} = 2.1909,
        \end{equation*}
        \begin{equation*}
        \frac{I|_{V_D = 0.8\:V}}{I|_{V_D = 0.5\:V}} = 2.1218.
        \end{equation*}
        %
        In both comparisons we see that $I$ changes in proportion to $[f_1(E + U) - f_2(E + U)]|_{E = \varepsilon}$, as predicted by Equation (9).
        
        \begin{center}
        \begin{tabular}{l | c | c}
                
                $V_D$ [V]       & $[f_1(E + U) - f_2(E + U)]|_{E = \varepsilon}$    & I [nA]          \\
                \hline
                0.0         & 0.000                                                   & 0.0           \\
                0.2         & 0.015                                                   & 17.0          \\
                0.5         & 0.440                                                   & 271           \\
                0.8         & 0.964                                                   & 575           \\
                1.0         & 0.999                                                   & 604           \\
            
            \end{tabular}
            \captionof{table}{Differences in contact Fermi functions evaluated at $E=\varepsilon$ and current $I$ at different drain voltages $V_D$. 
            Values are taken from Figures \ref{fig:q1b_current} and \ref{fig:q1c_1} to \ref{fig:q1c_5}.}
            \label{table:q1ci}
        \end{center}


        \subsubsection*{(ii)}
            Figure \ref{fig:q1c_ii} has the Fermi functions marked at their "step points", or when they are equal to $0.5$.
            This can be used to find the self-consistant potential $U$, via the equation for the source Fermi function:

            \begin{equation*}
                f(E + U) = {\frac{1}{1 + e^{(E + U - \mu_1) / k_B T}}}
            \end{equation*}
            We could also use the drain Fermi function but since $V_S = 0$, it is simpler to use the source.
            %
            The function will "step" when $E = \mu_1 - U$, which from Figure \ref{fig:q1c_ii} occurs at $E = \mu_1 - U = 0.4$ for the source contact.
            Since $\mu_1 = \mu - q V_S = 0$ We can simply rearrange to find $U = - E = -0.4\:eV$.

            \begin{figure}[H]
                % Mess around with widths later
                \centering
                \includegraphics[width=0.7\textwidth]{q1c_ii.png}
                \caption{Marked step points of contact Fermi functions.}
                \label{fig:q1c_ii}
                % Note that the comment after \end{subfigure} is required for side by side figures
            \end{figure}%

        \subsubsection*{(iii)}

            At $V_D = 1V$, the source is trying to {\em fill} the channel level while the drain is trying to {\em empty} it.
            This is because the source has electrons at the channel level, and is filling these in while the drain does not have
            any and is attempting to bring the channel level back down to where it is.

        \subsubsection*{(iv)}

            Referring to figure \ref{fig:q1c_ii} again, we can see the areas where the difference between the Fermi functions of each contact
            are 1. Roughly, this means that the channel current $I$ would remain the same if the channel energy level was anywhere between
            $0.3eV$ and $-0.5eV$.

        \subsubsection*{(v)}

            There is no current when $V_D = 0$ because the source does not want to fill the channel. It has no electrons
            at the channel energy level, and thus there is no impetus to fill the channel. A similar story occurs with the
            drain, in that it has no electrons at the appropriate energy level, and there are none in the channel for it to pull out.

    \section*{Question 2}

        \subsection*{(a)}

        \begin{figure}[H]
            \centering
            \begin{subfigure}{0.7\textwidth}
                \centering
                \includegraphics[width=0.7\textwidth]{q2_I-V.png}
                \caption{I-V curves.}
                \label{fig:q2_iv}
            \end{subfigure}
            \begin{subfigure}{0.33\textwidth}
                \centering
                \includegraphics[width=\textwidth]{q2_0V.png}
                \caption{Fermi functions at $0V$.}
                \label{fig:q2_0v}
            \end{subfigure}%
            \begin{subfigure}{0.33\textwidth}
                \centering
                \includegraphics[width=\textwidth]{q2_0V05.png}
                \caption{Fermi functions at $0.05V$.}
                \label{fig:q2_0v}
            \end{subfigure}%
            \begin{subfigure}{0.33\textwidth}
                \centering
                \includegraphics[width=\textwidth]{q2_0V1.png}
                \caption{Fermi functions at $0.1V$.}
                \label{fig:q2_0v}
            \end{subfigure}
            \begin{subfigure}{0.33\textwidth}
                \centering
                \includegraphics[width=\textwidth]{q2_0V2.png}
                \caption{Fermi functions at $0.2V$.}
                \label{fig:q2_0v}
            \end{subfigure}%
            \begin{subfigure}{0.33\textwidth}
                \centering
                \includegraphics[width=\textwidth]{q2_0V3.png}
                \caption{Fermi functions at $0.3V$.}
                \label{fig:q2_0v}
            \end{subfigure}
            \caption{(a) depicts I-V curves at $V_G = 0.5$ V (higher curve) and $V_G = 0.25$ V. 
             (b) through (f) show Fermi functions and $D(E)$ at various drain voltages.}
        \end{figure}

        \subsection*{(b)}

        As the drain voltage increases, we see that the overlapping area
        under the "curve" of $f_1(E + U) - f_2(E+U)$ and $D(E)$ increases, up
        until the point where the drain only has electrons below the energy
        levels available in the channel. The difference function also reaches its maximum height and it widens downward rather than upward. 
        At this point, the overlapping area no longer
        increases, and therefore the current does not increase either,
        hitting {\em saturation current}.

        %% WRONG EXPLANATION - kept for posterity :)
        % We can see from these figures that the area under the $f_1(E + U) - f_2(E + U)$ curve
        % gradually moves down in energy. This means that a smaller and smaller portion of it is
        % above $E = 0eV$, which means there are fewer total energy levels through which current
        % can flow. From this we can see that with a high enough $V_D$ we will hit saturation, and
        % be unable to pass more current. \\

        \subsection*{(c)}

        At $V_D = 0.3V$, the energy levels from approximately $0$ to $0.2$ $eV$ are being
        used for electron transport. This is the range where the difference in the contact Fermi
        functions is more than 0 and the energy is more than 0. \\
        The difference between the two contacts is the greatest at $0eV$, because this is the
        point at which their Fermi functions have the greatest difference, and thus will be making
        the most "effort" to equalize the channel potential.

        \subsection*{(d)}

        Based on the supposed material changes, we would choose material A to maximize the
        drain current. The energy levels vanishing at $-0.4eV$ would mean that there is a
        greater number of energy levels that would be able to be used for conduction. There
        would be a larger {\em area} where there is a non-zero difference in the two Fermi
        functions at the contacts and there are energy levels available in the channel. \\
        This can be contrasted with material B, which would have {\em no} energy levels in
        this "conduction" zone and thus no current would be able to flow at all.

    \section*{Question 3}
    
        \subsection*{(a)}

        Below is our code. Note that some variable names are different from those in the example code.


        \begin{lstlisting}[language=Matlab]
% Thermo-electric current

% Physical constants
hbar = 1.054e-34;
q = 1.602e-19;

% Parameters (eV)
kBT_1 = 0.025;
kBT_2 = 0.026;
mu = 0;
gamma_1 = 0.005;
gamma_2 = gamma_1;
gamma_sum = gamma_1 + gamma_2;

% Channel energy levels, varying between -0.25eV and 0.25eV
epsilon = linspace(-0.25, 0.25, 101);
depsilon = epsilon(2) - epsilon(1);

% Energy grid
E = linspace(-1, 1, 501);
dE = E(2) - E(1);

% Contact fermi functions
f_1 = 1 ./ (1 + exp((E - mu)./kBT_1));
f_2 = 1 ./ (1 + exp((E - mu)./kBT_2));

% Iterate through channel energy levels
for n = 1:length(epsilon)
    % Compute energy level density functions - integral normalized to unity
    D = (gamma_sum./(2*pi))./((E-epsilon(n)).^2+((gamma_sum./2).^2));
    D = D./(dE*sum(D));
    
    % Compute number of channel electrons
    
    N(n) = dE*sum( ((gamma_1./gamma_sum).*f_1 + (gamma_2./gamma_sum).*f_2).*D );

    % Compute the current in Amps; factor of q to resolve units

    I(n) = q*(q/hbar)*dE*sum((f_1 - f_2).*D.*gamma_1.*gamma_2./gamma_sum);

    %plot f_1 - f_2 and D/2

    if (abs(epsilon(n) + 0.05) <= depsilon / 2) & (epsilon(n) <= 0)
        figure(3);
        h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-');
        set(gca, 'Fontsize', [18]);
        axis([-0.01 0.02 -1 1]);
        xlabel('f1(E) - f2(E), D(E)/2500');
        ylabel('ENERGY  [eV]');
        legend('f1-f2', 'D(E)/2500');
        title('CHANNEL LEVEL = -0.05 eV');
    elseif (abs(epsilon(n)) <= depsilon / 2) & (epsilon(n) <= 0)
        figure(4);
        h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-');
        set(gca, 'Fontsize', [18]);
        axis([-0.01 0.02 -1 1]);
        xlabel('f1(E) - f2(E), D(E)/2500');
        ylabel('ENERGY  [eV]');
        legend('f1-f2', 'D(E)/2500');
        title('CHANNEL LEVEL = 0 eV');
    end

end


% Final plots
figure(1);
        \end{lstlisting}

        \subsection*{(b)}

        \begin{figure}[H]
            \centering
            \begin{subfigure}{0.5\linewidth}
                \centering
                \includegraphics[width=\textwidth]{q3_ne.png}
                \caption{Fermi Functions at $0.3V$.}
                \label{fig:q3_ne}
            \end{subfigure}%
            \begin{subfigure}{0.5\linewidth}
                \centering
                \includegraphics[width=\textwidth]{q3_current.png}
                \caption{Fermi Functions at $0.3V$.}
                \label{fig:q3_current}
            \end{subfigure}
        \end{figure}

        \subsection*{(c)}

        \begin{figure}[H]
            \centering
            \begin{subfigure}{0.5\linewidth}
                \centering
                \includegraphics[width=\textwidth]{q3_0V05.png}
                \caption{Fermi Functions at $0.3V$.}
                \label{fig:q3_0V05}
            \end{subfigure}%
            \begin{subfigure}{0.5\linewidth}
                \centering
                \includegraphics[width=\textwidth]{q3_0V.png}
                \caption{Fermi Functions at $0.3V$.}
                \label{fig:q3_0V}
            \end{subfigure}
        \end{figure}

        \subsection*{(d)}
        
        The difference in temperature causes a difference in the sharpness of the contact Fermi functions.
        This in turn leads to the behaviour in $f_1 - f_2$ seen in Figures \ref{fig:q3_0V05} and \ref{fig:q3_0V}.
        As seen in the previous questions, the current is proportional to the area under the curve of the $D(E) * [f_1 - f_2]$.
        When the channel level is $\varepsilon = -0.05$ eV, the positive region of $[f_1 - f_2]$ overlaps the non-zero part of $D(E)$, giving a positive current.
        Alternatively, when $\varepsilon = 0$ eV, half of the area under $D(E)$ overlaps with the negative part of $[f_1 - f_2]$, while half overlaps the positive part.
        Thus the areas cancel out completely, and the resultant current is 0. A plot of $\varepsilon = +0.05$ eV would show overlap of $D(E)$ with the negative part of 
        $[f_1 - f_2]$, explaining why we get a reverse current flow at that channel level.
        The maximum current occurs at $\varepsilon = \pm0.4$ eV (relative to $\mu$).

    \section*{Question 4}

        \subsection*{(a)}

            \subsubsection*{(i)}

                \begin{equation*}
                    I = 0, N = 0
                \end{equation*}
                %
                This is because the only allowed energy level in the channel is at a higher energy
                level than exists in either of the contacts, thus there are no electrons that would
                flow into the channel from either contact, thus no current {\em and} no electrons in the channel.

            \subsubsection*{(ii)}

                \begin{equation*}
                    I = 608\:nA, N = 0.5
                \end{equation*}
                %
                Given that we are operating with a single energy level in the channel, we can use
                the equations 9 and 10 (provided in the assignment) directly.

                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005eV \cdot 0.005eV}{0.005eV + 0.005eV}
                        \cdot [1 - 0] = 608\:nA
                \end{equation*}

                \begin{equation*}
                    N = \frac{0.005eV \cdot 1 + 0.005eV \cdot 0}{0.005eV + 0.005eV} = 0.5
                \end{equation*}

            \subsubsection*{(iii)}

                \begin{equation*}
                    I = 0\:A, N = 1
                \end{equation*}
                %
                Given that we are operating with a single energy level in the channel, we can use
                the equations 9 and 10 (provided in the assignment) directly.

                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005eV \cdot 0.005eV}{0.005eV + 0.005eV}
                        \cdot [1 - 1] = 0\:A 
                \end{equation*}

                \begin{equation*}
                    N = \frac{0.005eV \cdot 1 + 0.005eV \cdot 1}{0.005eV + 0.005eV} = 1
                \end{equation*}


        \subsection*{(b)}

            \subsubsection*{(i)}

            For $f_1(E + U)$ the step point occurs at $E = \mu_1 - U = 0.25$ eV. Since $U = -0.25$ eV, it follows that 
            $\mu_1 = 0$ eV. Since $\mu_1 = \mu - qV_S$ and $\mu = 0$ eV, $V_S = 0$ V.
            For $f_2(E + U)$ the step point occurs at $E = \mu_2 - U = -0.25$ eV. Since $U = -0.25$ eV, it follows that 
            $\mu_2 = -0.5$ eV. Since $\mu_1 = \mu - qV_D$ and $\mu = 0$ eV, $V_D = 0.5$ V.

            \subsubsection*{(ii)}

            For $f_1(E)$ the step point occurs at $E = \mu_1 = 0.25$ eV. Since $\mu_1 = \mu - qV_S$ and $\mu = 0$ eV, $V_S = -0.25$ V.
            For $f_2(E)$ the step point occurs at $E = \mu_2 = -0.25$ eV. Since $\mu_1 = \mu - qV_D$ and $\mu = 0$ eV, $V_D = 0.25$ V.
            This assumes $U = 0$ eV.

        \subsection*{(c)}

            \subsubsection*{(i)}

            \begin{figure}[H]
                \centering
                \includegraphics[width=0.7\textwidth, angle=90, origin=c]{q4c.jpg}
                \caption{Visualisation of energy levels.}
                \label{fig:q4c}
            \end{figure}

            \subsubsection*{(ii)}
                Starting with equation 2 in the assignment we get:

                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot \int_{-0.1}^{0.2} [1 - 0] \cdot 10^4 dE \\
                \end{equation*}
                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot [10^4 \cdot 0.2 - 10^4 \cdot -0.1]
                \end{equation*}
                \begin{equation*}
                    I = 1.8\:mA
                \end{equation*}

            \subsubsection*{(iii)}

                The graph for the fermi functions at the saturation potentials will look much the same,
                except that they will be shifted up. This means that there is more area "under" the
                difference between the fermi functions and thus the bounds of integration can shift and
                become $[-0.1, 0.3]$.

                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot \int_{-0.1}^{0.3} [1 - 0] \cdot 10^4 dE \\
                \end{equation*}
                \begin{equation*}
                    I = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot [10^4 \cdot 0.3 - 10^4 \cdot -0.1]
                \end{equation*}
                \begin{equation*}
                    I = 2.4mA
                \end{equation*}

        \subsection*{(d)}

        Using equation (5) in the assignment and plugging in the given values we obtain:


        \begin{align*}
            U &= -q[\frac{C_SV_S+C_GV_G+C_DV_D}{C_E}] + \frac{q^2(N-N_0)}{C_E} \\
            U &= -\alpha_SV_S - \alpha_GV_G - \alpha_DV_D + U_0(N - N_0)\:\:eV \\
            U &= -0 \cdot 0 - 0.5 \cdot 0\:eV - 0.5 \cdot 0.6\:eV + 0.25\:eV\cdot(0.325 - 0)\:\:eV \\
            U &= -0.21875\:eV
        \end{align*}
        %
        Then, starting with Equation \ref{eq:N_new} and plugging in $D(E-U) = \delta(E - \varepsilon)$, we obtain

        \begin{equation*}
            N = \frac{\gamma_1f_1(\varepsilon + U) + \gamma_2f_2(\varepsilon + U)}{\gamma_1 + \gamma_2}
        \end{equation*}
        Recalling that the expressions for the contact Fermi functions are (with $\mu = 0$ eV):

        \begin{align*}
            f_1(\varepsilon + U) = \frac{1}{1 + e^{(\varepsilon + U + qV_S)/k_BT}} = \frac{1}{1 + e^{(0.2 - 0.21875 + 0)/0.025}} = 0.679 \\
            f_2(\varepsilon + U) = \frac{1}{1 + e^{()\varepsilon + U + qV_D)/k_BT}} = \frac{1}{1 + e^{(0.2 - 0.21875 + 0.6)/0.025}} \simeq  0
        \end{align*}
        Then, solving for $\gamma_2$:

        \begin{align*}
            \gamma_2 = \gamma_2\frac{N - f_1(\varepsilon + U)}{f_2(\varepsilon + U) - N} = 5.45 \times 10^{-3} \: eV
        \end{align*}


        \subsection*{(e)}
        
        Assuming the density of states for each molecule can be modeled by $D(E) = \delta(E-\varepsilon)$,
        Equation (12) in the assignment is valid here. Thus the current will be maximized when $[f_1(E) - f_2(E)]|_{E=\varepsilon}$ is maximized. 
        Solving $[f_1(E) - f_2(E)]|_{E=\varepsilon}$ for each molecule, we obtain:\\
        Molecule A:
        \begin{equation*}
            [f_{1}(E) - f_{2}(E)]|_{E=\varepsilon_A} = 
            \frac{1}{1 + e^{\frac{\varepsilon_A}{k_BT_1}}} - \frac{1}{1 + e^{\frac{\varepsilon_A}{k_BT_2}}}
            = \frac{1}{1 + e^{\frac{0}{0.024}}} - \frac{1}{1 + e^{\frac{0}{0.027}}} = 0
        \end{equation*}
        Molecule B:
        \begin{equation*}
            [f_{1}(E) - f_{2}(E)]|_{E=\varepsilon_B} = 
            \frac{1}{1 + e^{\frac{\varepsilon_B}{k_BT_1}}} - \frac{1}{1 + e^{\frac{\varepsilon_B}{k_BT_2}}}
            = \frac{1}{1 + e^{\frac{-0.05}{0.024}}} - \frac{1}{1 + e^{\frac{-0.05}{0.027}}} = 0.02493
        \end{equation*}
        Molecule C:
        \begin{equation*}
            [f_{1}(E) - f_{2}(E)]|_{E=\varepsilon_C} = 
            \frac{1}{1 + e^{\frac{\varepsilon_C}{k_BT_1}}} - \frac{1}{1 + e^{\frac{\varepsilon_C}{k_BT_2}}}
            = \frac{1}{1 + e^{\frac{-0.1}{0.024}}} - \frac{1}{1 + e^{\frac{-0.1}{0.027}}} = 0.00877
        \end{equation*}
        Molecule B should be chosen.

        \subsection*{(f)}
        
            Referring again to equation 2 in the assignment, we know that $D(E)$ is valid for all $E$, however
            since $\gamma_1$ is dependant on the energy level, we can use it to reduce our limits. For material
            A, this means we only care about $E$ above $0eV$, and for B, $E$ below $0eV$.
            %
            Since we know the voltages on the contacts, we can calculate the effective fermi level at each.

            \begin{equation*}
                \mu_1 = \mu_0 - V_S = 0 - (-2V) = 2\:eV
            \end{equation*}
            \begin{equation*}
                \mu_2 = \mu_0 - V_D = 0 - 1V = -1\:eV
            \end{equation*}
            %
            With the fermi levels of each contact, we can adjust the limits of integration further for each
            material. Material A can be evaluated on $[0eV, 2eV]$, and material B can be evaluated on $[-1eV, 0eV]$.

            \begin{equation*}
                I_A = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot \int_{0}^{2} [1 - 0] \cdot 10^4 dE
            \end{equation*}
            \begin{equation*}
                I_A = 12.2\:mA
            \end{equation*}

            \begin{equation*}
                I_B = \frac{q}{\hbar} \cdot \frac{0.005^2}{0.01} \cdot \int_{0}^{2} [1 - 0] \cdot 10^4 dE
            \end{equation*}
            \begin{equation*}
                I_A = 6.1\:mA
            \end{equation*}
            %
            Material A should be chosen.

        \subsection*{(g)}
        
        Using Ohm's law we find the corresponding conductance:
        \begin{equation*}
            \sigma_{max} = \frac{I_{max}}{V} = \frac{500\:nA}{4.3\:mV} = 116\:\mu S.
        \end{equation*} 
        %
        We expect this result to be an integer multiple of $G_0 = 38.76$ $\mu S$, the quantum of conductance. We find

        \begin{equation*}
            \frac{\sigma_{max}}{G_0} = 3.
        \end{equation*}
        %
        We conclude there are 3 levels.



\appendix

    \newpage
    \section{Question 1b Code}

    \begin{lstlisting}[language=Matlab]
clear all;

%% Constants

% Physical constants
hbar = 1.052e-34;
q = 1.602e-19;
%epsilon_0 = 8.854e-12;
%epsilon_r = 4;
%mstar = 0.25 * 9.11e-31;

% Single-charge coupling energy (eV)
U_0 = 0.25;
% (eV)
kBT = 0.025;
% Contact coupling coefficients (eV)
gamma_1 = 0.005;
gamma_2 = gamma_1;
gamma_sum = gamma_1 + gamma_2;
% Capacitive gate coefficient
a_G = 0.5;
% Capacitive drain coefficient
a_D = 0.5;
a_S = 1 - a_G - a_D;

% Central energy level
mu = 0;

% Energy grid, from -1eV to 1eV
NE = 501;
E = linspace(-1, 1, NE);
dE = E(2) - E(1);
% TODO name this better
cal_E = 0.2;

% Lorentzian density of states, normalized so the integral is 1
D = (gamma_sum / (2*pi)) ./ ( (E-cal_E).^2 + (gamma_sum/2).^2 );
D = D ./ (dE*sum(D));

% Reference no. of electrons in channel
N_0 = 0;

voltages = linspace(0, 1, 101);
dV = voltages(2) - voltages(1);

% Terminal Voltages
V_G = 0;
V_S = 0;

for n = 1:length(voltages)
    % Set varying drain voltage
    V_D = voltages(n);

    % Shifted energy levels of the contacts
    mu_1 = mu - V_S;
    mu_2 = mu - V_D;

    % Laplace potential, does not change as solution is found (eV)
    % q is factored out here, we are working in eV
    U_L = - (a_G*V_G) - (a_D*V_D) - (a_S*V_S);

    % Poisson potential must change, assume 0 initially (eV)
    U_P = 0;

    % Assume large rate of change
    dU_P = 1;

    % Run until we get close enough to the answer
    while dU_P > 1e-6
        % source Fermi function
        f_1 = 1 ./ (1 + exp((E + U_L + U_P - mu_1) ./ kBT));
        % drain Fermi function
        f_2 = 1 ./ (1 + exp((E + U_L + U_P - mu_2) ./ kBT));

        % Update channel electrons against potential
        N(n) = dE * sum( ((gamma_1/gamma_sum) .* f_1 + (gamma_2/gamma_sum) .* f_2) .* D);

        % Re-update Poisson portion of potential
        tmpU_P = U_0 * ( N(n) - N_0);
        dU_P = abs(U_P - tmpU_P);

        % Unsure why U_P is updated incrementally, perhaps to avoid oscillations?
        %U_P = tmpU_P;
        U_P = U_P + 0.1 * (tmpU_P - U_P);
    end

    % Calculate current based on solved potential.
    % Note: f1 is dependent on changes in U but has been updated prior in the loop
    I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f_1-f_2).*D);

    if (abs(V_D-0.0) <= dV/2) 
        figure(3); title('VD = 0.0 V');
        subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
        subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
        subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
        subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
    elseif (abs(V_D-0.2) <= dV/2) 
        figure(4); title('VD = 0.2 V');
        subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
        subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
        subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
        subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
    elseif (abs(V_D-0.5) <= dV/2) 
        figure(5); title('VD = 0.5 V');
        subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.5 V');
        subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.5 V');
        subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.5 V');
        subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.5 V');
    elseif (abs(V_D-0.8) <= dV/2) 
        figure(6); title('VD = 0.8 V');
        subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.8 V');
        subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.8 V');
        subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.8 V');
        subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 0.8 V');
    elseif (abs(V_D-1.0) <= dV/2) 
        figure(7); title('VD = 1.0 V');
        subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 1.0 V');
        subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 1.0 V');
        subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 1.0 V');
        subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]);
        xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY  [eV]'); title('VD = 1.0 V');
    end

end


%%Plotting commands

figure(1); 
h = plot(voltages, N,'k'); 
grid on;
set(h,'linewidth',[2.0]); 
set(gca,'Fontsize',[18]);
xlabel('Drain voltage [V]'); 
ylabel('Number of electrons');

figure(2); 
h = plot(voltages, I,'k'); 
grid on;
set(h,'linewidth',[2.0]); 
set(gca,'Fontsize',[18]);
xlabel('Drain voltage [V]'); 
ylabel('Current [A]');
    \end{lstlisting}


\end{document}