In this report, we are interested in extending the non-relativistic scattering theory into the relativistic realm using the Di­rac Ham­il­to­ni­an, building upon the material established in the course notes.^[2]

The Di­rac Ham­il­to­ni­an

The free, single-electron Di­rac Ham­il­to­ni­an is given by

$$\begin{array}{}\text{(1)}& {\hat{H}}_0={\hat{\alpha }}_{0}m{c}^2+\sum _{j=1}^3{\hat{\alpha }}_j{\hat{p}}_{j}c={\hat{\alpha }}_{0}m{c}^2+\hat{\mathit{\alpha }}\cdot \hat{\mathbf{p}}c\end{array}$$

where ${\alpha }_{\mu }$ are the standard Di­rac matrices, $\hat{\mathit{\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 Di­rac Ham­il­to­ni­an in the position representation, which will be useful later, is

$$\begin{array}{}\text{(2)}& {{\hat{H}}_0}^2=(m^2{c}^4+{\left|\hat{\mathbf{p}}\right|}^2{c}^2)\hat{I}=(m^2{c}^4-{\hslash }^2{c}^2{\mathrm{\nabla }}^2)\hat{I}\end{array}$$

Plane Wave Eigenstates

From here, operators are to be understood as given in the position representation. Considering the dimensionality of the Di­rac matrices, the wave­func­tion $\psi (\mathbf{r},t)$ is a 4-column vector.

The plane wave eigenstates are^§

$$\begin{array}{}\text{(3)}& \psi (\mathbf{r},t)=\varphi e^{i(\mathbf{k}\cdot \mathbf{r}-\omega t)}=\left[\begin{array}{c}{\varphi }_{+}\\ {\varphi }_{-}\end{array}\right]e^{i(\mathbf{k}\cdot \mathbf{r}-\omega t)}\end{array}$$

with $\omega =\frac{E}{\hslash }$ and

$$\begin{array}{}\text{(4)}& {\varphi }_{+}& =\frac{\hslash c}{E-m{c}^2}\hat{\mathit{\sigma }}\cdot \mathbf{k}{\varphi }_{-}(\text{if }E<0)& {\varphi }_{-}& =\frac{\hslash c}{E+m{c}^2}\hat{\mathit{\sigma }}\cdot \mathbf{k}{\varphi }_{+}(\text{if }E>0)\end{array}$$

Scattering

The scattering formulation from Chapter 0 of PH4401 still stands for the Di­rac Ham­il­to­ni­an without any modification needed, as no “illegal” moves have been made with regards to the wave­func­tion now being a four-component vector.

Therefore, given a Ham­il­to­ni­an $\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,

1. $\hat{H}\psi (\mathbf{r})=E\psi (\mathbf{r})$
2. $\psi (\mathbf{r})={\psi }_i(\mathbf{r})+{\psi }_s(\mathbf{r})$
3. ${\hat{H}}_0{\psi }_i(\mathbf{r})=E{\psi }_i(\mathbf{r})$
4. ${\psi }_s(\mathbf{r})$ must be “outgoing”

Green's Function for a Di­rac Particle

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ö­ding­er case is

$$\begin{array}{}\text{(5)}& [E_s-\frac{{\hat{\mathbf{p}}}^2}{2m}]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{2m}[p^2+{\hslash }^2{\mathrm{\nabla }}^2]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })={\delta }^d(\mathbf{r}-{\mathbf{r}}^{\prime })\end{array}$$

we can try swapping out the free-particle Ham­il­to­ni­an with the appropriate Di­rac operators:

$$\begin{array}{}\text{(6)}& [E\hat{I}-{\hat{H}}_0]{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })=[E\hat{I}-{\hat{\alpha }}_{0}m{c}^2-\hat{\mathit{\alpha }}\cdot \hat{\mathbf{p}}]{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })={\delta }^d(\mathbf{r}-{\mathbf{r}}^{\prime })\hat{I}\end{array}$$

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 Di­rac equation. That is, given ${\hat{G}}_0$ (possibly a 4-by-4 matrix) that satisfies the above equation,

$$\begin{array}{}\text{(7)}& \psi (\mathbf{r})={\psi }_i(\mathbf{r})+\underset{{\psi }_s(\mathbf{r})}{\underbrace{\int d{\mathbf{r}}^{\prime }{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })\hat{V}({\mathbf{r}}^{\prime })\psi ({\mathbf{r}}^{\prime })}}\end{array}$$

