The code for this project can be found here, under the final folder.
One definition of the unitarity of a matrix $U$ is that $U^\dagger U = I$. Equation 2.2
$$ \phi^{(n+1)} = \Bigg( \frac{1-a}{1+a} \Bigg) \phi^{(n)}, $$where $a \equiv iH\Delta t/2$, can be derived by applying the Hamiltonian to the average of the current and next state vector, $(\phi_{i} + \phi_{i+1})/2$ (which we define as $\chi/2$ below), and making the approximaton $\partial \phi / \partial t \rightarrow \Delta \phi / \Delta t = (\phi_{i} + \phi_{i+1})/ \Delta t$. Since this updating scheme is a linear approximation in time, $A \equiv (1+a)^{-1}(1-a)$ should unitarity as $a \rightarrow 0$, which occurs as $\Delta t \rightarrow 0$. To show this, we start by observing that
$$ (1+a)(1-a) = 1 - a^{2} $$giving, upon left multiplication by $(1+a)^{-1}$,
$$ (1+a)^{-1} = (1-a) - (1+a)^{-1}a^2 \simeq (1-a) $$to first order. Thus
$$ A = (1+a)^{-1}(1-a) \simeq (1-a)(1-a) = 1 - 2a + a^2. $$We'd like to show that $A^{T}A \rightarrow 1$ in the small $\Delta t$ limit in order to demonstrate unitarity:
$$ A^{T}A \simeq (1 - 2a + a^2)^{T} (1 - 2a + a^2) = (1 - 2a + a^2)^{T} (1 - 2a + a^2) = (1 - 2(a^{T}+a) + O(a^2)) $$which clearly converges to the identity matrix as $\Delta t \rightarrow 0$.
The operators in the Schrodinger equation can be made dimensionless by dividing through by units of energy. Choose the rest energy of the particle, $mc^2$, giving
$$ \frac{i\hbar}{mc^2}\frac{\partial}{\partial t} \phi = -\Big(\frac{\hbar}{mc}\frac{\partial}{\partial x}\Big)^2 \phi + \frac{V}{mc^2} \phi $$which can be turned into
$$ i\frac{\partial}{\partial \tilde{t}} \phi = \frac{\partial^2}{\partial \tilde{x}^2} \phi + \tilde{V} \phi $$by substituting $t$, $x$, and $V$ for their dimensionless counterparts,
$$ \frac{mc^2}{\hbar}t \rightarrow \tilde{t} $$$$ \frac{mc}{\hbar}x \rightarrow \tilde{x} $$$$ \frac{1}{mc^2}V \rightarrow \tilde{V}. $$The dimensionality of $\phi$ can similarly be removed by multiplying it by units of square root length,
$$ \frac{1}{mc^2}V \rightarrow \tilde{V}, $$and $\epsilon$ must be handled the same as $x$ in order to remove its dimensionality,
$$ \frac{mc}{\hbar}\epsilon \rightarrow \tilde{\epsilon}. $$Note that the factors $m\epsilon^2/(\hbar\Delta{t})$ and $m\epsilon^2 V/\hbar^2$ are already dimensionless but now simplify to $\tilde{\epsilon}^2/\Delta{\tilde{t}}$ and $\tilde{\epsilon}^2\tilde{V}$. If we choose to express our dimensionless time step $\Delta{\tilde{t}}$ in units of $\tilde{\epsilon}$, the $\chi$ equation (eq. 6) simplifies to
$$ \chi^{(n)}_{j+1} + \Bigg( -2 + 4i - 2\epsilon^2 V_{j} \Bigg)\chi^{(n)}_{j} + \chi^{(n)}_{j-1} = 8i\phi^{(n)}_{j}. $$where all quantities are implicitly dimensionless. In addition to simplifying the formula, this approach makes it easy to find $\Delta{\tilde{t}}$ by simply squaring $\tilde{\epsilon}$.
For the reverse time evolution operator, we get
$$ \chi^{(n)}_{j+1} + \Bigg( -2 - 4i - 2\epsilon^2 V_{j} \Bigg)\chi^{(n)}_{j} + \chi^{(n)}_{j-1} = 8i\phi^{(n)}_{j}. $$simply by defining $\Delta t$ in units of $m\epsilon^2/\hbar$, which now gives -1 since the time change is negative.
We need $\epsilon$ and $\Delta t$ to be small enough for the discretized Laplacian to work well. We're working in a box of length $L = 2$ (dimensionless length units), so $\epsilon = L/N = 2/N$ where $N+1$ is defined to be the number of discrete positions (endpoints included). In order to simulate wave numbers of order 100 with high accuracy, we should pick $N$ to be on the order of 1,000.
The probability that a particle will be found in the range $[x_{min}, x_{max}]$ can be computed via the integral
$$ P(x \in [x_{min}, x_{max}]) = \int_{x_{min}}^{x_{max}} dx \Psi^{*}\Psi $$which can be calculated in our discrete model (where $\phi$ is the discretized wavefunction) as
epsilon * dot(phi[i:j], phi[i:j])
where the indices $i, j$ correspond to $x_{min}, x_{max}$ (respectively). This is implemented via the function prob, which takes a Wave object, xmin, and xmax as arguments.
Create a wavepacket function constructor that takes $k_{0}$ and $\sigma$ as parameters and have it return a wavefunction $\phi (x)$:
$$ \Phi (k_{0}, \sigma, x_{0}) = \phi(x) = e^{ik_{0}x} e^{-(x-x_{0})^2/(2\sigma^2)} $$Using the Kronecker delta functions from Laplacian module made for the last problem set, construct the $\chi$ evolution matrix from Eq. 6, rewriting it as and use it (along with equation 5) to update the state vector, $\phi$. Implement both the $\chi$ evolution matrix and the $\phi$ evolution function in a new module called Evolution. Rewrite Eq. 6 as
where
$$ A_{i,j} \equiv -\frac{i}{8}\Big( \delta_{i+1,j} - 2 \big( 1+\epsilon^2 V_{j} -2i \big) \delta_{i,j} + \delta_{i-1,j} \Big) = -\frac{i}{8}\delta_{i+1,j} +\Big(\frac{(1 + \epsilon^2 V_{j})}{4}i + 0.5\Big)\delta_{i,j} -\frac{i}{8}\delta_{i-1,j} $$which can be efficiently solved by LU factorization and plugged into Eq. 5:
$$ \phi^{(n+1)} = \chi^{(n)} - \phi^{(n)}. $$We can preserve the boundary conditions $\phi = 0$ by setting our initial $\phi$ to zero at those points and then setting the first and last rows and columns of A to zero, except for A[1,1] and A[end,end], which should both equal 1 in order to satisfy the Thomas algorithm stability requirement that A be diagonally dominant.
For the reverse time evolution operator, we have
$$ B_{i,j} = -\frac{i}{8}\delta_{i+1,j} +\Big(\frac{(1 + \epsilon^2 V_{j})}{4}i - 0.5\Big)\delta_{i,j} -\frac{i}{8}\delta_{i-1,j}. $$Our starting wavefunction (including the normalization factor of $(1/\pi\sigma^2)^{1/4}$) is given by
$$ \phi(x) = \Big( \frac{1}{\pi\sigma^2} \Big)^{1/4} e^{i k_{0} x} e^{-(x-x_{0})^2 / 2\sigma^2}, $$from which we can extract the spread in $k$ values, $\sigma_{k}$, via Fourier transform,
$$ \tilde{\phi}(k) = \frac{1}{\sqrt{2\pi}} \Big( \frac{1}{\pi\sigma^2} \Big)^{1/4} \int_{-\infty}^{\infty} dx e^{i k_{0} x} e^{-(x-x_{0})^2 / 2\sigma^2} e^{-ikx} $$giving
$$ \tilde{\phi}(k) = \Big( \frac{\sigma^2}{\pi} \Big)^{1/4} e^{-\sigma^2(k-k_{0})^2/2}. $$The standard deviation of $k$ is $\sigma_{k} = 1 / \sqrt{2}\sigma$, giving $\sigma = 1 / \sqrt{2}\sigma_{k}$ (note that $\sigma$ itself is not the standard deviation of $x$; that's $\sigma/\sqrt{2}$). Now, we want to make sure that $\sigma_{k} << k_{0}$ in order to slow dispersion. We want the wavelength to be long enough that it doesn't suffer from the discretization of the underlying graph; pick $\lambda = 2\pi/k_{0} > 10\epsilon$, giving $k_{0} > \pi/5\epsilon$, which is approximately 600 for $N=2,000$ (this lower bound on $k_{0}$ can obviously be raised by increasing the number of grid subdivisions).
We also want a relatively narrow wavepacket so that the probability is spatially well confined, so put an upper limit on $\sigma$ of 0.05. Picking $\sigma = 0.05$, $k_{0} = 200$, then, gives $\sigma_{k}/k_{0} \approx 0.07$, which is small enough that spreading will be fairly negligible (on the order of one unit spread for every 10 units of forward propagation).
I tried the simulation with a few combinations of wavenumber $(k)$ and wavepacket width $(\sigma)$. The simulation function runs for 16,000 steps and samples a frame every 50 steps for the animations shown below. Other parameters are as noted above. The simulation is implemented in a way that makes it quite straightforward to run simulations with different combinations of parameters.
The wave packet with $k = 200$, $\sigma = 0.05$ showed little dispersion as it evolved, so I picked these parameters for part 3. Run plot2.m to see better versions of the graphs in MATLAB. You can restart the GIFs by refreshing the page.
I also tested the stability of the algorithm by running for 2,000 steps before using the time evolution operator with $\Delta{t} \rightarrow -\Delta{t}$. The effect is the same as taking the complex conjugate of $\phi$ (thus reversing its momentum) and simulating for another 2,000 steps to return it to its original state.
The average absolute error between initial and final states was -3.7646e-17, which is plenty stable for our purposes. Taking mean(abs(y2(2:end-1)).^2) - mean(abs(y1(2:end-1)).^2) (where y2 and y1 are the final and initial states, respectively) returned precisely zero up to machine precision, confirming that the matrix preserves area and is unitary to a high degree of accuracy (since $U^\dagger U$ leaves the wave vector unchanged).
The theoretical result for scattering off of a step potential boundary gives reflection and transmission coefficients
$$ T = \frac{4k_{1}k_{2}}{(k_{1}+k_{2})^2} $$$$ R = \frac{(k_{1}-k_{2})^2}{(k_{1}+k_{2})^2} $$where
$$ k_{1} = \sqrt{2E} $$$$ k_{2} = \sqrt{2(E-V_{0})} = \sqrt{k_{1}^2-2V_{0}} $$with $T=0, R=1$ in the case where $E=k^2/2<V_{0}$. Scattering a wave packet with momentum $k_{0}$ centered about $k_{1}$ should produce a similar result without the sharp cutoff at $E=k^2/2<V_{0}$, since half of the wave's momentum spectrum lies above $k_{0}$.
Setting $k_{0}=200$ as before gives $E = 20,000$. I started the wave packet at $x=-0.5$ so that the right hand side of the field would be approximately zero initially. I set up my simulation to stop running as soon as the magnitude of the wave near the right edge became greater than $10^{-3}$, at which point it would record the average value of the square of the wavefunction on both left and right hand sides of the step, and then normalize both of those averages to get the transmission and reflection probabilities. The graph below shows $T$ and $R$ vs. step potential $V_{0}$, with the theoretically predicted results also shown for comparison. There is a very steep drop-off in transmission probability between $V = 10,000$ and $V = 30,000$ and generally tight tracking of the predicted results, but with the previously noted smoothing effect due to the continuum of values in the wave packet's momentum spectrum.
Even when the step energy equals that of the central $k_{0}$ at $V = 20,000$, the higher energy components of the wave packet are still able to pass through to the right hand side. By the time $V = 30,000$, though, the contributions from higher energies is negligible, and the transmission probability drops nearly to zero.
In the limit as $h \rightarrow \infty$, the pocket potential becomes an infinite square well potential, whose energy eigenstates are sinusoids centered at zero and with wavenumbers of the form $n\pi/w$. The energy eigenvalues in that limit are $E_{n} = \hbar^2 k^2 / 2m - V_{0}$, which translates to $E_{n} = n^2\pi^2/2w^2$ using our dimensionless units.
I tried the test with the ground state and the first excited state,
$$ \phi_{ground} = \sqrt{\frac{2}{w}} cos\Big(\frac{\pi x}{w}\Big) $$$$ \phi_{1} = \sqrt{\frac{2}{w}} sin\Big(\frac{2\pi x}{w}\Big). $$I varied the thickness of the pocket walls, choosing
$$ s \in \{0.001, 0.01, 0.1\} $$and setting
$$ w = 0.2 $$$$ b = 1 $$for the width of the pocket and energy inside the pocket, giving energies
$$ E_{ground} = \frac{\pi^2}{2w^2} + b = 124.37 $$$$ E_{1} = \frac{2\pi^2}{w^2} + b = 494.48, $$which, reassuringly, is the same value calculated by $\langle \phi | H | \phi \rangle$, where $H$ is the Hamiltonian matrix calculated for the discrete simulation in the limit where the pocket walls become infinite. I chose the height of the potential walls to be roughly twice the first excited energy:
$$ h = 1,000 $$For both the excited state and the ground state, it took very little time for the "leaking" part of the wavefunction to reach the outer walls; unsurprisingly, the . The amount of leakage depended appreciably on the thickness of the wall $s$; for $s = 0.1$, half the width of the pocket itself, the leakage is miniscule:
For a thinner wall, with $s = 0.01$, the transmission probability is clearly visible for both the ground and excited starting states. There is also a visible distortion in $\phi$ at the borders of the pocket:
As expected, an extremely thin wall with $s = 0.001$ (equal to $\epsilon$, and therefore the smallest $s$ we can try simulating) produces no visible distortion in $\phi$:
As the graphs show, the leakage probability deviates significantly from being perfectly exponential. We can nonetheless clearly see that the first excited state leaks approximately four times as quickly as the ground state.
We can reuse our approach from part 3 to scatter off of the potential from part 4. We send in a wave packet with $k_0 = 2\pi/w$, the wavenumber of the first excited resonant state, followed by more wavepackets at nearby frequencies. We need $k_0$ to be larger than the resonant wave number from part 4, though, since the wave packet will disperse too quickly if it is small. So, divide the pocket width by 10:
$$ w = 0.02, $$giving a $k_{0}$ that's ten times larger. The energy is now 100 times larger
$$ E_{ground} = 12,336, $$so multiply the barrier height by 100 as well. For this run, set the wall thickness to
$$ s = 0.005 $$so that it isn't thicker than the width of the pocket. Run the simulation for multiple $k_{0}$ values near the resonance frequency and find the transmission probabilities for each.
We find that the reflection probability peaks whenever we send in a wavepacket whose central wavenumber $k_{0}$ is near an integer multiple of the ground resonant state's wavenumber.Put another way, we (locally) maximize reflection by sending in wave packets whose energies are the same as the allowed energies of the infinite square well of width $2$. I tried this for several other combinations of pocket width $w$ and barrier thickness $s$, and across several ranges of $k_{0}$, and found the same result. Also visible is a roughly exponential rise in reflection probability as $k$ increases, which can be seen pushing the reflection peak leftwards in the vicinity of $k_{0} = 3k_{ground}$.