\documentclass{article}

\usepackage{graphicx}
\usepackage{caption}
\usepackage{setspace}
\usepackage{listings}
\usepackage{color}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{float}
\usepackage[justification=centering]{caption}
\usepackage{subcaption}
\usepackage[margin=0.75in]{geometry}
\usepackage[parfill]{parskip} 
\usepackage{fancyhdr}
\usepackage{titlesec}
\usepackage{lipsum}
\usepackage{enumitem}
\usepackage{cancel}
\usepackage{siunitx}
\usepackage{adjustbox}

\titleformat{\subsection}[runin]{\normalfont \large \bfseries}
    {\thesubsection}{1em}{}

\captionsetup[figure]{width=0.75\textwidth,labelfont=normalfont,font=it,labelsep=period}

\DeclareMathOperator*{\sinc}{sinc}
\DeclareMathOperator*{\rect}{rect}
\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=2em
    }
    %\lstset{basicstyle=\small,
    %        keywordstyle=\color{mauve},
    %        identifierstyle=\color{dkgreen},
    %        stringstyle=\color{gray},
    %        numbers=left,
    %        xleftmargin=5em
    %        }


\newcommand{\angstrom}{\textup{\AA}}

\title{ECE 456 - Problem Set 2}
\date{2021-03-01}
\author{David Lenfesty \\ lenfesty@ualberta.ca
    \and Phillip Kirwin \\ pkirwin@ualberta.ca}
    
\pagestyle{fancy}

\fancyhead[L]{\textbf{ECE 456} - Problem Set 2}
\fancyhead[R]{David Lenfesty and Phillip Kirwin}
\fancyfoot[C]{Page \thepage}
\renewcommand{\headrulewidth}{1pt}
\renewcommand{\footrulewidth}{1pt}

\begin{document}
    \doublespacing
    \maketitle
    \thispagestyle{empty}
    \singlespacing
    \newpage

    \section*{Problem 1}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)] 
            \item %1a
            Code: 
            \begin{lstlisting}[language=Matlab]
clear all;
%physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;
m = 9.110e-31;

%generate lattice

N  = 100;                   %number of lattice points
n  = [1:N];                 %lattice points
a  = 1e-10;                 %lattice constant
x  = a * n;                 %x-coordinates
t0 = (hbar^2)/(2*m*a^2)/q;  %encapsulating factor
L  = a * (N+1);             %total length of consideration

%set up Hamiltonian matrix 

U = 0*x; %0 potential at all x
main_diag  = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix

H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix

[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each 
                                %column is an eigenvector
                                %"E_diag" is a diagonal matrix where the
                                %corresponding eigenvalues are on the
                                %diagonal.
                                
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues

% return eigenvectors for the 1st and 50th eigenvalues

phi_1  = eigenvectors(:,1);
phi_50 = eigenvectors(:,50);

% find the probability densities of position for 1st and 50th eigenvectors

P_1 = phi_1 .* conj(phi_1);
P_50 = phi_50 .* conj(phi_50);

% Find first N analytic eigenvalues
E_col_analytic = (1/q) * (hbar^2 * pi^2 * n.*n) / (2*m*L^2);

% Plot the probability densities for 1st and 50th eigenvectors

figure(1); clf; h = plot(x,P_1,'kx',x,P_50,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
legend('n=1','n=50');

% Plot numerical eigenvalues
figure(2); clf; h = plot(n,E_col,'kx'); grid on;
set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('EIGENVALUE NUMBER'); ylabel('ENERGY [eV]');
axis([0 100 0 40]);

% Add analytic eigenvalues to above plot

hold on;
plot(n,E_col_analytic,'k-');
legend({'Numerical','Analytical'},'Location','northwest');

            \end{lstlisting}
            %
            \begin{minipage}[t]{\linewidth}
                
            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.5\linewidth}
                    \centering
                    \includegraphics[width=\textwidth]{q1a_1and50.png}
                    \caption{}
                    \label{fig:q3_ne}
                \end{subfigure}%
                \begin{subfigure}{0.5\linewidth}
                    \centering
                    \includegraphics[width=\textwidth]{q1a_eigenvals.png}
                    \caption{}
                    \label{fig:q3_current}
                \end{subfigure}
                \caption{(a) Probability densities for \(n=1\) and \(n=50\). 
                            (b) Comparison of first 101 numerical and analytic eigenvalues.}
            \end{figure}
            \end{minipage}
            \item %1b
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)] 
                \item %1bi
                The analytical solution is:

                \begin{equation}
                    \phi(x) = A \sin{\left( \frac{n \pi}{L} x \right)}. 
                    \label{eq:analytical}
                \end{equation}

                In order to normalise this equation it must conform to the following:

                \begin{equation}
                    \int_{0}^{L} \left| \phi(x) \right| ^2 dx = 1.
                    \label{eq:normalise}
                \end{equation}

                We use the following identity:

                \begin{equation}
                    \int\sin^2{(a x)}\:dx = \frac{1}{2} x -\frac{1}{4a}\sin{(2a x)}.
                    \label{eq:integral_ident}
                \end{equation}

                Given that the sine of a real value is always real, we can disregard the norm operation, and directly relate (\ref{eq:analytical}) 
                to the above identity.
                Evaluating the integral gives us the following relationship:

                \begin{equation*}
                    \frac{1}{A^2} = \frac{1}{2} L - \cancel{\frac{L}{4n \pi} \sin{\left( \frac{2n \pi}{L} L  \right)} }
                     - \cancel{\frac{1}{2} \cdot 0} + \cancel{\frac{L}{4n \pi} \sin{(0)}}.
                \end{equation*}

                From this, we find:

                \begin{equation*}
                    \boxed{A = \sqrt{\frac{2}{L}}.}
                \end{equation*}

                \item %1bii

                Starting with the normalization condition for the numerical case:
                \begin{align}
                    a\sum_{\ell=1}^N\left|\phi_l\right|^2  &= a \\
                    a\sum_{\ell=1}^N\left|B\sin{\left(\frac{n\pi}{L}x_\ell\right)}\right|^2 &= a, \nonumber
                    \label{eq:numerical_norm}
                \end{align}
                recalling that \(x = a\ell\), and allowing \(a \to 0\), while holding \(L\) 
                constant, implies that \(N \to \infty\), since \(a=\frac{L}{N}\). 
                An integral is defined as the limit of a Riemann sum as follows:
                \begin{equation}
                    \int_c^df(x)\:dx \equiv \lim_{n \to \infty}\sum_{i=1}^n\Delta x\cdot f(x_i),
                    \label{eq:Riemann_sum}
                \end{equation}
                where \(\Delta x = \frac{d-c}{n}\) and \(x_i = c + \Delta x \cdot i\). In our case,
                \(n = N\), \(i = \ell\), \(c = 0\), \(d = L\), and \(\Delta x = a\), \(x_i = x_\ell\), 
                \(f(x) = \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\). Therefore we can write
                \begin{equation*}
                    \int_0^L \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\:dx 
                    = \lim_{N \to \infty}\sum_{\ell=1}^N a \cdot \left|B\sin{\left(\frac{n\pi}{L}x_l\right)}\right|^2
                    = a.
                \end{equation*}
                
                Using (\ref{eq:integral_ident}), we have
                \begin{equation*}
                    \int_0^L \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\:dx 
                    = \frac{1}{2}L - \cancel{\frac{L}{4n\pi}\sin{\left(\frac{2n\pi}{L}L\right)}}
                        - 0 + 0
                    = \frac{a}{B^2}.
                \end{equation*}

                This means that \(B\) must be
                \begin{equation*}
                    \boxed{B = \sqrt{\frac{2 a}{L}} = \sqrt{a} \times A.}
                \end{equation*}

            \end{enumerate}

            \item %1c
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %1ci

                % TODO check that B is expressed correctly in all these equations (alpha, not ell, wait for pdf gen)
                From the base form of $\phi _\ell = B\sin{\left( \frac {n \pi}{L} a \ell \right)} $, we can see that
                $\phi _ {\ell + 1}$ and $\phi _ {\ell - 1}$ correspond to the trigonometric identities
                $\sin{(a + B)} = \sin{(a)}\cos{(B)} + \cos{(a)}\sin{(B)}$ and
                $\sin{(a + B)} = \sin{(a)}\cos{(B)} + \cos{(a)}\sin{(B)}$, respectively, where
                $a = \frac{n \pi a \ell}{L}$ and $B = \frac{n \pi a}{L}$.

                Plugging these identities into equation (7) from the assignment and simplifying, we get to this equation:

                \begin{equation*}
                    -t_0 B \sin{\left( \frac{n \pi a \ell}{L} \right)}  + 2 t_0 \phi _{\ell} - t_0 \sin{\left( \frac{n \pi a \ell}{L} \right)}. 
                \end{equation*}

                At this point, we notice that $\phi _ \ell = B \sin{\left( \frac {n \pi}{L} a \ell \right)} $,
                so we can factor it out.

                With some minor rearranging, this leaves us with the final expression for $E$:

                \begin{equation}
                    \boxed{E = 2t_0 \left( 1 - \cos{\left( \frac{n \pi a}{L} \right)}  \right).}
                    \label{eq:numerical_E}
                \end{equation}
                
                \item %1cii
                
                \begin{minipage}[t]{\linewidth}

                    \centering
                    \adjustbox{valign=t}{
                        \includegraphics[width=0.5\textwidth]{q1cii.png}
                    }
                    \captionof{figure}{Analytic solution to numeric system, plotted.}
                \end{minipage}

                We can see here that the "predicted" numerical response matches nearly exactly the actual
                calculated numerical solution.

                \item %1ciii

                Applying the approximation $\cos{(\theta)} \approx 1 - \frac{\theta^2}{2}$ for small $\theta$
                on equation (\ref{eq:numerical_E}), we get the following expression:

                \begin{equation*}
                    E = 2 t_0 \left( \frac{n^2 \pi^2 a^2}{2 L^2} \right).
                \end{equation*}

                 We can get our final analytical expression for $E$ by fully substituting the explicit form of $t_0$:

                \begin{equation*}
                    \boxed{E = \frac{ \hbar^2 n^2 \pi^2 }{ 2 m L^2 }}
                \end{equation*}

                \item %1civ
                With the decreased lattice spacing and increased number of points we can see the numerical solution
                more closely matches the analytical solution. As well, the \(n=50\) case is 
                now a constant-amplitude wave, which corresponds to the expected analytic result, in contrast to the 
                plot in section (a), which has a low-frequency envelope around it.
                
                \begin{minipage}[t]{\linewidth}
                    
                \begin{figure}[H]
                    \centering
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1civ_fig1.png}
                        \caption{}
                    \end{subfigure}%
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1civ_fig2.png}
                        \caption{}
                    \end{subfigure}

                    \caption{(a) Probability densities for \(n=1\) and \(n=50\). (b) Comparison of first 101 numerical and analytic eigenvalues.}

                \end{figure}
                \end{minipage}
                

            \end{enumerate}
            \item %1d
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %1di

                In order to modify the computations to those for a particle in a "ring"
                we simply had to add $-t_0$ elements as the "corner" elements of the
                hamiltonian operator array:
                \begin{lstlisting}[language=Matlab]