Instead of deriving the above, we will instead show that this is correct by applying the Ham­il­to­ni­an to both sides of the equation:

$$\begin{array}{rl}\hat{H}\psi (\mathbf{r})& =[{\hat{H}}_0+\hat{V}(\mathbf{r})]({\psi }_i(\mathbf{r})+\int d{\mathbf{r}}^{\prime }{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })\hat{V}({\mathbf{r}}^{\prime })\psi ({\mathbf{r}}^{\prime }))\\ E\psi (\mathbf{r})& =[E+\hat{V}(\mathbf{r})]{\psi }_i(\mathbf{r})+[{\hat{H}}_0-E]\int d{\mathbf{r}}^{\prime }{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })\hat{V}({\mathbf{r}}^{\prime })\psi ({\mathbf{r}}^{\prime })+(\hat{V}(\mathbf{r})+E){\psi }_s(\mathbf{r})\\ & =E({\psi }_i(\mathbf{r})+{\psi }_s(\mathbf{r}))+\hat{V}(\mathbf{r})({\psi }_i(\mathbf{r})+{\psi }_s(\mathbf{r}))-\int d{\mathbf{r}}^{\prime }\underset{{\delta }^d(\mathbf{r}-{\mathbf{r}}^{\prime })\hat{I}}{\underbrace{(E-{\hat{H}}_0){\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })}}\hat{V}({\mathbf{r}}^{\prime })\psi ({\mathbf{r}}^{\prime })\\ & =E\psi (\mathbf{r})+\hat{V}(\mathbf{r})\psi (\mathbf{r})-\hat{V}(\mathbf{r})\psi (\mathbf{r})\\ \text{(8)}& & =E\psi (\mathbf{r})\end{array}$$

Relation with the Schrö­ding­er Green's Function

As the one-, two-, and three-dimensional Green's functions in the Schrö­ding­er case have already been found in the course notes, it makes sense to equate Equation $\text{5}$ & Equation $\text{6}$ to relate them with the Di­rac Green's function so as to use our prior result.

$$\begin{array}{}\text{(9)}& [E\hat{I}-{\hat{H}}_0]{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{2m}[p^2+{\hslash }^2{\mathrm{\nabla }}^2]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{I}\end{array}$$

Since the relativistic energy of the incident particle is related to its momentum by

$$\begin{array}{}\text{(10)}& E^2=m^2{c}^4+p^2{c}^2⟹p^2=\frac{E^2-m^2{c}^4}{c^2}\end{array}$$

we place this back into Equation $\text{9}$.

$$\begin{array}{rl}[E\hat{I}-{\hat{H}}_0]{\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })& =\frac{1}{2m{c}^2}[E^2-m^2{c}^4+{\hslash }^2{c}^2{\mathrm{\nabla }}^2]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{I}\\ & =\frac{1}{2m{c}^2}[E^2\hat{I}-(m^2{c}^4-{\hslash }^2{c}^2{\mathrm{\nabla }}^2)\hat{I}]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{I}\\ & =\frac{1}{2m{c}^2}[{\left(E\hat{I}\right)}^2-{{\hat{H}}_0}^2]G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{I}\\ & =[E\hat{I}-{\hat{H}}_0]\left([E\hat{I}+{\hat{H}}_0]\frac{G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })}{2m{c}^2}\hat{I}\right)\\ \text{(11)}& \therefore {\hat{G}}_0(\mathbf{r},{\mathbf{r}}^{\prime })& =[E\hat{I}+{\hat{H}}_0]\frac{G_{0,s}(\mathbf{r},{\mathbf{r}}^{\prime })}{2m{c}^2}\hat{I}\end{array}$$

Here Equation $\text{2}$ is used in the second line, and the fact that $E\hat{I}$ commutes with ${\hat{H}}_0$ is used in the third.

Born series expansion

By repeatedly placing Equation $\text{7}$ back into itself, the scattered wave­func­tion is

$$\begin{array}{rl}{\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)\\ & =\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)+\dots \\ \text{(12)}& & +\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)+\dots \end{array}$$

This is no different from the Born series, and we find that

