Scattering is the process where a wave or particle collides with some object. An understanding of the scattering process allows one to look at objects that may not be easily observed: for example, in the Rutherford scattering experiment, the behaviour of particles scattering off gold foil led him to realise that most of the positive charges were concentrated in the centre of the atom, rather than being spread throughout the whole atom, as was the consensus in the 1900s.[1]
As such, the framework of scattering theory is particularly important in quantum mechanics, where direct observation is difficult – at times impossible. Rather, by throwing known particles at an unknown object, and studying how it interacts with them, one reconstructs the properties of the object itself.
In the first two chapters of Quantum Mechanics III, we were introduced to a formulation of the scattering theory under the Schrödinger framework. This formulation is non-relativistic: it takes the energy of a free particle to be \(E = \frac{p^2}{2m}\), and hence is only appropriate for slow-moving particles.
Here, we looked at the Dirac equation, a relativistic framework that marries quantum mechanics with special relativity – now \(E^2 = {\left(mc^2\right)}^2 + {\left(pc\right)}^2\), as is more accurate for particles moving closer to the speed of light. The scattering theory is then reformulated with the Dirac equation in place of the Schrödinger equation to acquire a scattering theory that is also relativistic.
In this report, we are interested in extending the non-relativistic scattering theory into the relativistic realm using the Dirac Hamiltonian, building upon the material established in the course notes.[2]
The free, single-electron Dirac Hamiltonian is given by
\begin{equation} \hat{H}_0 = \hat{\alpha}_0mc^2 + \sum_{j=1}^3\hat{\alpha}_j\hat{p}_jc = \hat{\alpha}_0mc^2 + \hat{\boldsymbol{\alpha}}\cdot\hat{\mathbf{p}}c\label{eq:H} \end{equation}
where \(\alpha_\mu\) are the standard Dirac matrices, \(\hat{\boldsymbol{\alpha}} = {\left[\hat{\alpha}_1\;\hat{\alpha}_2\;\hat{\alpha}_3\right]}\), and \(\hat{\mathbf{p}} = {\left[\hat{p}_x\;\hat{p}_y\;\hat{p}_z\right]}\).
The square of the Dirac Hamiltonian in the position representation, which will be useful later, is
\begin{equation} {\hat{H}_0}^2 = \left(m^2c^4 + {\left\lvert\hat{\mathbf{p}}\right\rvert}^2c^2 \right)\hat{I} = \left(m^2c^4 - \hbar^2c^2\nabla^2\right)\hat{I} \label{eq:H2} \end{equation}
From here, operators are to be understood as given in the position representation. Considering the dimensionality of the Dirac matrices, the wavefunction \(\psi(\mathbf{r}, t)\) is a 4-column vector.
The plane wave eigenstates are§
\begin{equation} \psi(\mathbf{r}, t) = \phi e^{i\left(\mathbf{k}\cdot\mathbf{r}-\omega t\right)} = \begin{bmatrix}\phi_+ \\ \phi_-\end{bmatrix}e^{i\left(\mathbf{k}\cdot\mathbf{r}-\omega t\right)} \end{equation}
with \(\omega = \frac{E}{\hbar}\) and
\begin{align} \phi_+ &= \frac{\hbar c}{E-mc^2}\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_-\quad(\text{if }E < 0) & \phi_- &= \frac{\hbar c}{E+mc^2}\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_+\quad(\text{if }E > 0) \label{eq:PMrelation} \end{align}
The scattering formulation from Chapter 0 of PH4401 still stands for the Dirac Hamiltonian without any modification needed, as no “illegal” moves have been made with regards to the wavefunction now being a four-component vector.
Therefore, given a Hamiltonian \(\hat{H} = \hat{H}_0 + \hat{V}(\mathbf{r})\), where \(\hat{V}(\mathbf{r})\) is a 4-by-4 matrix describing the scattering potential, \(\psi_i(\mathbf{r})\) describes an incident particle of a particular energy \(E\), and \(\psi(\mathbf{r})\) describes the total wave function,
Since the form of the Green's function (denoted \(G_{0,s}\), with the non-relativistic energy \(E_s = \frac{p^2}{2m}\), where \(p\) is the momentum of the incident particle) in the Schrödinger case is
\begin{equation} \left[E_s - \frac{\hat{\mathbf{p}}^2}{2m} \right]G_{0,s}(\mathbf{r}, \mathbf{r'}) = \frac{1}{2m}\left[p^2 + \hbar^2\nabla^2 \right]G_{0,s}(\mathbf{r}, \mathbf{r'}) = \delta^d\left(\mathbf{r} - \mathbf{r'}\right)\label{eq:greenS} \end{equation}
we can try swapping out the free-particle Hamiltonian with the appropriate Dirac operators:
\begin{equation} \left[E\hat{I} - \hat{H}_0\right]\hat{G}_0(\mathbf{r}, \mathbf{r'}) = \left[E\hat{I} - \hat{\alpha}_0mc^2 - \hat{\boldsymbol{\alpha}}\cdot\hat{\mathbf{p}} \right]\hat{G}_0(\mathbf{r}, \mathbf{r'}) = \delta^d\left(\mathbf{r} - \mathbf{r'}\right)\hat{I}\label{eq:greenD} \end{equation}
Note that the \(E\) here is the relativistic energy of the incident particle. We postulate that the Green's function plays the same role with the Dirac equation. That is, given \(\hat{G}_0\) (possibly a 4-by-4 matrix) that satisfies the above equation,
\begin{equation} \psi(\mathbf{r}) = \psi_i(\mathbf{r}) + \underbrace{\int d\mathbf{r'}\,\hat{G}_0(\mathbf{r}, \mathbf{r'})\hat{V}(\mathbf{r'})\psi(\mathbf{r'})}_{\psi_s(\mathbf{r})}\label{eq:preBorn} \end{equation}
Instead of deriving the above, we will instead show that this is correct by applying the Hamiltonian to both sides of the equation:
\begin{align} \hat{H}\psi(\mathbf{r}) &= \left[\hat{H}_0 + \hat{V}(\mathbf{r})\right]\left(\psi_i(\mathbf{r}) + \int d\mathbf{r'}\,\hat{G}_0(\mathbf{r}, \mathbf{r'})\hat{V}(\mathbf{r'})\psi(\mathbf{r'})\right)\nonumber\\ E\psi(\mathbf{r}) &= \left[E + \hat{V}(\mathbf{r})\right]\psi_i(\mathbf{r}) + \left[\hat{H}_0 - E\right]\int d\mathbf{r'}\,\hat{G}_0(\mathbf{r}, \mathbf{r'})\hat{V}(\mathbf{r'})\psi(\mathbf{r'}) + \left(\hat{V}(\mathbf{r}) + E\right)\psi_s(\mathbf{r})\nonumber \\ &= E\left(\psi_i(\mathbf{r}) + \psi_s(\mathbf{r})\right) + \hat{V}(\mathbf{r})\left(\psi_i(\mathbf{r}) + \psi_s(\mathbf{r})\right) - \int d\mathbf{r'}\,\smash{\underbrace{\left(E - \hat{H}_0\right)\hat{G}_0(\mathbf{r}, \mathbf{r'})}_{\delta^d\left(\mathbf{r} - \mathbf{r'}\right)\hat{I}}}\hat{V}(\mathbf{r'})\psi(\mathbf{r'})\nonumber\\ &= E\psi(\mathbf{r}) + \hat{V}(\mathbf{r})\psi(\mathbf{r}) - \hat{V}(\mathbf{r})\psi(\mathbf{r}) \nonumber\\[5pt] &= E\psi(\mathbf{r}) \end{align}
As the one-, two-, and three-dimensional Green's functions in the Schrödinger case have already been found in the course notes, it makes sense to equate Equation \ref{eq:greenS} & Equation \ref{eq:greenD} to relate them with the Dirac Green's function so as to use our prior result.
\begin{equation} \left[E\hat{I} - \hat{H}_0\right]\hat{G}_0(\mathbf{r}, \mathbf{r'}) = \frac{1}{2m}\left[p^2 + \hbar^2\nabla^2 \right]G_{0,s}(\mathbf{r}, \mathbf{r'})\hat{I}\label{eq:greenTot} \end{equation}
Since the relativistic energy of the incident particle is related to its momentum by
\begin{equation} E^2 = m^2c^4 + p^2c^2 \implies p^2 = \frac{E^2-m^2c^4}{c^2} \end{equation}
we place this back into Equation \ref{eq:greenTot}.
\begin{align} \left[E\hat{I} - \hat{H}_0\right]\hat{G}_0(\mathbf{r}, \mathbf{r'}) &= \frac{1}{2mc^2}\left[E^2 - m^2c^4 + \hbar^2c^2\nabla^2 \right]G_{0,s}(\mathbf{r}, \mathbf{r'})\hat{I}\nonumber\\ &= \frac{1}{2mc^2}\left[E^2\hat{I} - \left(m^2c^4 - \hbar^2c^2\nabla^2\right)\hat{I} \right]G_{0,s}(\mathbf{r}, \mathbf{r'})\hat{I} \nonumber\\ &= \frac{1}{2mc^2}\left[{\left(E\hat{I}\right)}^2 - {\hat{H}_0}^2 \right]G_{0,s}(\mathbf{r}, \mathbf{r'})\hat{I}\nonumber\\ &= \left[E\hat{I} - {\hat{H}_0} \right]\left(\left[E\hat{I} + {\hat{H}_0} \right]\frac{G_{0,s}(\mathbf{r}, \mathbf{r'})}{2mc^2}\hat{I}\right)\nonumber\\ \therefore \hat{G}_0(\mathbf{r}, \mathbf{r'}) &= \left[E\hat{I} + {\hat{H}_0} \right]\frac{G_{0,s}(\mathbf{r}, \mathbf{r'})}{2mc^2}\hat{I}\label{eq:greenRel} \end{align}
Here Equation \ref{eq:H2} is used in the second line, and the fact that \(E\hat{I}\) commutes with \(\hat{H}_0\) is used in the third.
By repeatedly placing Equation \ref{eq:preBorn} back into itself, the scattered wavefunction is
\begin{align} \psi_s(\mathbf{r}) &= \int d\mathbf{r}_1\;\hat{G}_0(\mathbf{r},\mathbf{r}_1)\hat{V}(\mathbf{r}_1)\;\psi(\mathbf{r}_1) \nonumber\\ &= \int d\mathbf{r}_1\;\hat{G}_0(\mathbf{r},\mathbf{r}_1)\hat{V}(\mathbf{r}_1)\;\psi_i(\mathbf{r}_1) + \int d\mathbf{r}_1 \int d\mathbf{r}_2 \;\hat{G}_0(\mathbf{r},\mathbf{r}_1)\hat{V}(\mathbf{r}_1)\;\hat{G}_0(\mathbf{r}_1,\mathbf{r}_2)\hat{V}(\mathbf{r}_2)\;\psi_i(\mathbf{r}_2) + … \nonumber\\&\quad+ \int d\mathbf{r}_1 \int d\mathbf{r}_2 \cdots \int d\mathbf{r}_n\; \hat{G}_0(\mathbf{r},\mathbf{r}_1)\hat{V}(\mathbf{r}_1)\;\hat{G}_0(\mathbf{r}_1,\mathbf{r}_2)\hat{V}(\mathbf{r}_2) \cdots \hat{G}_0(\mathbf{r}_{n-1},\mathbf{r}_n)\hat{V}(\mathbf{r}_n) \; \psi_i(\mathbf{r}_n) + … \end{align}
This is no different from the Born series, and we find that
\begin{equation} \psi_s(\mathbf{r}) = \sum_{n=1}^{\infty} \underbrace{\int d\mathbf{r}_0\;\mathbf{\delta}\left(\mathbf{r} - \mathbf{r}_0\right)\prod_{m=1}^{n}\left[\int d\mathbf{r}_m\; \hat{G}_0(\mathbf{r}_{m-1},\mathbf{r}_m)\hat{V}(\mathbf{r}_m)\right] \psi_i(\mathbf{r}_n)}_{\text{call this }\Psi_n(\mathbf{r})} = \sum_{n=1}^{\infty}\Psi_n(\mathbf{r})\label{eq:born} \end{equation}
Similar to the notes, we are interested in investigating the 1D potential barrier. We begin by writing the explicit forms of the equations we will need.
Unlike the general 3D case, only two anti-commutating matrices are needed. Therefore, we can stick to the Pauli matrices, and choose \(\hat{\alpha}_0 = \hat{\sigma}_3\), \(\hat{\alpha}_1 = \hat{\sigma}_1\). Therefore, the 1D Dirac Hamiltonian is
\begin{equation} \hat{H}_0 = \hat{\sigma}_3mc^2 + \hat{\sigma}_1\hat{p} = \begin{bmatrix} mc^2 & -i\hbar\partial_x \\ -i\hbar\partial_x & -mc^2 \end{bmatrix} \end{equation}
Using Equation \ref{eq:greenRel}, and substituting \(G_{0,s}(x, x')\) from the notes,
\begin{align} \hat{G}_0(x, x') &= \frac{1}{2ik\cdot\hbar^2c^2}\left[E\hat{I} + \hat{\sigma}_3mc^2 - \hat{\sigma}_1i\hbar c\partial_x \right]e^{ik|x-x'|}\hat{I} \nonumber\\[10pt] &= \frac{e^{ik|x-x'|}}{2ik\cdot\hbar^2c^2}\left[E\hat{I} + \hat{\sigma}_3mc^2 + \hat{\sigma}_1\hbar kc\,\partial_x|x-x'| \right] \nonumber\\[10pt] &= \frac{e^{ik|x-x'|}}{2ik\cdot\hbar^2c^2}\left[E\hat{I} + \hat{\sigma}_3mc^2 + \hat{\sigma}_1\sgn(x-x')\,\hbar kc \right]\nonumber\\[10pt] &= \frac{e^{ik|x-x'|}}{2i\cdot\hbar c} \begin{bmatrix} \frac{E + mc^2}{\hbar kc} & \sgn(x-x')\\ \sgn(x-x') & \frac{E - mc^2}{\hbar kc} \end{bmatrix} \end{align}
Far from the scatterer, we use the approximation \(e^{ik|x-x'|} \approx e^{ik|x|}e^{-ik\sgn(x)x'}\)† for \(x\gg x'\) and \(\sgn(x-x') = \sgn(x)\) for \(|x|>|x'|\) to get
\begin{equation} \hat{G}_0(x, x') \approx \frac{e^{ik\left(|x| - \sgn(x)x'\right)}}{2i\cdot\hbar c} \begin{bmatrix} \frac{E + mc^2}{\hbar kc} & \sgn(x)\\ \sgn(x) & \frac{E - mc^2}{\hbar kc} \end{bmatrix} \end{equation}
Notice that apart from the oscillatory factor, the only terms with any dependence on \(x\) are factors of \(\sgn(x)\). This means that in regions where \(x\) is always either positive or negative, we can write
\begin{equation} \hat{G}_0(x, x'; \pm x>0) = \hat{G}_0(0_\pm, x')\,e^{ik|x|}\label{eq:1Dg} \end{equation}
This can be further applied to Equation \ref{eq:born}
\begin{align} \Psi_n(x;\pm x>0) &= \int dx_0\;\delta\left(x-x_0\right)\prod_{m=1}^{n} \left[\int dx_m\; \hat{G}_0(x_{m-1},x_m)\hat{V}(x_m)\right] \psi_i(x_n)\nonumber\\ &= e^{ik|x|}\;\int dx_1\;\hat{G}_0(0_\pm,x_1)\hat{V}(x_1)\; \prod_{m=2}^{n} \left[\int dx_m\; \hat{G}_0(x_{m-1},x_m)\hat{V}(x_m)\right] \psi_i(x_n)\nonumber\\ &= \Psi_n(0_\pm)\;e^{ik|x|}\nonumber\\[1ex] \therefore \psi_s(x;\pm x > 0) &= \sum_{n=1}^{\infty}\left[\Psi_n(0_\pm)\right]e^{ik|x|}\label{eq:1Dborn} \end{align}
To avoid complications in this matter, we will stick to the situation where a positive energy particle incident from the left scatters off a potential.
Consider an incident particle from the left:
\begin{equation} \psi_i = \phi\,e^{ikx} = \frac{1}{\sqrt{E+mc^2}}\begin{bmatrix}E+mc^2\\\hbar kc\end{bmatrix}e^{ikx} \end{equation}
note that it is normalised§ so that \({|\phi|}^2=2E\). we can write our wavefunction in terms of the reflection (\(r\)) and transmission (\(t\)) coefficients:
\begin{equation} \psi(x) = \psi_i(x) + \psi_s(x) = \begin{cases} \phi\,e^{ikx} + r\phi\,e^{-ikx} & x < 0\\[1ex] t\phi\,e^{ikx} & x > 0 \end{cases} \end{equation}
The reflectance (\(R\)) and transmittance (\(T\)) are then
\begin{align} R &= |r|^2 = {\left\lvert\frac{\psi_s(x;x<0)}{\phi\,e^{-ikx}}\right\rvert}^2 & T &= |t|^2 = {\left\lvert\frac{\psi_i(x) + \psi_s(x; x > 0)}{\phi\,e^{ikx}}\right\rvert}^2 \nonumber\\[7.5pt] &= \frac{1}{2E}\;{\left\lvert \left(\sum_{n=1}^{\infty}\Psi_n(0_-)\right)e^{-ikx} \right\rvert}^2 & &= \frac{1}{2E}\;{\left\lvert\phi\,e^{ikx} + \left(\sum_{n=1}^{\infty}\Psi_n(0_+)\right)e^{ikx}\right\rvert}^2\nonumber\\[7.5pt] &= \frac{1}{2E}\;{\left\lvert\sum_{n=1}^{\infty}\Psi_n(0_-)\right\rvert}^2 & &= \frac{1}{2E}\;{\left\lvert\phi + \sum_{n=1}^{\infty}\Psi_n(0_+)\right\rvert}^2 \label{eq:RT} \end{align}
Where we have used Equation \ref{eq:1Dborn} to extract the oscillatory dependence of the scattered wavefunction.
The expressions in Equation \ref{eq:RT} are a convenient form for computation as we can truncate the series by calculating \(\Psi_n(\mathbf{r})\) up to some finite \(n\) to get the \(n^{\text{th}}\) estimate of the reflectance and transmittance.
In order to perform numerical analysis and cross-check the validity of our Green's function, we will choose the step potential barrier \(\hat{V} = V_0\,\Theta(a-|x|)\) with the strength \(V_0\) and width \(2a\). Primarily this is chosen due to the existence of analytical solutions.[3]
Computational units \(\hbar = m = c = 1\) are used, and Equation \ref{eq:RT} is programatically calculated — see Appendix B for the source code. The analytical solutions are plotted with a solid line.
order estimate of reflectance and transmittance with barrier size \(a=0.05\)
We notice the same behaviour as in the notes, where the approximation is better with smaller potentials, and it converges into the true value as the energy becomes larger. Now, if we increase the barrier size to \(a=0.5\),
order estimate of reflectance and transmittance with barrier size \(a=0.5\)
The approximation becomes worse, and convergence takes longer (note the difference in the range of \(E\)). In fact, it does not seem to converge for \(V_0=1.75\). If we now look at the dependence of the approximation on \(a\):
order estimate of reflectance and transmittance energy \(E=1.5\) and barrier size \(ka \ll 1\)
order estimate of reflectance and transmittance energy \(E=1.5\) and barrier size \(ka \sim 1\)
It can be seen that the approximation holds for \(ka \ll 1\), but starts to diverge when \(ka\) is larger. Hence, the Born approximation tends to converge quickly with a large \(E\), a small \(V\), and a small \(ka\).
We conclude this report by noting the limitations of this approach: firstly, by focusing on elastic scattering, we have mechanisms like pair-production. This means that phenomena like the Klein paradox, where the potential is large enough to create electron-positron pairs, causing the reflectance and transmittance to deviate from the elastic scattering case.[4]
For an eigenstate of the free Dirac Hamiltonian, we find that it has the same form as the Schrödinger equation of an eigenstate (four Schrödinger equations for each component of \(\psi\)):
\begin{equation}i\hbar\partial_t\psi(\mathbf{r}, t) = \hat{H}_0\psi(\mathbf{r}, t) = E\psi(\mathbf{r}, t)\end{equation}
So, the wavefunction is similarly seperable into \(\psi(\mathbf{r}, t) = \psi(\mathbf{r}) e^{-i\frac{E}{\hbar}t}\) with \(\psi(\mathbf{r}) = \psi(\mathbf{r}, t=0)\).
Furthermore, the eigenstates of the Dirac Hamiltonian should also share a similar \(\mathbf{r}\)-dependence as the Schrödinger case, which means
\begin{equation} \psi(\mathbf{r}) = \phi e^{i\mathbf{k}\cdot\mathbf{r}} = \begin{bmatrix}\phi_+ \\ \phi_-\end{bmatrix}e^{i\mathbf{k}\cdot\mathbf{r}} \end{equation}
where \(\phi_+\) and \(\phi_-\) are 2-column vectors to be found. The action of the momentum operator on the eigenstate is \(\hat{\mathbf{p}}\psi(\mathbf{r}) = \hbar\mathbf{k}\psi(\mathbf{r})\), as expected.
The eigenequation is now
\begin{align} 0 = \left[\hat{H}_0 - E\hat{I}\right]\psi(\mathbf{r}) &= \begin{bmatrix} mc^2 - E & \hat{\boldsymbol{\sigma}}\cdot\hat{\mathbf{p}}c \\ \hat{\boldsymbol{\sigma}}\cdot\hat{\mathbf{p}}c & -mc^2-E \end{bmatrix} \begin{bmatrix} \phi_+ \\ \phi_- \end{bmatrix} e^{i\mathbf{k}\cdot\mathbf{r}} \nonumber \\[10pt] &= \begin{bmatrix} -\left(E-mc^2\right)\phi_+ + \hbar c\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_- \\ \hbar c\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_+ - \left(mc^2+E\right)\phi_- \end{bmatrix} e^{i\mathbf{k}\cdot\mathbf{r}} \end{align}
Hence we have
\begin{align} \phi_+ &= \frac{\hbar c}{E-mc^2}\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_-\quad(\text{if }E<0) & \phi_- &= \frac{\hbar c}{E+mc^2}\hat{\boldsymbol{\sigma}}\cdot\mathbf{k}\,\phi_+\quad(\text{if }E>0) \end{align}
The conditions are required to prevent singular values for the denominator \(E \pm mc^2\) since \(|E| = mc^2\) at the particle's rest frame where \(p=0\).
Denoting a right/left moving plane wave with positive energy as \(\phi e^{\pm ikx}\), and with our choice of \(\alpha_1 = \sigma_1\), Equation \ref{eq:PMrelation} is
\begin{equation} \phi = \begin{bmatrix}\phi_+\\\phi_-\end{bmatrix} = \phi_+\begin{bmatrix}1\\\frac{\hbar k c}{E+mc^2}\end{bmatrix} \end{equation}
while the normalisation of a plane wave \(|\psi(x)|^2 = 2E\) (by convention) gives
\begin{align} 2E &= {|\phi_+|}^2\left(1 + \frac{\hbar^2k^2c^2}{{\left(E+mc^2\right)}^2}\right)\nonumber\\ &= {|\phi_+|}^2\;\frac{E^2 + 2Emc^2 + \overbrace{m^2c^4 + \hbar^2k^2c^2}^{E^2}}{{\left(E+mc^2\right)}^2}\nonumber\\[10pt] &= {|\phi_+|}^2\;\frac{2E\cancel{\left(E + mc^2\right)}}{{\left(E+mc^2\right)}^{\cancel{2}}}\nonumber\\[10pt] \implies \phi_+ &= \sqrt{E^2 + mc^2}\\[10pt] \therefore \phi &= \frac{1}{\sqrt{E+mc^2}}\begin{bmatrix}E+mc^2\\\hbar kc\end{bmatrix} \end{align}
Our initial approach was to code Equation \ref{eq:RT} as-is: define the Green's function and initial wave function, and use scipy.integrate.nquad to perform the \(n^{\text{th}}\)-order integration. This approach is still given in the Numerical Integration Approach, with the function name RT_int.
While this worked, the computation time increased rapidly with each order. We came up with an alternative approach, where the integration domain is split according to its polarity, and the integral in each subdomain worked out analytically. The intermediate steps for this calculation is documented in Non-numerical Integration Of Reflectance/Transmittance, while the function RT_int2 is laid out in the Non-numerical Integration Approach. RT_int and RT_int2 are drop-in replacements of each other.
The transmittance and reflectance is the analytical solution given in the reference[5].
from scipy import *
# When having several arrays as arguments,
# resizes them so they have the same length
def uniformise(*args):
args = list(args)
lengths = unique([1 if isscalar(arg) else len(arg) for arg in args])
# Valid only if there are at most two unique lengths,
# one of them being 1
if len(lengths) <= 2 and min(lengths) == 1:
N = max(lengths)
for i in range(len(args)):
if isscalar(args[i]) or len(args[i]) == 1:
args[i] = repeat(args[i], N)
else:
raise ValueError("multiple independent variables")
return tuple(args)
# Define the Green's function
def G0(x0, x1, E):
# We hijack infinity as an infinitesimal,
# we we keep the sign, then set it to 0
sgn = sign(x0)
if isinf(x0):
x0 = 0
k = sqrt(E**2 - 1)
return exp(
1j*k*(abs(x0)-sgn*x1)
)/2j * array([
[(E+1)/k, sgn],
[sgn, (E-1)/k]
])
# Define the initial wave function
def psi_i(x, E):
k = sqrt(E**2 - 1)
return exp(1j*k*x)/sqrt(E+1)*array([E+1, k])
# Define the reflectance/transmittance
def RT_theory(E, V, a, isT = True):
E, V, a = uniformise(E, V, a)
K = sqrt((V-E+1)/(V-E-1)*(E+1)/(E-1))
A = 4*K**2
B = (1-K**2)**2*sin(2*sqrt( (V-E)**2 - 1 )*a)**2
return ( (A if isT else B)/(A + B) ).real
RT_int and RT_int2 In the definitions of both functions, we store the result of each integrand in the global variable calcedIntegrand. This is due to the repeated integrals in each order (the fourth order includes the third order integral includes the second order integral…).
from scipy.integrate import nquad
def RT_int(E, V, a, order=1, isT=True):
global calcedIntegrand
def integrand(*x):
# The argument comes in the form
# (x0, x1, x2,… xn, E, V, j, real, n),
# so we reassign them to forms that make sense
x = list(x)
n = x.pop()
if len(x) != n + 4:
raise ValueError("invalid number of arguments")
else:
real = x.pop()
j = x.pop()
V = x.pop()
E = x.pop()
# Include the "0th-order" infinitesimal term
x = [-inf, inf][isT] + x
# Since V(x) = V, we'll just stack the
# V's in front and the Green's functions together
Gtot = eye(2)
for k in range(n):
Gtot = Gtot @ G0(x[k], x[k+1], E)
_integrand = V**n * Gtot @ psi_i(x[n], E)
return _integrand[j].real if real else _integrand[j].imag
# Initialise the given arguments, and an empty array
E, V, a = uniformise(E, V, a)
RT = zeros(len(E))
# Loop through all unique values of E
for i in range(len(E)):
calcedIntegral = 0.
# nquad only allows functions with float (real) outputs,
# hence we need to split the wavefunction up into
# its two components (j=0,1), and further into its
# real and imaginary parts
for j in range(2):
for real in [True, False]:
integral = 0.
# There's an extra initial wavefunction to be
# included for the transmittance
if isT:
if real:
integral = psi_i(0., E[i])[j].real
else:
integral = psi_i(0., E[i])[j].imag
for n in arange(order)+1:
# Unique key for the inputs of the
# integrands we're calculating
keyIntegrand = f"{['R', 'T'][isT]}\
E{E[i]}V{V[i]}a{a[i]}n{n},{2*j + real}"
if not "calcedIntegrand" in globals():
calcedIntegrand = {}
# Proceed to do numerical integration
# only if a cached value doesn't exist
if not keyIntegrand in calcedIntegrand:
# The limits are (-a, a), (-a, a)…
# n times.
calcedIntegrand[keyIntegrand] = nquad(
integrand,
tile([-a[i], a[i]], [n,1]),
args = (
E[i],
V[i],
j,
real,
n
)
)[0]
integral += calcedIntegrand[keyIntegrand]
# \({|[a+ib, c+id]|}^2 = a^2 + b^2 + c^2 + d^2\), so we'll
# sum up the squares of the individual
# components (which we already know to
# be real since nquad only returns floats)
calcedIntegral += integral**2
RT[i] += calcedIntegral
return RT/E/2.
from functools import reduce
import itertools
def RT_int2(E, V, a, order=1, isT=True):
# Similar to above
global calcedIntegrand
E, V, a = uniformise(E, V, a)
k = sqrt(E**2 - 1.)
RT = zeros(len(E))
for i in range(len(E)):
calcedIntegral = psi_i(0., E[i]) if isT else 0.
for n in arange(order)+1:
keyIntegrand = f"{['R', 'T'][isT]}E{E[i]}V{V[i]}a{a[i]}n{n}"
if not "calcedIntegrand" in globals():
calcedIntegrand = {}
if not keyIntegrand in calcedIntegrand:
calcedIntegrand[keyIntegrand] = 0.
# Loop through all possible domains, given by \({\{+, -\}}^n\)
for sgnList in itertools.product([-1,1], repeat=n):
# Include the sign of the "0th-order" term
sgnList = ([-1,1][isT],) + sgnList
# Calculate the diff to determine the integrals,
# adjusting the last sign accordingly
sgnDiff = diff(sgnList)
sgnDiff[-1] = int(sgnList[-2] == -1)
# Where the signs are the same, we get \(a\)
# Where they are different, we get \(\frac{e^{2ika}-1}{2ik}\)
# Finally, reduce multiplies all the terms together
integrandTerms = reduce(
lambda x, y: x * y,
where(
sgnDiff == 0,
a[i],
(exp(2.j*k[i]*a[i]) - 1)/(2.j*k[i])
)
)
# Additional factor if last two signs are negative
if (sgnList[-2] == sgnList[-1] == -1):
integrandTerms *= exp(-2.j*k[i]*a[i])
#
# Calculate the matrix \(\prod_{i=0}^{n-1} G_0(0_{\sgn{(x_i)}}, 0)\)
#
integrandMatrices = reduce(
lambda x, y: x @ y,
[G0(sgn*inf, 0., E[i]) for sgn in sgnList[:-1]]
)
calcedIntegrand[keyIntegrand] += V[i]**n * integrandTerms \
* integrandMatrices @ psi_i(0., E[i])
calcedIntegral += calcedIntegrand[keyIntegrand]
RT[i] = dot(integratedTerms, integratedTerms.conj()).real
return RT/E/2.
The below examples uses RT_int2, but this can be simply swapped out with RT_int. However, the orders should probably be reduced greatly as RT_int is comically slower.
matplotlib Setup
import matplotlib.pyplot as plt
# As we're intending to plot three potentials on one graph
# and have the nth-order of each potential share the same colour
colours = ["#509abe", "#ea7643", "#c21236"]
plt.rcParams["mathtext.fontset"] = "dejavuserif"
# ordinal indicator for labelling the order
def th(x):
return "th" if x != int(x) or x < 0 or x%10 > 3 else ["st", "nd", "rd"][int(x%10-1)]
# Three groups of lines according to their potential
Vs = [.1, .4, 1.75]
# And two sets of \(a\) and \(\{E\}\)
# \(E=1\) causes a divide-by-zero situation,
# so we begin very slightly away from it
As = [.05, .5]
Eses = [ [1. + 1E-10, 1.1], [1. + 1E-10, 5] ]
# We plot order = {1, 3, 5}
orders = range(1, 6, 2)
# Produce a graph each for transmittance, reflectance,
# for each set of values of \(a\).
for isT in [False, True]:
for j in range(2):
a = As[j]
Es = Eses[j]
ylims = array([inf, -inf])
for i in range(len(Vs)):
V = Vs[i]
c = colours[i]
E = linspace(Es[0], Es[1], 500)
# Plot the analytical value, using it
# to demarcate the limits of the figure
RT = RT_theory(E, V, a, isT=isT)
ylims = array([
min(ylims[0], min(RT)),
max(ylims[1], max(RT))
])
plt.plot(
E,
RT,
label=f"\(V_0 = \){V}",
c = c
)
# Plot each order
for order in orders:
plt.plot(
E,
RT_int2(E, V, a, isT=isT, order=order),
label = (f"{order}{th(order)} order"),
c = c,
# If it is the highest order, we add
# markers to highlight its values
marker = "x" if order == max(orders) else "",
markevery = 50,
markeredgewidth = 1.5,
# Each order made to have \(n\) dots,
# and higher orer lines made to be thicker
linestyle = (0, (5, 1, *((1,1)*order))),
linewidth = (.5 + order/max(orders))/1.15
)
# Plot the limits with a bit of breather space
plt.ylim(ylims + array([-1., 1.]) * diff(ylims)/10.)
plt.title(f"\(a=\){a}")
plt.xlabel("\(E\)")
plt.ylabel(["\(R\)", "\(T\)"][isT])
# Since the linestyles are consistent among every potential,
# remove all the other orders from the legend apart from
# those belonging to the last potential to reduce clutter
handles, labels = plt.gca().get_legend_handles_labels()
handles = handles[0::len(orders)+1] + handles[-len(orders):]
labels = labels[0::len(orders)+1] + labels[-len(orders):]
legends = plt.legend(
handles,
labels,
fontsize = "small",
ncol = 1,
loc = ["upper right", "lower right"][isT],
bbox_to_anchor = [(0.975, 0.975), (0.975, 0.025)][isT],
frameon = False
).get_texts()
for legend in legends[0:len(Vs)]:
legend.set_fontproperties(legend.get_fontproperties())
legend.set_size("large")
plt.show()
In this approach, we calculate \(\Psi_n\) by splitting the limits of each integral into \((-a,0)\), \((0,a)\):
\begin{align} \Psi_n(0_{\pm}) &= \dots \times \int_{-a}^{a}dx_{m-1}\dots \times \int_{-a}^{a}dx_{m}G_0(x_{m-1}, x_m)V_0 \times \dots G_0(x_m, x_m+1) \dots \nonumber\\ &= \dots \times\int_{0}^{a}dx_{m-1}\dots \times\int_{0}^{a}dx_{m}G_0(x_{m-1}, x_m)V_0 \times \dots G_0(x_m, x_m+1) \dots \nonumber\\ &\hspace{.25em} + \dots \times\int_{0}^{a}dx_{m-1}\dots \times\int_{-a}^{0}dx_{m}G_0(x_{m-1}, x_m)V_0 \times \dots G_0(x_m, x_m+1) \dots \nonumber\\ &\hspace{.25em} + \dots \times\int_{-a}^{0}dx_{m-1}\dots \times\int_{0}^{a}dx_{m}G_0(x_{m-1}, x_m)V_0 \times \dots G_0(x_m, x_m+1) \dots \nonumber\\ &\hspace{.25em} + \dots \times\int_{-a}^{0}dx_{m-1}\dots \times\int_{-a}^{0}dx_{m}G_0(x_{m-1}, x_m)V_0 \times \dots G_0(x_m, x_m+1) \dots \end{align}
Effectively, splitting it in such a manner means that \(x_m\) has a consistent sign within each term. To evaluate each case, Equation \ref{eq:1Dg} can be first further simplified:
\begin{equation} \hat{G}_0(x, x'; \pm x>0) = e^{\pm ikx}\,\hat{G}_0(0_\pm, 0)\,e^{\mp ikx'} \end{equation}
Leaving \(\hat{G}_0(0_\pm, 0)\) as a matrix that doesn't depend on \(x\) or \(x'\). Finally, factoring out \(V_0\), then denoting the domain \(\int_0^a dx_m\) as \([\dots,\,p_m = +,\,\dots ]\) and \(\int_{-a}^0 dx_m\) as \([\dots,\,p_m = -,\,\dots ]\),
\begin{align} \text{Case }[\dots ,p_{m-1} = +, p_m = +,\dots ]:&\phantom{=}\;\;\dots \int_0^a dx_m\;G_0(x_{m-1},x_m)\dots G_0(x_m, x_{m+1})\dots \nonumber\\ &=\dots \int_0^a dx_m\;e^{ikx_{m-1}}\;G_0(0_+,0)e^{-ikx_m}\dots e^{ikx_m}\,G_0(0_+, x_{m+1})\dots \nonumber\\ &= \dots e^{ikx_{m-1}}\dots \int_0^a dx_m \; G_0(0_+,0)\dots \,G_0(0_+, x_{m+1})\dots \nonumber\\ &= \dots a \times G_0(0_+,0)\dots\\[10pt] \text{Case }[\dots ,p_{m-1} = -, p_m = -,\dots ]:&\phantom{=}\;\;\dots \int_{-a}^0 dx_m\;G_0(x_{m-1},x_m)\dots G_0(x_m, x_{m+1})\dots \nonumber\\ &= \dots e^{-ikx_{m-1}}\dots \int_{-a}^0 dx_m \; G_0(0_-,0)\dots \,G_0(0_-, x_{m+1})\dots \nonumber\\ &= \dots a \times G_0(0_-,0)\dots\\[10pt] \text{Case }[\dots ,p_{m-1} = +, p_m = -,\dots ]:&\phantom{=}\;\;\dots \int_{-a}^0 dx_m\;G_0(x_{m-1},x_m)\dots G_0(x_m, x_{m+1})\dots \nonumber\\ &=\dots \int_0^a dx_m\;e^{ikx_{m-1}}\;G_0(0_+,0)e^{-ikx_m}\dots e^{-ikx_m}\,G_0(0_-, x_{m+1})\dots \nonumber\\ &= \dots e^{ikx_{m-1}}\dots \int_{-a}^0 dx_m\;e^{-2ikx_m} \; G_0(0_+,0)\dots \,G_0(0_-, x_{m+1})\dots \nonumber\\ &= \dots \frac{e^{2ika}-1}{2ik} \times G_0(0_+,0)\dots\\[10pt] \text{Case }[\dots ,p_{m-1} = -, p_m = +,\dots ]:&\phantom{=}\;\;\dots \int_{0}^a dx_m\;G_0(x_{m-1},x_m)\dots G_0(x_m, x_{m+1})\dots \nonumber\\ &= \dots e^{-ikx_{m-1}}\dots \int_{0}^a dx_m\;e^{2ikx_m} \; G_0(0_-,0)\dots \,G_0(0_+, x_{m+1})\dots \nonumber\\ &= \dots \frac{e^{2ika}-1}{2ik} \times G_0(0_-,0)\dots \end{align}
Succinctly, we find:
\begin{equation} \begin{gathered} \text{Case }[\dots ,p_{m-1}, p_m,\dots ]:\phantom{=} \dots\times\mathcal{I}(p_{m-1}, p_m)\,G_0(0_{p_{m-1}},0)\times\dots\\[10pt] \text{where } \mathcal{I}(p_{m-1}, p_m) = \begin{cases} a & p_{m-1}\cdot p_m = 1\\[10pt] \frac{e^{2ika}-1}{2ik} & p_{m-1}\cdot p_m = -1 \end{cases} \end{gathered} \end{equation}
We take care to consider the edge case:
\begin{align} \text{Case }[\dots ,p_{n-1} = +, p_n = +]:&\phantom{=}\;\;\dots \int_0^a dx_n\;G_0(x_{n-1},x_n)\;\psi_i(x_n) \nonumber\\ &= \dots \int_0^a dx_n\;e^{ikx_{n-1}}\;G_0(0_+,0)e^{-ikx_n} \phi\,e^{ikx_n} \nonumber\\ &= \dots e^{ikx_{n-1}}\dots \int_0^a dx_n \; G_0(0_+,0)\,\phi\nonumber\\ &= \dots\mathcal{I}(+, +)\,G_0(0_+,0)\,\phi\\[10pt] \text{Case }[\dots ,p_{n-1} = +, p_n = -]:&\phantom{=}\;\;\dots \int_{-a}^0 dx_n\;G_0(x_{n-1},x_n)\;\psi_i(x_n) \nonumber\\ &= \dots \int_{-a}^0 dx_n\;e^{ikx_{n-1}}\;G_0(0_+,0)e^{-ikx_n} \phi\,e^{ikx_n} \nonumber\\ &= \dots\mathcal{I}(+, +)\,G_0(0_+,0)\,\phi\\[10pt] \text{Case }[\dots ,p_{n-1} = -, p_n = +]:&\phantom{=}\;\;\dots \int_0^a dx_n\;G_0(x_{n-1},x_n)\;\psi_i(x_n) \nonumber\\ &=\dots \int_0^a dx_n\;e^{-ikx_{n-1}}\;G_0(0_-,0)e^{ikx_n}\phi\,e^{ikx_n}\nonumber\\ &= \dots e^{-ikx_{m-1}}\dots \int_0^a dx_n\;e^{2ikx_n} \; G_0(0_-,0)\,\phi\nonumber\\ &= \dots\mathcal{I}(-, +)\,G_0(0_-,0)\,\phi\\[10pt] \text{Case }[\dots ,p_{n-1} = -, p_n = -]:&\phantom{=}\;\;\dots \int_{-a}^0 dx_n\;G_0(x_{n-1},x_n)\;\psi_i(x_n) \nonumber\\ &=\dots \int_{-a}^0 dx_n\;e^{-ikx_{n-1}}\;G_0(0_-,0)e^{ikx_n}\phi\,e^{ikx_n}\nonumber\\ &= \dots e^{-ikx_{m-1}}\dots \int_{-a}^0 dx_n\;e^{2ikx_n} \; G_0(0_-,0)\,\phi\nonumber\\ &= \dots\frac{1-e^{-2ika}}{2ik}\times G_0(0_-,0)\,\phi\nonumber\\ &= \dots e^{-2ika}\times\frac{e^{2ika}-1}{2ik}\times G_0(0_-,0)\,\phi \nonumber\\ &= \dots e^{-2ika}\times\mathcal{I}(-, +)\,G_0(0_-,0)\,\phi \end{align}
This is summarised by,
\begin{equation} \text{Case }[\dots ,p_{n-1}, p_n]:\;\;\dots \mathcal{I}(p_{n-1}, +)\,G_0(0_{p_{n-1}},0)\,\phi \times \left\lbrace e^{-2ika} \quad\text{if }p_{n-1}=p_n=-\right\rbrace \end{equation}
\begin{align} \therefore \Psi_n(0_{\pm}) = {V_0}^n\hspace{-1em}\sum_{(p_1, p_2, ...)\,\in\,{\left\{ -, + \right\}}^n}\!\!\!\!\!\! &\delta_{\pm,\,p_0}\prod_{m=1}^{n-1}\left[\mathcal{I}(p_{m-1}, p_m)\,G_0(0_{p_{m-1}}, 0)\right]\times\mathcal{I}(p_n, +)\,G_0(0_{p_{n-1}},0)\,\phi \nonumber\\ &\times \left\lbrace e^{-2ika} \quad\text{if }\,p_{n-1}=p_n=-\right\rbrace \end{align}