% Modify hamiltonian for circular boundary conditions
H(1, N) = -t0;
H(N, 1) = -t0;
                \end{lstlisting}

                \begin{minipage}[t]{\linewidth}

                \begin{figure}[H]
                    \centering
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1di_fig1.png}
                        \caption{}
                    \end{subfigure}%
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1di_fig2.png}
                        \caption{}
                    \end{subfigure}
                    \caption{(a) Probability densities for \(n=4\) and \(n=5\). (b) Comparison of first 101 numerical and analytic eigenvalues.}
                \end{figure}
                \end{minipage}

                \item %1dii

                The energy levels for eigenvalues number 4 and 5 are both \boxed{0.06\:\mathrm{eV}}.
                These eigenstates are degenerate because they both have the same 
                eigenvalue/energy.

                \item %1diii

                \begin{minipage}[t]{\linewidth}
                    \centering
                    \adjustbox{valign=t}{
                        \includegraphics[width=0.5\textwidth]{q1div.jpg}
                    }
                    \captionof{figure}{Sketch of degenerate energy levels.}
                \end{minipage}

                \item %1div

                Plugging the valid levels for $n$ into equation (10) from the assignment (and dividing
                by the requisite $q$), we get energy levels of \(0\) eV, \(0.0147\) eV, and \(0.0589\) eV for the
                indices \(n=0\), \(1\), and \(2\), respectively. These match very closely, to within
                acceptable margin of the numerical results from part (ii).


            \end{enumerate}

        \end{enumerate}

    \newpage
    \section*{Problem 2}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)] 
           
            \item %2a
            Since \(t_0 = \frac{\hbar^2}{2ma^2}\), and \(U_l = -\frac{q^2}{4\pi\varepsilon_0r_l}
            + \frac{l_o(l_o + 1)\hbar^2}{2mr_l^2}\), the middle diagonal elements will have values
            \begin{equation*}
                \boxed{\hat{H}_{ll} = \frac{\hbar^2}{ma^2}-\frac{q^2}{4\pi\epsilon_0r_l} + \frac{l_o(l_o + 1)\hbar^2}{2mr_l^2},}
            \end{equation*}
            and the upper and lower diagonals elements will have values
            \begin{equation*}
                \boxed{\hat{H}_{l(l\pm1)} = -\frac{\hbar^2}{2ma^2}.}
            \end{equation*}
            \item %2b
            Homogenous boundary conditions imply that the corner entries of \(\hat{H}\) will be \boxed{0}.
            \item %2c
            Code:
            \begin{lstlisting}