$$\begin{array}{}\text{(13)}& {\psi }_s(\mathbf{r})=\sum _{n=1}^{\mathrm{\infty }}\underset{\text{call this }{\mathrm{\Psi }}_n(\mathbf{r})}{\underbrace{\int d{\mathbf{r}}_0\delta (\mathbf{r}-{\mathbf{r}}_0)\prod _{m=1}^n[\int d{\mathbf{r}}_m{\hat{G}}_0({\mathbf{r}}_{m-1},{\mathbf{r}}_m)\hat{V}({\mathbf{r}}_m)]{\psi }_i({\mathbf{r}}_n)}}=\sum _{n=1}^{\mathrm{\infty }}{\mathrm{\Psi }}_n(\mathbf{r})\end{array}$$

1D Scattering

The 1D Equations

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.

Di­rac Ham­il­to­ni­an

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 Di­rac Ham­il­to­ni­an is

$$\begin{array}{}\text{(14)}& {\hat{H}}_0={\hat{\sigma }}_{3}m{c}^2+{\hat{\sigma }}_1\hat{p}=\left[\begin{array}{cc}m{c}^2& -i\hslash {\partial }_x\\ -i\hslash {\partial }_x& -m{c}^2\end{array}\right]\end{array}$$

Green's function

Using Equation $\text{11}$, and substituting $G_{0,s}(x,x^{\prime })$ from the notes,

$$\begin{array}{rl}{\hat{G}}_0(x,x^{\prime })& =\frac{1}{2ik\cdot {\hslash }^2{c}^2}[E\hat{I}+{\hat{\sigma }}_{3}m{c}^2-{\hat{\sigma }}_{1}i\hslash c{\partial }_x]e^{ik|x-x^{\prime }|}\hat{I}\\ & =\frac{e^{ik|x-x^{\prime }|}}{2ik\cdot {\hslash }^2{c}^2}[E\hat{I}+{\hat{\sigma }}_{3}m{c}^2+{\hat{\sigma }}_1\hslash kc{\partial }_x|x-x^{\prime }|]\\ & =\frac{e^{ik|x-x^{\prime }|}}{2ik\cdot {\hslash }^2{c}^2}[E\hat{I}+{\hat{\sigma }}_{3}m{c}^2+{\hat{\sigma }}_1\text{sgn}(x-x^{\prime })\hslash kc]\\ \text{(15)}& & =\frac{e^{ik|x-x^{\prime }|}}{2i\cdot \hslash c}\left[\begin{array}{cc}\frac{E+m{c}^2}{\hslash kc}& \text{sgn}(x-x^{\prime })\\ \text{sgn}(x-x^{\prime })& \frac{E-m{c}^2}{\hslash kc}\end{array}\right]\end{array}$$

Far from the scatterer, we use the approximation $e^{ik|x-x^{\prime }|}\approx e^{ik|x|}{e}^{-ik\text{sgn}(x)x^{\prime }}$^† for $x\gg x^{\prime }$ and $\text{sgn}(x-x^{\prime })=\text{sgn}(x)$ for $|x|>|x^{\prime }|$ to get

$$\begin{array}{}\text{(16)}& {\hat{G}}_0(x,x^{\prime })\approx \frac{e^{ik(|x|-\text{sgn}(x)x^{\prime })}}{2i\cdot \hslash c}\left[\begin{array}{cc}\frac{E+m{c}^2}{\hslash kc}& \text{sgn}(x)\\ \text{sgn}(x)& \frac{E-m{c}^2}{\hslash kc}\end{array}\right]\end{array}$$

Notice that apart from the oscillatory factor, the only terms with any dependence on $x$ are factors of $\text{sgn}(x)$. This means that in regions where $x$ is always either positive or negative, we can write

$$\begin{array}{}\text{(17)}& {\hat{G}}_0(x,x^{\prime };\pm x>0)={\hat{G}}_0(0_{\pm },x^{\prime })e^{ik|x|}\end{array}$$

This can be further applied to Equation $\text{13}$

$$\begin{array}{rl}{\mathrm{\Psi }}_n(x;\pm x>0)& =\int d{x}_0\delta (x-x_0)\prod _{m=1}^n[\int d{x}_m{\hat{G}}_0(x_{m-1},x_m)\hat{V}(x_m)]{\psi }_i(x_n)\\ & =e^{ik|x|}\int d{x}_1{\hat{G}}_0(0_{\pm },x_1)\hat{V}(x_1)\prod _{m=2}^n[\int d{x}_m{\hat{G}}_0(x_{m-1},x_m)\hat{V}(x_m)]{\psi }_i(x_n)\\ & ={\mathrm{\Psi }}_n(0_{\pm })e^{ik|x|}\\ \text{(18)}& \therefore {\psi }_s(x;\pm x>0)& =\sum _{n=1}^{\mathrm{\infty }}[{\mathrm{\Psi }}_n(0_{\pm })]e^{ik|x|}\end{array}$$

