302 lines
14 KiB
TeX
302 lines
14 KiB
TeX
\documentclass[../schrodinger_simulation.tex]{subfiles}
|
|
|
|
\begin{document}
|
|
\section{Methods}\label{sec:methods} %
|
|
\subsection{The Schrödinger equation}\label{ssec:schrodinger} %
|
|
% Add something that takes Planck to Schrödinger
|
|
% In classical mechanics, we have Newton laws and conservation of energy. In quantum
|
|
% mechanics, we have Schrödinger equation.
|
|
Erwin Schrödinger wanted to find a mathematical description of the wave characteristics
|
|
of matter, supporting the wave-particle idea. He postulated a wave function which varies
|
|
with position, where the function squared can be interpreted as the probability
|
|
of finding an electron at a given position. This resulted in the Schrödinger equation,
|
|
a wave eqution of the energy levels for a hydrogen atom. It also shows how a quantum
|
|
state evolves with time \cite[p. 81]{wu:2023:quantum}. The Schrödinger equation
|
|
has a general form
|
|
\begin{align}
|
|
i \hbar \frac{\partial}{\partial t} | \Psi \rangle &= \hat{H} | \Psi \rangle \ ,
|
|
\label{eq:schrodinger_general}
|
|
\end{align}
|
|
where $i$ is the imaginary unit, and $\hbar$ is Plancks constant. $\hat{H}$ is
|
|
a Hamiltonian operator, which represents the energy for the system, and $| \Psi \rangle$
|
|
is the quantum state. In two-dimensional position space, the quantum state can
|
|
be expressed using the time-dependent complex-valued wave function $\Psi (x, y, t)$.
|
|
Using Born rule, the square modulus of the wave function is proportional to the
|
|
probability density of detecting a particle at position $(x, y)$ at time $t$. The
|
|
relation is given by
|
|
\begin{align}
|
|
p(x, y \ | \ t) &= |\Psi(x, y, t)|^{2} = \Psi^{*}(x, y, t) \Psi(x, y, t) \ ,
|
|
\label{eq:born_rule}
|
|
\end{align}
|
|
where $\Psi^{*}$ denotes the complex conjugated wave function.
|
|
% Add something about kinetic and potential energy, to introduce the potential V
|
|
When the potential is time-independent, and the particle is non-relativistic,
|
|
the Schrödinger equation can be expressed as
|
|
\begin{align*}
|
|
i \hbar \frac{\partial}{\partial t} \Psi (x, y, t) &= - \frac{\hbar^{2}}{2m} \bigg( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \bigg) \Psi (x, y, t) \\
|
|
& \quad + V(x, y, t) \Psi (x, y, t) \numberthis \ .
|
|
\label{eq:schrodinger_special}
|
|
\end{align*}
|
|
The partial derivatives are expressions of the kinetic energy, and the potental $V$
|
|
encodes the external environment. In this experiment we will only consider the case where
|
|
the potential is time-independent, resulting in $V = V(x, y)$
|
|
|
|
When we scale Schrödinger equation by the dimensionful variables, we are left with
|
|
the wave function $u$ and the potential $v$. The dimensionless equation is given by
|
|
\begin{align}
|
|
i \frac{\partial u}{\partial t} &= - \frac{\partial^{2} u}{\partial x^{2}} - \frac{\partial^{2} u}{\partial y^{2}} + v(x, y) u \ .
|
|
\label{eq:schrodinger_dimensionless}
|
|
\end{align} %
|
|
As a result of working in position space, the Born rule is given by
|
|
\begin{align}
|
|
p(x, y \ | \ t) &= |u(x, y, t)|^{2} = u^{*}(x, y, t) u(x, y, t) \ ,
|
|
\label{eq:born_rule_scaled}
|
|
\end{align}
|
|
where we assume a normalized wave function $u(x, y, t)$. We will initialize the wave
|
|
function, using a Gaussian wavepacket, given by
|
|
\begin{align*}
|
|
u(x, y, t=0) &= e^{- \frac{(x-x_{c})^{2}}{2 \sigma_{x}^{2}} - \frac{(y-y_{c})^{2}}{2 \sigma_{y}^{2}} + ip_{x}x + ip_{y}y} \ .
|
|
\end{align*}
|
|
$x_{c}$ and $y_{c}$ are the coordinates of the center of the wavepacket, $\sigma_{x}$
|
|
and $\sigma_{y}$ are the width of the wavepacket. The wave packet momenta are
|
|
given by $p_{x}$ and $p_{y}$.
|
|
|
|
|
|
\subsection{The Crank-Nicolson scheme}\label{ssec:crank_nicolson} %
|
|
When we evaluate a particles position, we have to consider partial differential
|
|
equations (PDE). To solve these numerically, we have to discretize Equation \eqref{eq:schrodinger_dimensionless}.
|
|
We use the $\theta$-rule \footnote{Using the $\theta$-rule, we can derive Forward Euler using $\theta = 1$, and Backward Euler using $\theta = 0$},
|
|
to combine the forward (explicit) and backward (implicit) finite difference methods.
|
|
The result is a linear combination of the explicit and implicit scheme, given by
|
|
\begin{align}
|
|
\frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \theta F_{\ivec, \jvec}^{n+1} + (1 - \theta) F_{\ivec, \jvec}^{n} \ ,
|
|
\label{eq:theta_rule}
|
|
\end{align} %
|
|
where $\theta \in [0, 1]$.
|
|
To simplify the notation, and avoid any confusion of the indices with the imaginary unit $i$,
|
|
we have used the notation $\ivec, \jvec$ in subscript to indicate the commonly named indices $i, j$
|
|
in x- and y-direction. In addition, the superscript $n, n+1$ indicate position in time.
|
|
We derive the Crank-Nicolson scheme (CN) by using $\theta = 1/2$, given by
|
|
\begin{align}
|
|
\frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \ .
|
|
\label{eq:crank_nicolson_scheme}
|
|
\end{align} %
|
|
Using CN, we derive the discretized Schrödinger equation, given by
|
|
\begin{align*}
|
|
& u_{\ivec, \jvec}^{n+1} - r \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
|
|
& - r \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
|
|
&= u_{\ivec, \jvec}^{n} + r \big[ u_{\ivec+1, \jvec}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec-1, \jvec}^{n} \big] \\
|
|
& \quad + r \big[ u_{\ivec, \jvec+1}^{n} - 2u_{\ivec, \jvec}^{n} + u_{\ivec, \jvec-1}^{n} \big] - \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n} \numberthis \ ,
|
|
\label{eq:schrodinger_discretized}
|
|
\end{align*} %
|
|
where $r$ is defined as
|
|
\begin{align*}
|
|
r \equiv \frac{i \Delta t}{2 \Delta h^{2}} \ .
|
|
\end{align*} %
|
|
The full derivation of both Equation \eqref{eq:crank_nicolson_scheme} and Equation \eqref{eq:schrodinger_discretized}
|
|
can be found in Appendix \ref{ap:crank_nicolson}.
|
|
|
|
|
|
\subsection{The double-slit experiment}\label{ssec:double_slit} %
|
|
Thomas Young first performed the double-slit experiment in 1801 to demonstrate the
|
|
principle of interference of light \cite{britannica:2023:young}, while postulating
|
|
light as waves rather than particles. The double-slit experiment results in a diffraction
|
|
pattern on a detector screen, where constructive interference of light result in
|
|
bright spots, and destructive interference result in dark spots. An illustration
|
|
of Thomas Young's setup can be found in Figure \ref{fig:youngs_double_slit}.
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=\linewidth]{images/youngs_double_slit.pdf}
|
|
\caption{The setup of Thomas Young's double slit experiment, where $S_{0}$ denotes
|
|
the light source, $S_{1}$ and $S_{2}$ denotes the slits in the barrier \cite[p. 4]{mit:2004:physics}.}
|
|
\label{fig:youngs_double_slit}
|
|
\end{figure}
|
|
|
|
After the wave passes through the barrier, the pattern observed is determined by
|
|
the path difference given by
|
|
\begin{align}
|
|
\delta = d \sin (\theta) = m \lambda \ ,
|
|
\label{eq:interference}
|
|
\end{align}
|
|
where $\lambda$ is the wavelength and $m$ is called the order number. $d$ is the
|
|
distance between the center of the two slits, while assuming that the distance between
|
|
the wall and the detector screen $L >> \delta$ \cite[p. 6]{mit:2004:physics}. In
|
|
this case, we observe constructive interference when
|
|
\begin{align*}
|
|
\delta = m \lambda && m = 0, \pm 1, \pm 2 \dots \ ,
|
|
\end{align*}
|
|
and destructive interference when
|
|
\begin{align*}
|
|
\delta = (m + \frac{1}{2}) \lambda && m = 0, \pm 1, \pm 2 \dots \ .
|
|
\end{align*}
|
|
% Something about Heisenberg uncertainty principle
|
|
|
|
|
|
\subsection{Implementation}\label{ssec:implementation} %
|
|
In this experiment, we set up the grid with an equal step size in x- and y-direction $h$,
|
|
and step size in t-direction $\Delta t$, such that
|
|
\begin{align*}
|
|
x \in [0, 1] && x \rightarrow x_{\ivec} = \ivec h && \ivec = 0, 1, \dots, M-1 \\
|
|
y \in [0, 1] && y \rightarrow y_{\jvec} = \jvec h && \jvec = 0, 1, \dots, M-1 \\
|
|
t \in [0, T] && t \rightarrow t_{n} = n \Delta t && n = 0, 1, \dots, N_{t}-1 \ .
|
|
\end{align*}
|
|
In addition, we simplified the indices such that
|
|
\begin{align*}
|
|
u(x, y, t) \rightarrow u(\ivec h, \jvec h, n \Delta t) \equiv u_{\ivec, \jvec}^{n} \\
|
|
v(x, y) \rightarrow u(\ivec h, \jvec h) \equiv v_{\ivec, \jvec} \ ,
|
|
\end{align*}
|
|
which results in a matrix $U^{n}$ that contains elements $u_{\ivec, \jvec}^{n}$, and
|
|
a matrix $V$ that contains elements $v_{\ivec, \jvec}$. We used Dirichlet boundary
|
|
conditions, given by
|
|
\begin{align*}
|
|
u(x=0, y, t) &= 0 & u(x=1, y, t) &= 0 \\
|
|
u(x, y=0, t) &= 0 & u(x, y=1, t) &= 0 \ ,
|
|
\end{align*}
|
|
which allowed us to express Equation \eqref{eq:schrodinger_discretized} as a matrix
|
|
equation
|
|
\begin{align}
|
|
A u^{n+1} = B u^{n} \ .
|
|
\end{align}
|
|
Here, both $u^{n+1}$ and $u^{n}$ are column vectors containing the internal points
|
|
of the $xy$ grid at time step $n+1$ and $n$, respectively. Since we have $M$ points
|
|
in $x$- and $y$-direction, we have $M-2$ internal points. Both $u$ vectors have
|
|
length $(M-2)^{2}$, and the matrices $A$ and $B$ have size $(M-2)^{2} \times (M-2)^{2}$.
|
|
The matrices are sparse and can be decomposed as submatrices of size $(M-2) \times (M-2)$,
|
|
with the following pattern
|
|
\begin{align*}
|
|
A, B =
|
|
\begin{bmatrix}
|
|
\begin{matrix}
|
|
\bullet & \bullet & \phantom{\bullet} \\
|
|
\bullet & \bullet & \bullet \\
|
|
\phantom{\bullet} & \bullet & \bullet
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\bullet & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \bullet & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \bullet
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet}
|
|
\end{matrix} \\
|
|
\hline
|
|
\begin{matrix}
|
|
\bullet & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \bullet & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \bullet
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\bullet & \bullet & \phantom{\bullet} \\
|
|
\bullet & \bullet & \bullet \\
|
|
\phantom{\bullet} & \bullet & \bullet
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\bullet & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \bullet & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \bullet
|
|
\end{matrix} \\
|
|
\hline
|
|
\begin{matrix}
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \phantom{\bullet}
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\bullet & \phantom{\bullet} & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \bullet & \phantom{\bullet} \\
|
|
\phantom{\bullet} & \phantom{\bullet} & \bullet
|
|
\end{matrix}
|
|
& \rvline &
|
|
\begin{matrix}
|
|
\bullet & \bullet & \phantom{\bullet} \\
|
|
\bullet & \bullet & \bullet \\
|
|
\phantom{\bullet} & \bullet & \bullet
|
|
\end{matrix}
|
|
\end{bmatrix} \ .
|
|
\end{align*}
|
|
To fill the matrices $A$ and $B$, we used
|
|
\begin{align*}
|
|
a_{k} &= 1 + 4r + \frac{i \Delta t}{2} v_{\ivec, \jvec} \\
|
|
b_{k} &= 1 - 4r - \frac{i \Delta t}{2} v_{\ivec, \jvec} \ .
|
|
\end{align*}
|
|
An example of filled matrices can be found in Appendix \ref{ap:matrix_structure}.
|
|
|
|
For the general setup of the barrier, we used the values in Table \ref{tab:barrier_setup},
|
|
and for the simulations, we used the parameter settings in Table \ref{tab:sim_settings}.
|
|
% Insert Heisenberg uncertainty here? Or refer to it?
|
|
\begin{table}[H]
|
|
\centering
|
|
\begin{tabular}{l r} % @{\extracolsep{\fill}}
|
|
\hline
|
|
Parameter & Value \\
|
|
\hline
|
|
Wall thickness & $0.02$ \\
|
|
Wall position & $0.5$ \\
|
|
Separator length & $0.05$ \\
|
|
Slit aperture & $0.05$ \\
|
|
\hline
|
|
\end{tabular}
|
|
\caption{Barrier parameters and values.}
|
|
\label{tab:barrier_setup}
|
|
\end{table}
|
|
\begin{table}[H]
|
|
\centering
|
|
\begin{tabular}{l r r} % @{\extracolsep{\fill}}
|
|
\hline
|
|
Parameter & Setting 1 & Setting 2 \\
|
|
\hline
|
|
$h$ & $0.005$ & $0.005$ \\
|
|
$\Delta t$ & $2.5 \times 10^{-5}$ & $2.5 \times 10^{-5}$ \\
|
|
$T$ & $0.008$ & $0.002$ \\
|
|
$x_{c}$ & $0.25$ & $0.25$ \\
|
|
$\sigma_{x}$ & $0.05$ & $0.05$ \\
|
|
$p_{x}$ & $200$ & $200$ \\
|
|
$y_{c}$ & $0.5$ & $0.5$ \\
|
|
$\sigma_{y}$ & $0.05$ & $0.20$ \\
|
|
$p_{y}$ & $0$ & $0$ \\
|
|
$v_{0}$ & $0$ & $1 \times10^{10}$ \\
|
|
\hline
|
|
\end{tabular}
|
|
\caption{Simulation settings used in the double slit experiment. Setting 1 is
|
|
used when the barrier is switched off and setting 2 is used when the barrier
|
|
switched on.}
|
|
\label{tab:sim_settings}
|
|
\end{table}
|
|
|
|
To check if the total probability is conserved over time, and that the implementation
|
|
was correct, we computed the deviation from $1.0$ given by
|
|
\begin{align*}
|
|
s^{n} &= |1.0 - \sum_{\ivec , \jvec} p_{\ivec , \jvec}^{n}| \\
|
|
&= |1.0 - \sum_{\ivec , \jvec} u_{\ivec , \jvec}^{n*} u_{\ivec , \jvec}^{n}| \ .
|
|
\end{align*}
|
|
|
|
\subsection{Tools}\label{ssec:tools} %
|
|
The double-slit experiment is implemented in C++. We use the Python library
|
|
\verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and
|
|
\verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures.
|
|
\end{document}
|
|
|
|
|
|
|
|
% \begin{table}[H]
|
|
% \centering
|
|
% \begin{tabular}{l r} % @{\extracolsep{\fill}}
|
|
% \hline
|
|
% Position & Value \\
|
|
% \hline
|
|
% $u(x=0, y, t)$ & $0$ \\
|
|
% $u(x=1, y, t)$ & $0$ \\
|
|
% $u(x, y=0, t)$ & $0$ \\
|
|
% $u(x, y=1, t)$ & $0$ \\
|
|
% \hline
|
|
% \end{tabular}
|
|
% \caption{Boundary conditions in the xy-plane, also known as Dirichlet boundary conditions.}
|
|
% \label{tab:boundary_conditions}
|
|
% \end{table} |