clear all;
%physical constants in MKS units


hbar = 1.054e-34;
q = 1.602e-19;
m = 9.110e-31;
epsilon_0 = 8.854e-12;

%generate lattice

N  = 100;                   %number of lattice points
n  = [1:N];                 %lattice points
a  = 0.1e-10;                 %lattice constant
r  = a * n;                 %x-coordinates
t0 = (hbar^2)/(2*m*a^2)/q;  %encapsulating factor
L  = a * (N+1);             %total length of consideration

%set up Hamiltonian matrix 

U = -q^2./(4*pi*epsilon_0.*r) * (1/q); %potential at r in [eV]
main_diag  = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix

H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix

[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each column is an eigenvector
                                %"E_diag" is a diagonal matrix where the
                                %corresponding eigenvalues are on the
                                %diagonal.
                                
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues

% return eigenvectors for the 1st and 50th eigenvalues

phi_1  = eigenvectors(:,1);
phi_2 = eigenvectors(:,2);

% find the probability densities of position for 1st and 50th eigenvectors

P_1 = phi_1 .* conj(phi_1);
P_2 = phi_2 .* conj(phi_2);

% Plot the probability densities for 1st and 2nd eigenvectors

figure(1); clf; h = plot(r,P_1,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.02 0.04 0.06 0.08 0.10 0.12]);
legend('n=1');
axis([0 1e-9 0 0.12]);