Setting Up the Scattering Problem

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.

Reflection and Transmission Coefficients

Consider an incident particle from the left:

$$\begin{array}{}\text{(19)}& {\psi }_i=\varphi e^{ikx}=\frac{1}{\sqrt{E+m{c}^2}}\left[\begin{array}{c}E+m{c}^2\\ \hslash kc\end{array}\right]e^{ikx}\end{array}$$

note that it is normalised^§ so that ${|\varphi |}^2=2E$. we can write our wave­func­tion in terms of the reflection ($r$) and transmission ($t$) coefficients:

$$\begin{array}{}\text{(20)}& \psi (x)={\psi }_i(x)+{\psi }_s(x)=\{\begin{array}{ll}\varphi e^{ikx}+r\varphi e^{-ikx}& x<0\\ t\varphi e^{ikx}& x>0\end{array}\end{array}$$

The reflectance ($R$) and transmittance ($T$) are then

$$\begin{array}{rlrl}R& =|r{|}^2={\left|\frac{{\psi }_s(x;x<0)}{\varphi e^{-ikx}}\right|}^2& T& =|t{|}^2={\left|\frac{{\psi }_i(x)+{\psi }_s(x;x>0)}{\varphi e^{ikx}}\right|}^2\\ & =\frac{1}{2E}{\left|(\sum _{n=1}^{\mathrm{\infty }}{\mathrm{\Psi }}_n(0_{-}))e^{-ikx}\right|}^2& & =\frac{1}{2E}{|\varphi e^{ikx}+(\sum _{n=1}^{\mathrm{\infty }}{\mathrm{\Psi }}_n(0_{+}))e^{ikx}|}^2\\ \text{(21)}& & =\frac{1}{2E}{|\sum _{n=1}^{\mathrm{\infty }}{\mathrm{\Psi }}_n(0_{-})|}^2& & =\frac{1}{2E}{|\varphi +\sum _{n=1}^{\mathrm{\infty }}{\mathrm{\Psi }}_n(0_{+})|}^2\end{array}$$

Where we have used Equation $\text{18}$ to extract the oscillatory dependence of the scattered wave­func­tion.

The expressions in Equation $\text{21}$ are a convenient form for computation as we can truncate the series by calculating ${\mathrm{\Psi }}_n(\mathbf{r})$ up to some finite $n$ to get the $n^{\text{th}}$ estimate of the reflectance and transmittance.

Numerical Analysis

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\mathrm{\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 $\hslash =m=c=1$ are used, and Equation $\text{21}$ is programatically calculated — see Appendix B for the source code. The analytical solutions are plotted with a solid line.

1st2nd3rd4th5th6th7th8th9th10th11th12th13th14th15th16th17th18th19th20th 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$,

1st2nd3rd4th5th6th7th8th9th10th11th12th13th14th15th16th17th18th19th20th 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$:

1st2nd3rd4th5th6th7th8th9th10th11th12th13th14th15th16th17th18th19th20th order estimate of reflectance and transmittance energy $E=1.5$ and barrier size $ka\ll 1$

1st2nd3rd4th5th6th7th8th9th10th11th12th13th14th15th16th17th18th19th20th 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]

Programming Approach

Our initial approach was to code Equation $\text{21}$ 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.

General Setup

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

Definition of 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…).

Numerical Integration Approach

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.

Non-numerical Integration Approach

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.

Plotting

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)]

Energy as an independent variable

# 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()

Non-numerical Integration Of Reflectance/Transmittance

In this approach, we calculate ${\mathrm{\Psi }}_n$ by splitting the limits of each integral into $(-a,0)$, $(0,a)$:

$$\begin{array}{rl}{\mathrm{\Psi }}_n(0_{\pm })& =\cdots \times {\int }_{-a}^{a}d{x}_{m-1}\cdots \times {\int }_{-a}^{a}d{x}_m{G}_0(x_{m-1},x_m)V_0\times \dots G_0(x_m,x_m+1)\dots \\ & =\cdots \times {\int }_0^{a}d{x}_{m-1}\cdots \times {\int }_0^{a}d{x}_m{G}_0(x_{m-1},x_m)V_0\times \dots G_0(x_m,x_m+1)\dots \\ & +\cdots \times {\int }_0^{a}d{x}_{m-1}\cdots \times {\int }_{-a}^{0}d{x}_m{G}_0(x_{m-1},x_m)V_0\times \dots G_0(x_m,x_m+1)\dots \\ & +\cdots \times {\int }_{-a}^{0}d{x}_{m-1}\cdots \times {\int }_0^{a}d{x}_m{G}_0(x_{m-1},x_m)V_0\times \dots G_0(x_m,x_m+1)\dots \\ \text{(29)}& & +\cdots \times {\int }_{-a}^{0}d{x}_{m-1}\cdots \times {\int }_{-a}^{0}d{x}_m{G}_0(x_{m-1},x_m)V_0\times \dots G_0(x_m,x_m+1)\dots \end{array}$$

Effectively, splitting it in such a manner means that $x_m$ has a consistent sign within each term. To evaluate each case, Equation $\text{17}$ can be first further simplified:

$$\begin{array}{}\text{(30)}& {\hat{G}}_0(x,x^{\prime };\pm x>0)=e^{\pm ikx}{\hat{G}}_0(0_{\pm },0)e^{\mp ik{x}^{\prime }}\end{array}$$

Leaving ${\hat{G}}_0(0_{\pm },0)$ as a matrix that doesn't depend on $x$ or $x^{\prime }$. Finally, factoring out $V_0$, then denoting the domain ${\int }_0^{a}d{x}_m$ as $[\dots ,p_m=+,\dots ]$ and ${\int }_{-a}^{0}d{x}_m$ as $[\dots ,p_m=-,\dots ]$,

$$\begin{array}{rl}\text{Case }[\dots ,p_{m-1}=+,p_m=+,\dots ]:& \phantom{=}\dots {\int }_0^{a}d{x}_m{G}_0(x_{m-1},x_m)\dots G_0(x_m,x_{m+1})\dots \\ & =\dots {\int }_0^{a}d{x}_m{e}^{ik{x}_{m-1}}{G}_0(0_{+},0)e^{-ik{x}_m}\dots e^{ik{x}_m}{G}_0(0_{+},x_{m+1})\dots \\ & =\dots e^{ik{x}_{m-1}}\dots {\int }_0^{a}d{x}_m{G}_0(0_{+},0)\dots G_0(0_{+},x_{m+1})\dots \\ \text{(31)}& & =\dots a\times G_0(0_{+},0)\dots \text{Case }[\dots ,p_{m-1}=-,p_m=-,\dots ]:& \phantom{=}\dots {\int }_{-a}^{0}d{x}_m{G}_0(x_{m-1},x_m)\dots G_0(x_m,x_{m+1})\dots \\ & =\dots e^{-ik{x}_{m-1}}\dots {\int }_{-a}^{0}d{x}_m{G}_0(0_{-},0)\dots G_0(0_{-},x_{m+1})\dots \\ \text{(32)}& & =\dots a\times G_0(0_{-},0)\dots \text{Case }[\dots ,p_{m-1}=+,p_m=-,\dots ]:& \phantom{=}\dots {\int }_{-a}^{0}d{x}_m{G}_0(x_{m-1},x_m)\dots G_0(x_m,x_{m+1})\dots \\ & =\dots {\int }_0^{a}d{x}_m{e}^{ik{x}_{m-1}}{G}_0(0_{+},0)e^{-ik{x}_m}\dots e^{-ik{x}_m}{G}_0(0_{-},x_{m+1})\dots \\ & =\dots e^{ik{x}_{m-1}}\dots {\int }_{-a}^{0}d{x}_m{e}^{-2ik{x}_m}{G}_0(0_{+},0)\dots G_0(0_{-},x_{m+1})\dots \\ \text{(33)}& & =\dots \frac{e^{2ika}-1}{2ik}\times G_0(0_{+},0)\dots \text{Case }[\dots ,p_{m-1}=-,p_m=+,\dots ]:& \phantom{=}\dots {\int }_0^{a}d{x}_m{G}_0(x_{m-1},x_m)\dots G_0(x_m,x_{m+1})\dots \\ & =\dots e^{-ik{x}_{m-1}}\dots {\int }_0^{a}d{x}_m{e}^{2ik{x}_m}{G}_0(0_{-},0)\dots G_0(0_{+},x_{m+1})\dots \\ \text{(34)}& & =\dots \frac{e^{2ika}-1}{2ik}\times G_0(0_{-},0)\dots \end{array}$$

Succinctly, we find:

$$\begin{array}{}\text{(35)}& \begin{array}{c}\text{Case }[\dots ,p_{m-1},p_m,\dots ]:\phantom{=}\cdots \times \mathcal{I}(p_{m-1},p_m)G_0(0_{p_{m-1}},0)\times \dots \\ \text{where }\mathcal{I}(p_{m-1},p_m)=\{\begin{array}{ll}a& p_{m-1}\cdot p_m=1\\ \frac{e^{2ika}-1}{2ik}& p_{m-1}\cdot p_m=-1\end{array}\end{array}\end{array}$$

We take care to consider the edge case:

$$\begin{array}{rl}\text{Case }[\dots ,p_{n-1}=+,p_n=+]:& \phantom{=}\dots {\int }_0^{a}d{x}_n{G}_0(x_{n-1},x_n){\psi }_i(x_n)\\ & =\dots {\int }_0^{a}d{x}_n{e}^{ik{x}_{n-1}}{G}_0(0_{+},0)e^{-ik{x}_n}\varphi e^{ik{x}_n}\\ & =\dots e^{ik{x}_{n-1}}\dots {\int }_0^{a}d{x}_n{G}_0(0_{+},0)\varphi \\ \text{(36)}& & =\dots \mathcal{I}(+,+)G_0(0_{+},0)\varphi \text{Case }[\dots ,p_{n-1}=+,p_n=-]:& \phantom{=}\dots {\int }_{-a}^{0}d{x}_n{G}_0(x_{n-1},x_n){\psi }_i(x_n)\\ & =\dots {\int }_{-a}^{0}d{x}_n{e}^{ik{x}_{n-1}}{G}_0(0_{+},0)e^{-ik{x}_n}\varphi e^{ik{x}_n}\\ \text{(37)}& & =\dots \mathcal{I}(+,+)G_0(0_{+},0)\varphi \text{Case }[\dots ,p_{n-1}=-,p_n=+]:& \phantom{=}\dots {\int }_0^{a}d{x}_n{G}_0(x_{n-1},x_n){\psi }_i(x_n)\\ & =\dots {\int }_0^{a}d{x}_n{e}^{-ik{x}_{n-1}}{G}_0(0_{-},0)e^{ik{x}_n}\varphi e^{ik{x}_n}\\ & =\dots e^{-ik{x}_{m-1}}\dots {\int }_0^{a}d{x}_n{e}^{2ik{x}_n}{G}_0(0_{-},0)\varphi \\ \text{(38)}& & =\dots \mathcal{I}(-,+)G_0(0_{-},0)\varphi \text{Case }[\dots ,p_{n-1}=-,p_n=-]:& \phantom{=}\dots {\int }_{-a}^{0}d{x}_n{G}_0(x_{n-1},x_n){\psi }_i(x_n)\\ & =\dots {\int }_{-a}^{0}d{x}_n{e}^{-ik{x}_{n-1}}{G}_0(0_{-},0)e^{ik{x}_n}\varphi e^{ik{x}_n}\\ & =\dots e^{-ik{x}_{m-1}}\dots {\int }_{-a}^{0}d{x}_n{e}^{2ik{x}_n}{G}_0(0_{-},0)\varphi \\ & =\dots \frac{1-e^{-2ika}}{2ik}\times G_0(0_{-},0)\varphi \\ & =\dots e^{-2ika}\times \frac{e^{2ika}-1}{2ik}\times G_0(0_{-},0)\varphi \\ \text{(39)}& & =\dots e^{-2ika}\times \mathcal{I}(-,+)G_0(0_{-},0)\varphi \end{array}$$

This is summarised by,

$$\begin{array}{}\text{(40)}& \text{Case }[\dots ,p_{n-1},p_n]:\dots \mathcal{I}(p_{n-1},+)G_0(0_{p_{n-1}},0)\varphi \times \{e^{-2ika}\text{if }{p}_{n-1}=p_n=-\}\end{array}$$

$$\begin{array}{rl}\therefore {\mathrm{\Psi }}_n(0_{\pm })={V_0}^n\sum _{(p_1,p_2,...)\in {\{-,+\}}^n}& {\delta }_{\pm ,p_0}\prod _{m=1}^{n-1}[\mathcal{I}(p_{m-1},p_m)G_0(0_{p_{m-1}},0)]\times \mathcal{I}(p_n,+)G_0(0_{p_{n-1}},0)\varphi \\ \text{(41)}& & \times \{e^{-2ika}\text{if }{p}_{n-1}=p_n=-\}\end{array}$$