figure(2); clf; h = plot(r,P_2,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04]);
legend('n=2');
axis([0 1e-9 0 0.04]);
            \end{lstlisting}
            \begin{minipage}[t]{\linewidth}

            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q2c_fig1_1s.png}
                    \caption{}
                \end{subfigure}%
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q2c_fig2_2s.png}
                    \caption{}
                \end{subfigure}
                \caption{(a) 1s probability density. (b) 2s probability density.}
            \end{figure}
            \end{minipage}
            \item %2d
            For the 1s level, \boxed{E = -13.4978\:\mathrm{eV}}.
            \item %2e
            Beginning with equation (11) from the assignment, with \(l_o=0\):
            \begin{align*}
                \left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2} - \frac{q^2}{4\pi\epsilon_0r}\right]f(r) &= Ef(r)\\
                -\frac{\hbar^2}{2m}\frac{d^2}{dr^2}\left(\frac{2r}{a_0^{3/2}}e^{-r/a_0}\right) - \frac{q^2}{4\pi\epsilon_0r}f(r) &= Ef(r)\\
                -\frac{\hbar^2}{2m}\frac{2}{a_0^{3/2}}\frac{d}{dr}\left(e^{-r/a_0} - \frac{r}{a_0}e^{-r/a_0}\right) - \frac{q^2}{4\pi\epsilon_0r}f(r) &= Ef(r)\\
                -\frac{\hbar^2}{2m}\frac{2}{a_0^{3/2}}\left(-\frac{1}{a_0}e^{-r/a_0} - \frac{1}{a_0}e^{-r/a_0} + \frac{r}{a_0^2}e^{-r/a_0}\right) - \frac{q^2}{4\pi\epsilon_0r}f(r) &= Ef(r)\\
                -\frac{\hbar^2}{2m}\left(-\frac{2}{a_0r} + \frac{1}{a_0^2}\right)f(r) - \frac{q^2}{4\pi\epsilon_0r}f(r) &= Ef(r)\\
                -\frac{\hbar^2}{2m}\left(-\frac{2}{a_0r} + \frac{1}{a_0^2}\right) - \frac{q^2}{4\pi\epsilon_0r} &= E.
            \end{align*}
            Recalling that \(a_0 = 4\pi\epsilon_0\hbar^2/mq^2\), we can eliminate \(r\):
            \begin{align*}
                \cancel{\frac{\hbar^2}{2m}}\frac{\cancel{2m}q^2}{4\pi\epsilon_0\cancel{\hbar^2}}\frac{1}{r} - \frac{\hbar^2}{2m}\frac{1}{a_0^2} - \frac{q^2}{4\pi\epsilon_0}\frac{1}{r} &= E\\
              \cancel{\frac{q^2}{4\pi\epsilon_0}\frac{1}{r}} -  \frac{\hbar^2}{2m}\frac{1}{a_0^2} - \cancel{\frac{q^2}{4\pi\epsilon_0}\frac{1}{r}} &= E.
            \end{align*}
            We can then solve for \(E\):
            \begin{equation*}
                E\:\mathrm{[eV]} = -\frac{1}{q}\cdot\frac{\hbar^2}{2ma_0^2} = -\frac{1}{q}\cdot\frac{(\SI{1.054e-34}{\joule\cdot\second})^2}{2(\SI{9.110e-31}{\kilogram})(\SI{0.0529}{\nm})^2}
                    = \boxed{- 13.6\:\mathrm{eV}.}
            \end{equation*}
            This is very similar to the result in (d).

            \item %2f
            In the figure below we can see that the numerical and analytical results agree up to scaling by \(a\). The scale difference is expected,
            as discussed in Problem 1. From (d), we also expect agreement in the curve shapes because the numerical and analytical energies for the 1s
            level are very similiar. We can see that the peak value of the analytic result is very slightly higher than that of the numerical result,
            which corresponds to the analytical result for the energy being slightly greater in magnitude (\(-13.6\) eV versus \(-13.4978\) eV).
            \begin{minipage}[t]{\linewidth}
                \centering
                \includegraphics[width=0.5\textwidth]{q2f_fig.png}
                \captionof{figure}{Numerical result (black line) and analytical solution scaled by \(a\) (orange circles).}
            \end{minipage}
            


        \end{enumerate}

    \newpage    
    \section*{Problem 3}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)]
           
            \item %3a

            \begin{equation}
                \theta (x') = \sqrt{\frac{2}{L}} sin \left( \frac{\pi (x' + L/2)}{L} \right) = \sqrt{\frac{2}{L}} cos \left( \frac{\pi x'}{L} \right)
            \end{equation}

            \item %3b

            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %3bi

                Mapping the provided Fourier identities from \( t \) and \( \omega \) onto
                \(x'\) and \(k'\), we can evaluate the Fourier transform \(A(k')\):

                \begin{align*}
                    \mathcal{F}\left[rect(\frac{x'}{L}) \right] &= \frac{L}{\sqrt{2 \pi}} sinc \left( \frac{k' L}{2 \pi} \right) \\
                    \mathcal{F}\left[f(x')cos(\frac{\pi}{L}x')\right] &= \frac{1}{2} \left[ F(k' + k_1) + F(k' - k_1) \right] \\
                    A(k') &= \frac{1}{2} \sqrt{\frac{L}{\pi}} \left\{ \sinc \left( \frac{L}{2\pi} (k' + k_1) \right) + \sinc \left( \frac{L}{2\pi} (k' - k_1) \right) \right\}.
                \end{align*}

                \item %3bii
                Beginning with the result for \(A(k')\) above, and writing
                \begin{equation*}
                    \Phi(p') \equiv \frac{1}{\sqrt{\hbar}} A \left(\frac{p'}{\hbar}\right),
                \end{equation*}
                we can obtain
                \begin{equation*}
                    \Phi(p') = \frac{1}{2} \sqrt{\frac{L}{\pi\hbar}} \left\{ \sinc \left( \frac{L}{2\pi\hbar} (p' + p_1) \right) + \sinc \left( \frac{L}{2\pi\hbar} (p' - p_1) \right) \right\}.
                \end{equation*}
                

                \item %3biii
                \(\left|\Phi(p')\right|^2\) has units of [\si{\s.kg^{-1}.m^{-1}}], which are those of inverse momentum. Thus, multiplication 
                (or integration) by a differential of momentum results in a unitless probability, as we should expect. This holds in the 1D case
                and can easily be generalized to higher dimensions.
                

                \item %3biv

                \( sinc \) is a purely real function, so we can ignore the normalizing portion of the integral.
                As well, to simplify the interim equations we will assign the constants \( A = \frac{1}{2} \sqrt{\frac{L}{\pi \hbar}} \)
                and \( B = \frac{L}{2 \pi \hbar} \).

                \begin{equation*}
                    \int _{-\infty} ^{\infty} \Phi (p')^2 dp = \int _{-\infty} ^{\infty} A^2 \left[ sinc\left( B(p'+p_1) \right)^2 + 2 sinc \left( B(p'+p_1) \right) sinc \left( B(p'-p_1) \right) + sinc \left( B(p'-p_1) \right) \right] dp
                \end{equation*}

                Given property (26) of the \( sinc \) function, we can simplify the left and right terms.
                Using a change of variable \( p'' = p_1 - p' \), and properties (27) and (28), we can further evaluate the central term.

                \begin{align}
                    \int _{-\infty} ^{\infty} \Phi (p')^2 dp = \int _{-\infty} ^{\infty} A^2 \left[2 sinc \left( B(2p_1 - p'') \right) sinc \left( B(-p'') \right) \right] dp + A^2 \frac{2}{B} \\
                    &= \int _{-\infty} ^{\infty} A^2 \left[ 2 sinc \left(B(2p_1 - p'') \right) sinc(B p'') \right] dp + A^2 \frac{2}{B} \\
                    &= \frac{A^2}{B} \left[ sinc(2p_1) + 2 \right] \\
                    &= \frac{1}{2} \left[ sinc(2p_1) + 2 \right]
                \end{align}

                \item %3bv
                
                
            \begin{minipage}[t]{\linewidth}

            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q3bv_fig1.png}
                    \caption{}
                \end{subfigure}%
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q3bv_fig2.png}
                    \caption{}
                \end{subfigure}
                \caption{(a) momentum wavefunction versus normalized momentum. (b) Probability density versus normalized momentum.}
            \end{figure}
            \end{minipage}
                \item %3bvi
                The points of classical momentum are given by \( p_1 = \pm \sqrt{2mE}\). On the normalized plots, these occur at \(\pm 1\) on the \(p/p_1\) axis.
                Given that \(L = \SI{101}{\angstrom}\), we can find the velocity of the electron by taking \(v = \frac{p_1}{m_e}\). We find that
                \(v = \pm \SI{3.6e4}{m/s}\).
                
            \begin{minipage}[t]{\linewidth}

            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q3bvi_fig1.png}
                    \caption{}
                \end{subfigure}%
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q3bvi_fig2.png}
                    \caption{}
                \end{subfigure}
                \caption{(a),(b) Previous plots but with the classical momentum marked in red.}
            \end{figure}
            \end{minipage}
                \item %3bvii
                From the plot of the probability density, we can clearly see that the particle can take a continuum of momentum values.
                Thus the statement is false.
            \end{enumerate}

            \item %3c
            Because the probability density is even about \(p' = 0\), we can surmise that \(\left\langle p'\right\rangle = 0\).
            \newline
            To verify this, we find \(\left\langle p'\right\rangle\) 
            from \(\phi(x') = \sqrt{\frac{2}{L}}\sin{\left(\frac{\pi}{L}x'\right)}\times\rect{\left(\frac{x'}{L}\right)}\) according to
            \begin{align*}
                \left\langle p'\right\rangle\ &= \int_{-\infty}^\infty \phi^*(x')\:\hat{p}\:\phi(x')\:dx'\\
                \left\langle p'\right\rangle\ &= -i\hbar\int_{-L/2}^{L/2} \sqrt{\frac{2}{L}}\sin{\left(\frac{\pi}{L}x'\right)}\:\frac{d}{dx'}\left[\sqrt{\frac{2}{L}}\sin{\left(\frac{\pi}{L}x'\right)}\right]\:dx'\\
                \left\langle p'\right\rangle\ &= \frac{-2i\pi\hbar}{L^2}\int_{-L/2}^{L/2} \sin{\left(\frac{\pi}{L}x'\right)}\:\cos{\left(\frac{\pi}{L}x'\right)}\:dx'.
            \end{align*}
            Using Equation (31) in the assignment we can write

            \begin{align*}
                \left\langle p'\right\rangle\ &= \left.\frac{-i\hbar}{L}\sin^2{\left(\frac{\pi}{L}x'\right)}\right|_{-L/2}^{L/2}\\
                \left\langle p'\right\rangle\ &= \frac{-i\hbar}{L}\left(1 - 1\right) = 0,
            \end{align*}
            which verifies our above inference.
            \item %3d
            
            The momentum associated with the wavefunction \( \theta (x') = e^{ik'x'} \) is sharp,
            and the corresponding value is \( p' = \hbar k' \).
            
            \begin{equation*}
                \hat{p} = -i \hbar \frac{d}{dx'} e ^{i k' x'} = -i^2 \hbar k' e ^{i k' x'} = \hbar k' e^{i k' x'}
            \end{equation*}

        \end{enumerate}
    \newpage 
    \section*{Problem 4}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)]
           
            \item %4a
            \begin{minipage}[t]{\linewidth}

            \begin{figure}[H]
                \centering
                \begin{subfigure}[b]{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q4a_e1.png}
                    \caption{}
                \end{subfigure}%
                \begin{subfigure}[b]{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q4a_e2.png}
                    \caption{}
                \end{subfigure}
                \caption{Probability Density Plots for first two energy levels.}
            \end{figure}
            \end{minipage}

            An \( a \) value of \( 0.53 \angstrom \) was chosen in order to provide an adequately
            shaped graph without sacrificing too much computation time and ensure that the first two
            numerical energies correspond to the given experimental results. The experimental results are
            \(\SI{0.14395}{\electronvolt}\) and \(\SI{0.43185}{\electronvolt} \) for the first and second energy levels 
            respectively, and the numerical results with our chosen \(a\) are \(\SI{0.14386}{\electronvolt}\) and \(\SI{0.43140}{\electronvolt} \),
            which are in agreement.




            \item %4b
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %4bi

                The energies were used were \( 0.14395eV \) and \( 0.43185eV \) for the first and second energy levels,
                respectively.

                \begin{minipage}[H]{\linewidth}
                    \centering
                    \includegraphics[width=\textwidth]{q4bi_I-V.png}
                    \captionof{figure}{Current vs. Voltage.}
                \end{minipage}

                \item %4bii

                Between \( 0V \) and \( 0.25V \), only the first energy level is carrying any current.
                This current drops to 0 above \( 0.25V\) because the coupling between the contacts and that
                energy level drops to 0, meaning no electrons can transfer.

                Between \(0.4V\) and \(0.65V\), only the second energy level is carrying current.
                This energy level stops conducting current because it's shifted energy drops below
                the threshold where the contacts have any coupling with it.

                \item %4biii
                
                Negative Drain Resistance (NDR) is present in this design from a \(V_D\) of
                approximately \(0.27V\) to \(0.45V\), as well as from \(0.65V\) to \(0.8V\)

            \end{enumerate}

        \end{enumerate}
        
    
\end{document}