<figure><img src="image_1.png" /></figure>

HAL
open science

# Foldy-Wouthuysen WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap

Mark Behzad Doost

► To cite this version:

Mark Behzad Doost. Foldy-Wouthuysen WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap. 2022. hal-01389947v21

HAL Id: hal-01389947

https://hal.science/hal-01389947v21

Preprint submitted on 10 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

HAL
open science

# Foldy-Wouthuysen Green's function and WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap

Mark Behzad Doost

► To cite this version:

Mark Behzad Doost. Foldy-Wouthuysen Green's function and WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap. 2024. hal-01389947v23

HAL Id: hal-01389947

https://hal.science/hal-01389947v23

Preprint submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Foldy-Wouthuysen WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap

Mark Behzad Doost<sup>a,*</sup>

<sup>a</sup>Independent Researcher

## ARTICLE INFO

Keywords:

Dirac equation
Foldy-Wouthuysen transformation
WKB approximation
WKB connection formulae
Band engineering
Graphene

## ABSTRACT

I provide a transfer matrix method for the Foldy-Wouthuysen representation of the Dirac equation. I develop the relationship between the reflection and transmission coefficients of the Dirac spinors and the wavefunction in the transformed representation. I develop a WKB approximation for Dirac fermions which has the same elegant form as the WKB solution to Schrödinger's equation. My WKB approximation is to all orders and includes the semi-classical turning point. I provide an extension to fully 2-dimensional periodic structures by Fourier methods for band engineering. I verify my methods for all energies by comparison with analytic solutions developed in the Dirac spinor representation. I also include a rich appendix detailing my research into the Green's functions of Dirac fermions.

## 1. Introduction

In Refs. [1, 2, 3] the pseudo-relativistic dispersion of 2-dimensional (2D) materials was investigated theoretically for the emergence of extra Dirac points with 1D periodic barriers. The method employed was the application of Bloch's theorem to transfer matrices of Dirac spinors. These band engineering studies attracted considerable interest since the vanishing band gap of 2D materials is an Achilles heel for applications. Other methods for band engineering that open the way to creating 2D semiconductor materials include strain engineering [4], vertical stacking [5], and chemical doping [6].

Armed with the knowledge to create semiconductor 2D materials a theoretical study of the operational window for negative differential resistance was pursued by other authors [7]. The mathematical method employed in their study was the application of transfer matrices to the Dirac spinors of a rotated Dirac Hamiltonian derived in Ref. [8].

Motivated by the study of unconventional transmission and density of states, the modeling of Dirac fermions has moved beyond unrealistic square barrier systems. In experimental devices barriers are more often smoothly varying, semi-classical (WKB) methods are needed. An approach using separate Hamiltonians for electrons and holes in Ref. [9] gave way in Ref. [10] to an expansion in powers of $\hbar$ partly resembling the WKB method for the Schrödinger equation. However it is acknowledged by the authors of Refs. [9, 10] that their mathematical approaches to semi-classical analysis can be divergent at the classical turning point.

The contribution to the theory of tunneling in 2D materials that I present in this article is the application of the Foldy-Wouthuysen (FW) representation of the Dirac equation.

Whilst the FW transformation was first conceived as a semi-relativistic representation of the Dirac equation [11], I will show, through my rigorous derivations in Section 2

to Section 13 and my analytic and numerical verification of Section 14 and Appendix B, that the FW representation remains an exactly accurate representation of the Dirac equation for all energies.

I find that the FW representation allows a derivation of the relativistic WKB approximation, Section 7 to Section 10, in the same elegant form as the approximate WKB solution to Schrödinger's equation. Since at low energies the FW representation reduces to the Schrödinger equation, I have found, in Section 11, a simple method to mathematically describe the connections between WKB regions.

Fourier analysis is used in Section 12 to apply my approach to fully 2D periodic structures. Section 12 offers the prospect of extending the scope of the band engineering studies made in Refs. [1, 2, 3]. My numerical and analytic results showing that the FW equation is an exactly accurate representation of relativistic fermions, for all energies, will be important here. The Fourier transform method in Section 12 requires expansion of the wavefunction in plane waves of momentum beyond semi-relativistic energies.

In Appendix E I make the reader aware of a limit to appealing for solutions of the Dirac equation from its FW representation, I reveal and discuss a Green's function (GF) paradox arising when transforming between the two representations.

My article is organized as follows: Section 2 outlines the derivation of the FW representation from the Dirac equation. Section 3 gives the transmission of the Dirac equation in terms of the transmission of its FW representation. Section 4 derives the boundary conditions of the FW equation at sharp steps in potential. Section 5 and Section 6 develop the transfer matrices for square barriers and delta potentials. Section 7, Section 8, Section 9 and Section 10 give a 1st, iterative and 2nd order WKB approximations for two interpretations of the FW equation. Section 11 gives the connecting formulae between regions where the WKB approximation is appropriate. Section 12 extends my approach to include periodic structures. Section 13 details how I evaluated the boundary conditions. Section 14 gives numerical results for

<figure><img src="image_2.png" /></figure>

<footnote>*Corresponding author</footnote>

<footnote>doostmb@gmail.com, https://gofund.me/34294631 (M.B. Doost)</footnote>

<footnote>ORCID(s): 0000-0003-4682-2889 (M.B. Doost)</footnote>
# Foldy-Wouthuysen Green's function and WKB transfer matrix method for Dirac tunneling through monolayer graphene with a mass gap

Mark Behzad Doost<sup>a,*</sup>

<sup>a</sup>Independent Researcher, Hullbridge, Essex, England

## ARTICLE INFO

**Keywords:**

Dirac equation

Foldy-Wouthuysen transformation

Green's functions

WKB approximation

WKB connection formulae

Band-gap engineering

Graphene

## ABSTRACT

I provide a transfer matrix method for the Foldy-Wouthuysen representation of the Dirac equation. I derive the relationship between the reflection and transmission coefficients of the Dirac spinors and the wavefunction in the transformed representation. I develop a WKB approximation for Dirac fermions that has the same elegant form as the WKB solution to Schrödinger's equation. My WKB approximation is to all orders and includes the semi-classical turning point. I provide an extension to fully 2-dimensional periodic structures by Fourier methods for band-gap engineering. I verify my methods for all energies by comparison with analytic solutions developed in the Dirac spinor representation. Rich appendices detail my research into the Green's functions of Dirac fermions, where I rigorously derive the free space Green's functions for the Foldy-Wouthuysen representation of the Dirac equation.

## 1. Introduction

In Refs. [1, 2, 3] the pseudo-relativistic dispersion of 2-dimensional (2D) materials was investigated theoretically for the emergence of extra Dirac points with 1D periodic barriers. The method employed was the application of Bloch's theorem to transfer matrices of Dirac spinors. These band-gap engineering studies have attracted considerable interest since the vanishing band gap of 2D materials presents a significant challenge for practical applications. Other methods for band-gap engineering used to create 2D semiconductor materials include strain engineering [4], vertical stacking [5], and chemical doping [6].

With the necessary knowledge for creating semiconductor 2D materials, a theoretical study of the operational window for negative differential resistance was pursued by other authors [7]. The mathematical method employed in their study was the application of transfer matrices to the Dirac spinors of a rotated Dirac Hamiltonian derived in Ref. [8].

Motivated by the study of unconventional transmission and density of states, the modeling of Dirac fermions has advanced beyond idealized square barrier systems. In experimental devices, barriers are more often smoothly varying; semi-classical (WKB) methods are needed. An approach using separate Hamiltonians for electrons and holes in Ref. [9] gave way in Ref. [10] to an expansion in powers of $ħ$ partly resembling the WKB method for the Schrödinger equation. However, it is acknowledged by the authors of Refs. [9, 10] that their mathematical approaches to semi-classical analysis can be divergent at the classical turning point.

This article contributes to the theory of tunneling in 2D materials by applying the Foldy-Wouthuysen (FW) representation of the Dirac equation.

While the FW transformation was first conceived as a semi-relativistic representation of the Dirac equation [11],

I will show, through my rigorous derivations in Section 2 through Section 13 and my analytic and numerical verification of Section 14 and Appendix B, that the FW representation remains an exactly accurate representation of the Dirac equation for all energies.

I find that the FW representation allows a derivation of the relativistic WKB approximation, Section 7 to Section 10, in the same elegant form as the approximate WKB solution to Schrödinger's equation. Since at low energies, the FW representation reduces to the Schrödinger equation, I have found, in Section 11, a simple method to mathematically describe the connections between WKB regions.

Fourier analysis is used in Section 12 to extend my approach to fully 2D periodic structures. Section 12 offers the prospect of extending the scope of the band-gap engineering studies made in Refs. [1, 2, 3]. My numerical and analytic results showing that the FW equation is an exactly accurate representation of relativistic fermions, for all energies, will be important here. The Fourier transform method in Section 12 requires expansion of the wavefunction in plane waves of momentum associated with beyond semi-relativistic energies.

In Appendix E I make the reader aware of a limit to appealing for solutions of the Dirac equation from its FW representation, I reveal and discuss a Green's function (GF) paradox arising when transforming between the two representations. In spite of the paradox, I have rigorously derived the 1D, 2D and 3D free space FW GFs in Appendix D and Appendix F. The Dyson equations for the FW GFs are available for all types of Born approximation [12, 13].

My article is organized as follows: Section 2 outlines the derivation of the FW representation from the Dirac equation. Section 3 gives the transmission of the Dirac equation in terms of the transmission of its FW representation. Section 4 derives the boundary conditions of the FW equation at sharp steps in potential. Section 5 and Section 6 develop the transfer matrices for rectangular barriers and delta potentials. Section 7, Section 8, Section 9 and Section 10 give

<figure><img src="image_3.png" /></figure>

<footnote>*Corresponding author</footnote>

<footnote>doostmb@gmail.com, https://gofund.me/34294631 (M.B. Doost)</footnote>

<footnote>ORCID(s): 0000-0003-4682-2889 (M.B. Doost)</footnote>
examples of tunneling through barriers and resonant diodes.

The appendices are organized as follows: Appendix A analytically justifies my commutation of FW wavefunction operators. Appendix B analytically calculates the reflection at a step in both the Dirac spinor and FW representations for comparison. Appendix C examines the tunneling of a massless fermion through a magnetic delta barrier in the FW representation. Appendix D calculates GFs of the FW equation. Appendix E uncovers a relativistic GF paradox.

## 2. Dirac equation and its Foldy-Wouthuysen representation

The Dirac equation for a scalar potential V is given in Hamiltonian form by

$$[c\alpha \cdot \hat{p} + \beta mc^2] \Psi = H\Psi = (E + V)\Psi, \quad (1)$$

$$\gamma_t = \begin{pmatrix} +I & 0 \\ 0 & -I \end{pmatrix}, \quad \gamma_{x,y,z} = \begin{pmatrix} 0 & +\sigma_{x,y,z} \\ -\sigma_{x,y,z} & 0 \end{pmatrix}, \quad (2)$$

$$I = \begin{pmatrix} +1 & 0 \\ 0 & +1 \end{pmatrix}, \quad \sigma_x = \begin{pmatrix} 0 & +1 \\ +1 & 0 \end{pmatrix}, \quad (3)$$

$$\sigma_y = \begin{pmatrix} 0 & -i \\ +i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} +1 & 0 \\ 0 & -1 \end{pmatrix}, \quad (4)$$

with $\alpha_{x,y,z} = \gamma_t \gamma_{x,y,z}$ and $\beta = \gamma_t$.

$$\hat{p} = -i\hbar\nabla = -i\hbar \left[ \frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} \right] \quad (5)$$

is the momentum operator for the electron which is of mass m and energy E. The charge and current density are given by the following well known expressions:

$$\rho = \Psi^* \Psi, \quad (6)$$

$$j = \Psi^* \alpha \Psi. \quad (7)$$

There are two linearly independent solutions of the free particle Dirac equation,

$$\Psi_{\uparrow}^{(\pm)} = \begin{pmatrix} \mp mc^2 - E - V \\ 0 \\ 0 \\ c(p_x \pm ip_y) \end{pmatrix} \exp \left( \frac{ip_x x + ip_y y}{\hbar} \right), \quad (8)$$

$$\Psi_{\downarrow}^{(\pm)} = \begin{pmatrix} 0 \\ \mp mc^2 - E - V \\ c(p_x \mp ip_y) \\ 0 \end{pmatrix} \exp \left( \frac{ip_x x + ip_y y}{\hbar} \right), \quad (9)$$

each can be normalised to one electron per unit volume by Eq. (6).

By making use of the FW transformation I will decouple the differential equations of the Dirac equation to give four independent identical equations. I will show that by calculating the propagation for only one of these components I may obtain the full GF and transmission.

Under unitary transformation the Dirac equation becomes

$$e^{+iS} H e^{-iS} e^{+iS} \Psi = (E + V) e^{+iS} \Psi. \quad (10)$$

Foldy and Wouthuysen [11] found that their rotation

$$e^{\pm iS} = \cos \theta \pm \beta \alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta, \quad (11)$$

transformed the Dirac Hamiltonian into diagonal and anti-diagonal parts

$$\beta (mc^2 \cos (2\theta) + c |\hat{p}| \sin (2\theta)) + c\alpha \cdot \hat{p} \left( \cos (2\theta) - \frac{mc}{|\hat{p}|} \sin (2\theta) \right) = e^{+iS} H e^{-iS}. \quad (12)$$

Subsequently it was chosen [11]

$$\cos (2\theta) = \frac{mc^2}{\sqrt{c^2 \hat{p}^2 + m^2 c^4}}, \quad (13)$$

$$\sin (2\theta) = \frac{c |\hat{p}|}{\sqrt{c^2 \hat{p}^2 + m^2 c^4}}, \quad (14)$$

so that by direct substitution of Eq. (13) and Eq. (14) into Eq. (12)

$$[\beta \sqrt{\hat{p}^2 c^2 + m^2 c^4} - \hat{1} (V + E)] \tilde{\Psi} = 0. \quad (15)$$

Eq. (15) shows there are two possibilities for the FW equation

$$[\mp \sqrt{\hat{p}^2 c^2 + m^2 c^4} - (V + E)] \tilde{\Psi}^{(\pm)} = 0, \quad (16)$$

$\tilde{\Psi}^{(-)}$ corresponds to the electron, $\tilde{\Psi}^{(+)}$ to the positron. The rotated wavefunction $\tilde{\Psi}$ is related to $\Psi$ by

$$[\cos \theta + \beta \alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta] \Psi = \tilde{\Psi}. \quad (17)$$

Taylor expanding $\sqrt{\hat{p}^2 c^2 + m^2 c^4}$ in Eq. (16) gives

$$[\mp \mathcal{L} - (V + E \pm mc^2)] \tilde{\Psi}^{(\pm)} = 0, \quad (18)$$

where

$$\mathcal{L} = mc^2 \left[ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 - \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 + \frac{1}{16} \left[ \frac{\hat{p}}{mc} \right]^6 - \dots \right]. \quad (19)$$

Eq. (16) and Eq. (18) are two interpretations of the FW equation.

I demonstrate the unsuitability of my FW approach for calculating magnetic barrier tunneling in Appendix C. I discuss the GFs of the FW Eq. (18) in Appendix D. I uncover and discuss a GF paradox of the FW Eq. (18) in Appendix E.
a 1st, iterative and 2nd order WKB approximations for two
interpretations of the FW equation. Section 11 gives the con-
necting formulae between regions where the WKB approx-
imation is appropriate. Section 12 extends my approach to
include periodic structures. Section 13 details how I eval-
uated the boundary conditions. Section 14 gives numerical
results for examples of tunneling through barriers and reso-
nant diodes.

The appendices are organized as follows: Appendix A analytically justifies my commutation of FW wavefunction operators. Appendix B analytically calculates the reflection at a step in both the Dirac spinor and FW representations for comparison. Appendix C examines the tunneling of a massless fermion through a magnetic delta barrier in the FW representation. Appendix D calculates 1D GFs of the FW equation. Appendix E uncovers a relativistic GF paradox. Appendix F gives the correct 2D and 3D free space analytic FW GFs.

## 2. Dirac equation and its Foldy-Wouthuysen representation

The Dirac equation for a scalar potential V is given in
Hamiltonian form by

$$
[c\alpha \cdot \hat{p} + \beta mc^2] \Psi = H\Psi = (E + V)\Psi, \quad (1)
$$

$$
\gamma_t = \begin{pmatrix} +I & 0 \\ 0 & -I \end{pmatrix}, \quad \gamma_{x,y,z} = \begin{pmatrix} 0 & +\sigma_{x,y,z} \\ -\sigma_{x,y,z} & 0 \end{pmatrix}, \quad (2)
$$

$$
I = \begin{pmatrix} +1 & 0 \\ 0 & +1 \end{pmatrix}, \quad \sigma_x = \begin{pmatrix} 0 & +1 \\ +1 & 0 \end{pmatrix}, \quad (3)
$$

$$
\sigma_y = \begin{pmatrix} 0 & -i \\ +i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} +1 & 0 \\ 0 & -1 \end{pmatrix}, \quad (4)
$$

with $\alpha_{x,y,z} = \gamma_t \gamma_{x,y,z}$ and $\beta = \gamma_t$.

$$
\hat{p} = -i\hbar\nabla = -i\hbar \left[ \frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} \right] \quad (5)
$$

is the momentum operator for the electron which is of mass
m and energy E. The charge and current density are given
by the following well-known expressions:

$$
\rho = \Psi^* \Psi, \tag{6}
$$

$$
j = \Psi^* \alpha \Psi. \tag{7}
$$

There are two linearly independent solutions of the free par-
ticle Dirac equation,

$$
\Psi_{\uparrow}^{(\pm)} = \begin{pmatrix} \mp mc^2 - E - V \\ 0 \\ 0 \\ c(p_x \pm ip_y) \end{pmatrix} \exp\left(\frac{ip_x x + ip_y y}{\hbar}\right), \quad (8)
$$

$$
\Psi_{\downarrow}^{(\pm)} = \begin{pmatrix} 0 \\ \mp mc^2 - E - V \\ c(p_x \mp ip_y) \\ 0 \end{pmatrix} \exp\left(\frac{ip_x x + ip_y y}{\hbar}\right), \quad (9)
$$

each can be normalized to one electron per unit volume by
Eq. (6).

By making use of the FW transformation, I will decouple
the differential equations of the Dirac equation to give four
independent identical equations. I will show that by calcu-
lating the propagation for only one of these components I
may obtain the full GF and transmission.

Under unitary transformation the Dirac equation becomes

$$
e^{+iS} H e^{-iS} e^{+iS} \Psi = (E + V) e^{+iS} \Psi. \quad (10)
$$

Foldy and Wouthuysen [11] found that their rotation

$$
e^{\pm iS} = \cos \theta \pm \beta \alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta, \quad (11)
$$

transformed the Dirac Hamiltonian into diagonal and anti-
diagonal parts

$$
\begin{align}
& \beta (mc^2 \cos (2\theta) + c |\hat{p}| \sin (2\theta)) \nonumber \\
& + c\alpha \cdot \hat{p} \left( \cos (2\theta) - \frac{mc}{|\hat{p}|} \sin (2\theta) \right) \nonumber \\
& = e^{+iS} H e^{-iS}. \tag{12}
\end{align}
$$

Subsequently it was chosen [11]

$$
\cos (2\theta) = \frac{mc^2}{\sqrt{c^2 \hat{p}^2 + m^2 c^4}}, \quad (13)
$$

$$
\sin (2\theta) = \frac{c |\hat{p}|}{\sqrt{c^2 \hat{p}^2 + m^2 c^4}}, \quad (14)
$$

so that by direct substitution of Eq. (13) and Eq. (14) into
Eq. (12)

$$
[\beta \sqrt{\hat{p}^2 c^2 + m^2 c^4} - \hat{1} (V + E)] \tilde{\Psi} = 0. \quad (15)
$$

Eq. (15) shows there are two possibilities for the FW equa-
tion

$$
[\mp \sqrt{\hat{p}^2 c^2 + m^2 c^4} - (V + E)] \tilde{\Psi}^{(\pm)} = 0, \quad (16)
$$

$\tilde{\Psi}^{(-)}$ corresponds to the electron, $\tilde{\Psi}^{(+)}$ to the positron. The rotated wavefunction $\tilde{\Psi}$ is related to $\Psi$ by

$$
\left[ \cos \theta + \beta \alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta \right] \Psi = \tilde{\Psi}. \qquad (17)
$$

Taylor expanding $\sqrt{\hat{p}^2 c^2 + m^2 c^4}$ in Eq. (16) gives

$$
[\mp \mathcal{L} - (V + E \pm mc^2)] \tilde{\Psi}^{(\pm)} = 0, \quad (18)
$$
### 3. Relationship between the transmission of the Dirac equation and its Foldy-Wouthuysen representation

Throughout my numerical and analytic results, Section 14 and Appendix B, I will calculate the reflection coefficient $R$ to deduce the transmission coefficient $T$ as

$$T = 1 - R. \qquad (20)$$

Calculating $R$, rather than $T$ directly, does not require the normalisation of the wavefunction in both the incident and transmitted regions.

Let $\Psi$ be given by Eq. (8) and Eq. (9) as

$$\Psi = u \exp \left( \frac{i p_x x + i p_y y}{\hbar} \right), \qquad (21)$$

taking $\tilde{\Psi}$ of the form

$$\tilde{\Psi} = \tilde{u} \exp \left( \frac{i p_x x + i p_y y}{\hbar} \right), \qquad (22)$$

I find by application of Eq. (17) to Eq. (21)

$$\begin{aligned}
\tilde{u}_1 &= u_1 \cos \theta + (d_x - id_y) u_4 \sin \theta, \\
\tilde{u}_2 &= u_2 \cos \theta + (d_x + id_y) u_3 \sin \theta, \\
\tilde{u}_3 &= u_3 \cos \theta - (d_x - id_y) u_2 \sin \theta, \\
\tilde{u}_4 &= u_4 \cos \theta - (d_x + id_y) u_1 \sin \theta,
\end{aligned}
\qquad (23)$$

where $\hat{d} = \hat{p}/|\hat{p}|$.

Eq. (21), Eq. (22), and Eq. (23) show $\tilde{\Psi}$ in exponential planar form which is useful for calculating the fermion propagation through homogeneous space.

The incident (i) and reflected (r) waves for planar structures in the Dirac equation description are related to the reflection coefficient in the well known way:

$$R = \frac{(\Psi \cdot \Psi^*)^{(r)}}{(\Psi \cdot \Psi^*)^{(i)}}. \qquad (24)$$

For the two linearly independent solutions of the Dirac equation, Eq. (8) and Eq. (9), Eq. (24) gives rise to two degenerate expressions for the reflection coefficient:

$$R_{\uparrow} = \frac{(|\Psi_1|^2 + |\Psi_4|^2)^{(r)}}{(|\Psi_1|^2 + |\Psi_4|^2)^{(i)}}, \qquad (25)$$

$$R_{\downarrow} = \frac{(|\Psi_2|^2 + |\Psi_3|^2)^{(r)}}{(|\Psi_2|^2 + |\Psi_3|^2)^{(i)}}. \qquad (26)$$

I take note of Eqs. (23) linking $\tilde{u}$ to $u$ and calculate when $d$ is real

$$|\tilde{u}_1|^2 + |\tilde{u}_4|^2 = |u_1|^2 + |u_4|^2, \qquad (27)$$

$$|\tilde{u}_2|^2 + |\tilde{u}_3|^2 = |u_2|^2 + |u_3|^2, \qquad (28)$$

therefore $R$ described by the Dirac equation can be rewritten in terms of the FW wavefunctions:

$$R_{\uparrow} = \frac{\left( |\tilde{\Psi}_1|^2 + |\tilde{\Psi}_4|^2 \right)^{(r)}}{\left( |\tilde{\Psi}_1|^2 + |\tilde{\Psi}_4|^2 \right)^{(i)}}, \qquad (29)$$

$$R_{\downarrow} = \frac{\left( |\tilde{\Psi}_2|^2 + |\tilde{\Psi}_3|^2 \right)^{(r)}}{\left( |\tilde{\Psi}_2|^2 + |\tilde{\Psi}_3|^2 \right)^{(i)}}. \qquad (30)$$

Since the components of $\tilde{\Psi}$ are degenerate I am only required to evaluate the propagation of a fermion by the FW Eq. (18) once when calculating $R$ and $T$.

### 4. Boundary Conditions of the Foldy-Wouthuysen equation at a sharp step

Consider a sharp step in potential at $x = a$. I will integrate the FW Eq. (18) through the boundary, normal to the boundary

$$\lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} [\mp \mathcal{L} - (V + E \pm mc^2)] \tilde{\Psi}^{(\pm)} \frac{dx}{c\hbar} = 0. \qquad (31)$$

In this section I introduce the operator $\tilde{\mathcal{L}}$ which commutes with $c\hat{p}_x$ and is defined by

$$\mathcal{L} \tilde{\Psi} = - \tilde{\mathcal{L}} c \hat{p}_x \tilde{\Psi} = - c \hat{p}_x \tilde{\mathcal{L}} \tilde{\Psi}. \qquad (32)$$

I give a discussion of Eq. (32) in Appendix A.

Since

$$\lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} (V + E \pm mc^2) \tilde{\Psi}^{(\pm)} \frac{dx}{c\hbar} = 0, \qquad (33)$$

due to the function integrated being finite but the integration range tending to zero, I have

$$[i \tilde{\mathcal{L}} \tilde{\Psi}]_{a_-}^{a_+} = \lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} \mathcal{L} \tilde{\Psi} \frac{dx}{c\hbar} = 0. \qquad (34)$$

I can understand Eq. (34) when I consider that

$$\mathcal{L} \tilde{\Psi} = i c \hbar \frac{\partial (\tilde{\mathcal{L}} \tilde{\Psi})}{\partial x}. \qquad (35)$$

Eq. (34) shows that $\tilde{\Psi}$ is continuous with respect to the operation of $\tilde{\mathcal{L}}$.

In order to evaluate $\tilde{\mathcal{L}} \tilde{\Psi}$ for the mathematical description of the boundary conditions note

$$\pm \tilde{\mathcal{L}} c \hat{p}_x \tilde{\Psi}^{(\pm)} = (V + E \pm mc^2) \tilde{\Psi}^{(\pm)} = \pm c p_x \tilde{\mathcal{L}} \tilde{\Psi}^{(\pm)}, \qquad (36)$$
where

$$ \mathcal{L} = mc^2 \left[ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 - \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 + \frac{1}{16} \left[ \frac{\hat{p}}{mc} \right]^6 - \dots \right] . \quad (19) $$

Eq. (16) and Eq. (18) are two interpretations of the FW equation.

I demonstrate the unsuitability of my FW approach for calculating magnetic barrier tunneling in Appendix C. I rigorously derive the GFs of the FW Eq. (18) in Appendix D and Appendix F. I uncover and discuss a GF paradox of the FW Eq. (18) in Appendix E.

### 3. Relationship between the transmission of the Dirac equation and its Foldy-Wouthuysen representation

Throughout my numerical and analytic results, Section 14 and Appendix B, I will calculate the reflection coefficient $R$ to deduce the transmission coefficient $T$ as

$$ T = 1 - R. \qquad (20) $$

Calculating $R$, rather than $T$ directly, avoids the need for wavefunction normalization in both the incident and transmitted regions.

Let $\Psi$ be given by Eq. (8) and Eq. (9) as

$$ \Psi = u \exp \left( \frac{ip_x x + ip_y y}{\hbar} \right), \qquad (21) $$

taking $\tilde{\Psi}$ of the form

$$ \tilde{\Psi} = \tilde{u} \exp \left( \frac{ip_x x + ip_y y}{\hbar} \right), \qquad (22) $$

I find by application of Eq. (17) to Eq. (21)

$$ \begin{aligned} \tilde{u}_1 &= u_1 \cos \theta + (d_x - id_y) u_4 \sin \theta, \\ \tilde{u}_2 &= u_2 \cos \theta + (d_x + id_y) u_3 \sin \theta, \\ \tilde{u}_3 &= u_3 \cos \theta - (d_x - id_y) u_2 \sin \theta, \\ \tilde{u}_4 &= u_4 \cos \theta - (d_x + id_y) u_1 \sin \theta, \end{aligned} \qquad (23) $$

where $\hat{d} = \hat{p}/|\hat{p}|$.

Eq. (21), Eq. (22), and Eqs. (23) show $\tilde{\Psi}$ in exponential planar form which is useful for calculating the fermion propagation through homogeneous space.

The incident (i) and reflected (r) waves for planar structures in the Dirac equation description are related to the reflection coefficient in the well-known way:

$$ R = \frac{(\Psi \cdot \Psi^*)^{(r)}}{(\Psi \cdot \Psi^*)^{(i)}}. \qquad (24) $$

For the two linearly independent solutions of the Dirac equation, Eq. (8) and Eq. (9), Eq. (24) gives rise to two degenerate expressions for the reflection coefficient:

$$ R_{\uparrow} = \frac{\left( |\Psi_1|^2 + |\Psi_4|^2 \right)^{(r)}}{\left( |\Psi_1|^2 + |\Psi_4|^2 \right)^{(i)}}, \qquad (25) $$

$$ R_{\downarrow} = \frac{\left( |\Psi_2|^2 + |\Psi_3|^2 \right)^{(r)}}{\left( |\Psi_2|^2 + |\Psi_3|^2 \right)^{(i)}}. \qquad (26) $$

I take note of Eqs. (23) linking $\tilde{u}$ to $u$ and calculate when $d$ is real

$$ |\tilde{u}_1|^2 + |\tilde{u}_4|^2 = |u_1|^2 + |u_4|^2, \qquad (27) $$

$$ |\tilde{u}_2|^2 + |\tilde{u}_3|^2 = |u_2|^2 + |u_3|^2, \qquad (28) $$

therefore $R$ described by the Dirac equation can be rewritten in terms of the FW wavefunctions:

$$ R_{\uparrow} = \frac{\left( |\tilde{\Psi}_1|^2 + |\tilde{\Psi}_4|^2 \right)^{(r)}}{\left( |\tilde{\Psi}_1|^2 + |\tilde{\Psi}_4|^2 \right)^{(i)}}, \qquad (29) $$

$$ R_{\downarrow} = \frac{\left( |\tilde{\Psi}_2|^2 + |\tilde{\Psi}_3|^2 \right)^{(r)}}{\left( |\tilde{\Psi}_2|^2 + |\tilde{\Psi}_3|^2 \right)^{(i)}}. \qquad (30) $$

Since the components of $\tilde{\Psi}$ are degenerate I am only required to evaluate the propagation of a fermion by the FW Eq. (18) once when calculating $R$ and $T$.

### 4. Boundary Conditions of the Foldy-Wouthuysen equation at a sharp step

Consider a sharp step in potential at $x = a$. I will integrate the FW Eq. (18) through the boundary, normal to the boundary

$$ \lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} [\mp \mathcal{L} - (V + E \pm mc^2)] \tilde{\Psi}^{(\pm)} \frac{dx}{c\hbar} = 0. \qquad (31) $$

In this section, I introduce the operator $\hat{\mathcal{L}}$, which commutes with $c\hat{p}_x$ and is defined by

$$ \mathcal{L} \tilde{\Psi} = -\hat{\mathcal{L}} c \hat{p}_x \tilde{\Psi} = -c \hat{p}_x \hat{\mathcal{L}} \tilde{\Psi}. \qquad (32) $$

Appendix A discusses Eq. (32).

Since

$$ \lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} (V + E \pm mc^2) \tilde{\Psi}^{(\pm)} \frac{dx}{c\hbar} = 0, \qquad (33) $$

due to the function integrated being finite but the integration range tending to zero, I have

$$ \left[ i \hat{\mathcal{L}} \tilde{\Psi} \right]_{a_-}^{a_+} = \lim_{(a_+ - a_-) \to 0} \int_{a_-}^{a_+} \hat{\mathcal{L}} \tilde{\Psi} \frac{dx}{c\hbar} = 0. \qquad (34) $$
therefore

$$ \mathcal{L} \tilde{\Psi}^{(\pm)} = \pm \frac{V + E \pm mc^2}{c p_x} \tilde{\Psi}^{(\pm)} = \mp \alpha^{(\pm)} \tilde{\Psi}^{(\pm)}. \quad (37) $$

For my second boundary condition I note that $\tilde{p}\tilde{\Psi}$ always gives a finite observed value, therefore $\tilde{\Psi}$ is continuous everywhere.

My boundary conditions are in stark contrast to the continuity conditions of the Klein-Gordon equation, in that case continuity is with respect to momentum and wavefunction. I verify my boundary conditions numerically and analytically in Section 14 and Appendix B.

## 5. Transfer matrices for a square potential barrier

Inside a layer thickness $D$ the transfer matrix of $\tilde{\Psi}$ is given by

$$ \begin{pmatrix} \exp\left(-\frac{iDp_x}{\hbar}\right) & 0 \\ 0 & \exp\left(+\frac{iDp_x}{\hbar}\right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}, \quad (38) $$

where $\tilde{a}_0$ and $\tilde{a}_1$ denote forward traveling FW wave parts and $\tilde{b}_0$ and $\tilde{b}_1$ denote backward travelling FW wave parts, 0 denotes on the left, 1 denotes on the right. The fermion is refracted at the angle $\phi$ to the layer boundary normal, where

$$ p_x = p \cos \phi. \quad (39) $$

For a boundary at a sharp step, applying the wavefunction continuity condition set out in Section 4 gives

$$ \tilde{a}_0 + \tilde{b}_0 = \tilde{a}_1 + \tilde{b}_1. \quad (40) $$

Secondly, applying wavefunction continuity under the operation of $\mathcal{L}$, derived in Section 4, to Eq. (40) gives

$$ \alpha_0 \tilde{a}_0 - \alpha_0 \tilde{b}_0 = \alpha_1 \tilde{a}_1 - \alpha_1 \tilde{b}_1. \quad (41) $$

Eq. (40) and Eq. (41) can be written in matrix form as

$$ \frac{1}{2} \begin{pmatrix} \left( 1 + \frac{\alpha_1}{\alpha_0} \right) & \left( 1 - \frac{\alpha_1}{\alpha_0} \right) \\ \left( 1 - \frac{\alpha_1}{\alpha_0} \right) & \left( 1 + \frac{\alpha_1}{\alpha_0} \right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}. \quad (42) $$

In the event of tunneling, wave modes become an exponentially decaying and growing set. Transfer matrix multiplication combines these large and small values, this process is numerically unstable and inaccurate. Instead I suggest using the method of Refs. [12, 13] to combine my transfer matrices. The method of Refs. [12, 13] addresses these instability issues by separating the exponentially growing and decaying terms.

## 6. Transfer matrix for a delta potential barrier

Integrating Eq. (18), with $V(x) = g\delta(x)$ and with respect to $dx/c\hbar$, across $[0_-, 0_+]$ gives

$$ \left[ \mp i \mathcal{L} \left( \tilde{a}^{(\pm)} + \tilde{b}^{(\pm)} \right) \right]_{0^{-}}^{0^{+}} = \frac{g}{c\hbar} \left( \tilde{a}^{(\pm)} + \tilde{b}^{(\pm)} \right) \Big|_{x=0} . \quad (43) $$

Continuity of the wavefunction gives

$$ \tilde{a}_0 + \tilde{b}_0 = \tilde{a}_1 + \tilde{b}_1. \quad (44) $$

Eq. (43) and Eq. (44) can be written in the following matrix form

$$ \begin{pmatrix} \left(1 - \frac{ig}{2c\hbar\alpha_k}\right) & \left(-\frac{ig}{2c\hbar\alpha_k}\right) \\ \left(+\frac{ig}{2c\hbar\alpha_k}\right) & \left(1 + \frac{ig}{2c\hbar\alpha_k}\right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}, \quad (45) $$

where I define $\tilde{a}_0$, $\tilde{a}_1$, $\tilde{b}_0$, and $\tilde{b}_1$ in the same way as Section 5.

## 7. 1st order relativistic WKB approximation

Consider the following result from Appendix D derived from the correct fermion interpretation of the FW equation in terms of $\mathcal{L}$,

$$ T(k) = \frac{(2\alpha_k \alpha_q)^2}{\left(\alpha_q^2 - \alpha_k^2\right)^2 \sin^2\left(\frac{2m^*qa}{\hbar}\right) + (2\alpha_k \alpha_q)^2} \quad (46) $$

for a rectangular barrier given by

$$ V(x) = \begin{cases} V & \text{for } |x| < a, \\ 0 & \text{for } |x| > a. \end{cases} \quad (47) $$

Eq. (46) can be approximated as

$$ T(k) \approx \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)^2 \exp \left( -\frac{4im^*qa}{\hbar} \right) \quad (48) $$

by assuming exponential wavefunction decay within the barrier. Taking logarithms of Eq. (48), the non-logarithmic term dominates the logarithmic term so that

$$ \log_e T(k) \approx -\frac{4im^*qa}{\hbar}. \quad (49) $$

Assuming independence of transmission events, the total transmission probability for crossing an overall barrier is

$$ T \approx \prod_{i=1}^{\infty} T_i, \quad (50) $$

$$ \log_e \prod_{i=1}^{\infty} T_i \approx \sum_{i=1}^{\infty} \frac{2im_i^* q_i \Delta x}{\hbar} = -\frac{2i}{\hbar} \int p_x dx, \quad (51) $$
I can understand Eq. (34) when I consider that

$$ \mathcal{L}\tilde{\Psi} = ic\hbar \frac{\partial(\mathcal{L}\tilde{\Psi})}{\partial x}. \qquad (35) $$

Eq. (34) shows that $\tilde{\Psi}$ is continuous under the operation of $\tilde{\mathcal{L}}$.

In order to evaluate $\tilde{\mathcal{L}}\tilde{\Psi}$ for the mathematical description of the boundary conditions note

$$ \pm \mathcal{L}c\hat{p}_x \tilde{\Psi}^{(\pm)} = (V + E \pm mc^2) \tilde{\Psi}^{(\pm)} = \pm c p_x \mathcal{L} \tilde{\Psi}^{(\pm)}, \quad (36) $$

therefore

$$ \tilde{\mathcal{L}}\tilde{\Psi}^{(\pm)} = \pm \frac{V + E \pm mc^2}{c p_x} \tilde{\Psi}^{(\pm)} = \mp \alpha^{(\pm)} \tilde{\Psi}^{(\pm)}. \quad (37) $$

For my second boundary condition I note that $\hat{p}\tilde{\Psi}$ always gives a finite observed value, therefore $\tilde{\Psi}$ is continuous everywhere.

My boundary conditions are in stark contrast to the continuity conditions of the Klein-Gordon equation; in that case, continuity is with respect to momentum and wavefunction. I verify my boundary conditions numerically and analytically in Section 14 and Appendix B.

## 5. Transfer matrices for a rectangular potential barrier

Inside a layer thickness $D$ the transfer matrix of $\tilde{\Psi}$ is given by

$$ \begin{pmatrix} \exp\left(-\frac{iDp_x}{\hbar}\right) & 0 \\ 0 & \exp\left(+\frac{iDp_x}{\hbar}\right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}, \quad (38) $$

where $\tilde{a}_0$ and $\tilde{a}_1$ denote forward traveling FW wave parts and $\tilde{b}_0$ and $\tilde{b}_1$ denote backward traveling FW wave parts, 0 denotes on the left, 1 denotes on the right. The fermion is refracted at an angle $\phi$ to the layer boundary normal, where

$$ p_x = p \cos \phi. \qquad (39) $$

For a boundary at a sharp step, applying the wavefunction continuity condition set out in Section 4 gives

$$ \tilde{a}_0 + \tilde{b}_0 = \tilde{a}_1 + \tilde{b}_1. \qquad (40) $$

Secondly, applying wavefunction continuity under the operation of $\tilde{\mathcal{L}}$, derived in Section 4, to Eq. (40) gives

$$ \alpha_0 \tilde{a}_0 - \alpha_0 \tilde{b}_0 = \alpha_1 \tilde{a}_1 - \alpha_1 \tilde{b}_1. \qquad (41) $$

Eq. (40) and Eq. (41) can be written in matrix form as

$$ \frac{1}{2} \begin{pmatrix} \left( 1 + \frac{\alpha_1}{\alpha_0} \right) & \left( 1 - \frac{\alpha_1}{\alpha_0} \right) \\ \left( 1 - \frac{\alpha_1}{\alpha_0} \right) & \left( 1 + \frac{\alpha_1}{\alpha_0} \right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}. \quad (42) $$

During tunneling, wave modes either decay or grow exponentially. Transfer matrix multiplication combines these large and small values, this process is numerically unstable and inaccurate. Instead, I suggest using the method of Refs. [14, 15] to combine my transfer matrices. The method of Refs. [14, 15] addresses these instability issues by separating the exponentially growing and decaying terms.

## 6. Transfer matrix for a delta potential barrier

Integrating Eq. (18), with $V(x) = g\delta(x)$ and with respect to $dx/c\hbar$, across $[0_-, 0_+]$ gives

$$ [\mp i \tilde{\mathcal{L}} (\tilde{a}^{(\pm)} + \tilde{b}^{(\pm)})]_{0^{-}}^{0_{+}} = \frac{g}{c\hbar} (\tilde{a}^{(\pm)} + \tilde{b}^{(\pm)}) \Big|_{x=0}. \quad (43) $$

Continuity of the wavefunction gives

$$ \tilde{a}_0 + \tilde{b}_0 = \tilde{a}_1 + \tilde{b}_1. \qquad (44) $$

Eq. (43) and Eq. (44) can be written in the following matrix form

$$ \begin{pmatrix} \left(1 - \frac{ig}{2c\hbar\alpha_k}\right) & \left(-\frac{ig}{2c\hbar\alpha_k}\right) \\ \left(+\frac{ig}{2c\hbar\alpha_k}\right) & \left(1 + \frac{ig}{2c\hbar\alpha_k}\right) \end{pmatrix} \begin{pmatrix} \tilde{a}_1 \\ \tilde{b}_1 \end{pmatrix} = \begin{pmatrix} \tilde{a}_0 \\ \tilde{b}_0 \end{pmatrix}, \quad (45) $$

where I define $\tilde{a}_0$, $\tilde{a}_1$, $\tilde{b}_0$, and $\tilde{b}_1$ in the same way as Section 5.

## 7. 1st order relativistic WKB approximation

Consider the following result from Appendix D derived from the correct fermion interpretation of the FW equation in terms of $\mathcal{L}$,

$$ T(k) = \frac{(2\alpha_k \alpha_q)^2}{(\alpha_q^2 - \alpha_k^2)^2 \sin^2\left(\frac{2m^*qa}{\hbar}\right) + (2\alpha_k \alpha_q)^2} \quad (46) $$

for a rectangular barrier given by

$$ V(x) = \begin{cases} V & \text{for } |x| < a, \\ 0 & \text{for } |x| > a. \end{cases} \qquad (47) $$

Eq. (46) can be approximated as

$$ T(k) \approx \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)^2 \exp\left(-\frac{4im^*qa}{\hbar}\right) \quad (48) $$

by assuming exponential wavefunction decay within the barrier. Taking logarithms of Eq. (48), the non-logarithmic term dominates the logarithmic term so that

$$ \log_e T(k) \approx -\frac{4im^*qa}{\hbar}. \qquad (49) $$
where the $T_i$ are the transmission probabilities of the individual barriers. I see from Eq. (51) that inside a smoothly varying layer the FW wavefunction is approximately given by

$$ \tilde{\Psi} \approx \exp \left( \pm \frac{i}{\hbar} \int p_x dx \right) \prod_{i=1}^{N} \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)_i . \quad (52) $$

## 8. Iterative relativistic WKB approximation

Consider a wavefunction of the form

$$ \tilde{\Psi} = \exp \left( \frac{i\Phi}{\hbar} \right), \quad (53) $$

then the boson FW interpretation Eq. (16) becomes

$$ \left[ \sqrt{c^2 \left( \frac{\partial \Phi}{\partial x} \right)^2 - c^2 i \hbar \frac{\partial^2 \Phi}{\partial x^2} + m^2 c^4 - (V + E)} \right] \tilde{\Psi} = 0. \quad (54) $$

Eq. (54) can be exactly rewritten as

$$ c^2 \left( \frac{\partial \Phi}{\partial x} \right)^2 - c^2 i \hbar \frac{\partial^2 \Phi}{\partial x^2} + m^2 c^4 = (V + E)^2 . \quad (55) $$

I now make use of $(E + V)^2 = m^2 c^4 + p^2 c^2$ in Eq. (55) and find that as expected

$$ \hat{p}_x^2 \tilde{\Psi} = p_x^2 \tilde{\Psi}, \quad \frac{\partial \Phi}{\partial x} = \sqrt{p_x^2 + i\hbar \frac{\partial^2 \Phi}{\partial x^2}}. \quad (56) $$

I know from Section 7

$$ \pm \frac{\partial \Phi}{\partial x} \approx p_x, \quad (57) $$

suggesting a next approximation on from Eq. (52) is

$$ \Phi \approx \pm \int \sqrt{p_x^2 \pm i \hbar \frac{\partial p_x}{\partial x}} dx. \quad (58) $$

Eq. (58) follows from the integration of Eq. (56). This procedure can be iterated to an $n^{th}$ order WKB approximation in the boson FW interpretation.

## 9. 2nd order relativistic WKB approximation

Expanding $\Phi(x)$ in powers of $\hbar$

$$ \Phi(x) = \Phi_0(x) + \hbar\Phi_1(x) + \hbar^2\Phi_2(x) + \dots, \quad (59) $$

then substituting Eq. (59) into Eq. (56) gives

$$ (\Phi'_0 + \hbar\Phi'_1 + \dots)^2 - i\hbar(\Phi''_0 + \hbar\Phi''_1 + \dots) = p_x^2. \quad (60) $$

Equating terms in powers of $\hbar$, in Eq. (60), I see

$$ \pm \Phi'_0 = p_x, \quad 2\Phi'_0\Phi'_1 - i\Phi''_0 = 0, \quad (61) $$

so that

$$ \Phi_0 = \pm \int p_x dx, \quad \Phi_1 = \frac{i}{2} \log_e \Phi'_0. \quad (62) $$

Eqs. (62), Eq. (59), and Eq. (53) give

$$ \tilde{\Psi} \approx \frac{\exp \left( \pm \frac{i}{\hbar} \int p_x dx \right)}{\sqrt{p_x}}. \quad (63) $$

This result is identical to the Klein-Gordon WKB approximation because on the way I approximated FW Eq. (18) with Eq. (16).

## 10. Discussion of fermion and boson WKB approximations

When I made my WKB approximation in Section 7, from what I will show in Section 14 and Appendix B to be the correct form Eq. (18) of the FW equation, I saw in the non-exponential factor of Eq. (48) appearances of $p_i$ changed to appearances of $\alpha_i$ going from boson to fermion descriptions. I have noticed, that for tunneling through slowly varying potentials, Eq. (48), Eq. (50) and Eq. (63) give the implication:

$$ \prod_{i=1}^{N} \left( \frac{2kq}{k^2 - q^2} \right)_i^2 \approx \frac{p_1}{p_N} \Rightarrow \prod_{i=1}^{N} \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)_i^2 \approx \frac{\alpha_{p_1}^{(\pm)}}{\alpha_{p_N}^{(\pm)}}. \quad (64) $$

Eq. (64) suggests that the WKB approximation for fermions, in the correct FW representation, will take the form

$$ \tilde{\Psi} \approx \frac{\exp \left( \pm \frac{i}{\hbar} \int p_x dx \right)}{\sqrt{\alpha_{p_x}}} . \quad (65) $$

Going from Eq. (63) to Eq. (65) I swapped an appearance of $p_i$ to $\pm \alpha_i^{(\pm)}$, from boson to fermion descriptions, as suggested by Eq. (48), Eq. (50) and Eq. (63).

## 11. Connection formulae

It is suggested by Section 8 that for the WKB approximation to be appropriate

$$ |p_x|^2 \gg \left| \hbar \frac{\partial p_x}{\partial x} \right|. \quad (66) $$

Eq. (66) suggests that when $p_x \approx 0$ the WKB method fails. I need a connection formula between the regions where $p_x^2 \neq 0$.

I assume that when Eq. (66) fails

$$ mc \gg |p_x|, \quad (67) $$

so that I may treat the connecting region with the $c = \infty$ non-relativistic limit of the FW Eq. (18). The wavefunction
Assuming independence of transmission events, the total transmission probability for crossing an overall barrier is

$$T \approx \prod_{i=1}^{\infty} T_i, \qquad (50)$$

$$\log_e \prod_{i=1}^{\infty} T_i \approx \sum_{i=1}^{\infty} -\frac{2im_i^* q_i \Delta x}{\hbar} = -\frac{2i}{\hbar} \int p_x dx, \quad (51)$$

where the $T_i$ are the transmission probabilities of the individual barriers. I see from Eq. (51) that inside a smoothly varying layer the FW wavefunction is approximately given by

$$\tilde{\Psi} \approx \exp \left( \pm \frac{i}{\hbar} \int p_x dx \right) \prod_{i=1}^{N} \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)_i. \quad (52)$$

## 8. Iterative relativistic WKB approximation

Consider a wavefunction of the form

$$\tilde{\Psi} = \exp \left( \frac{i\Phi}{\hbar} \right), \qquad (53)$$

then the boson FW interpretation Eq. (16) becomes

$$\left[ \sqrt{c^2 \left( \frac{\partial \Phi}{\partial x} \right)^2 - c^2 i \hbar \frac{\partial^2 \Phi}{\partial x^2} + m^2 c^4 - (V + E)} \right] \tilde{\Psi} = 0. \qquad (54)$$

Eq. (54) can be exactly rewritten as

$$c^2 \left( \frac{\partial \Phi}{\partial x} \right)^2 - c^2 i \hbar \frac{\partial^2 \Phi}{\partial x^2} + m^2 c^4 = (V + E)^2. \quad (55)$$

I now make use of $(E+V)^2 = m^2c^4 + p^2c^2$ in Eq. (55) and find that as expected

$$\hat{p}_x^2 \tilde{\Psi} = p_x^2 \tilde{\Psi}, \quad \frac{\partial \Phi}{\partial x} = \sqrt{p_x^2 + i\hbar \frac{\partial^2 \Phi}{\partial x^2}}. \qquad (56)$$

I know from Section 7

$$\pm \frac{\partial \Phi}{\partial x} \approx p_x, \qquad (57)$$

suggesting a next approximation on from Eq. (52) is

$$\Phi \approx \pm \int \sqrt{p_x^2 \pm i\hbar \frac{\partial p_x}{\partial x}} dx. \qquad (58)$$

Eq. (58) follows from the integration of Eq. (56). This procedure can be iterated to an $n^{th}$ order WKB approximation in the boson FW interpretation.

## 9. 2nd order relativistic WKB approximation

Expanding $\Phi(x)$ in powers of $\hbar$

$$\Phi(x) = \Phi_0(x) + \hbar\Phi_1(x) + \hbar^2\Phi_2(x) + \dots, \qquad (59)$$

then substituting Eq. (59) into Eq. (56) gives

$$(\Phi'_{0} + \hbar\Phi'_{1} + \dots)^2 - i\hbar(\Phi''_{0} + \hbar\Phi''_{1} + \dots) = p_x^2. \quad (60)$$

Equating terms in powers of $\hbar$, in Eq. (60), I see

$$\pm \Phi'_{0} = p_x, \quad 2\Phi'_{0}\Phi'_{1} - i\Phi''_{0} = 0, \qquad (61)$$

so that

$$\Phi_0 = \pm \int p_x dx, \quad \Phi_1 = \frac{i}{2} \log_e \Phi'_0. \qquad (62)$$

Eqs. (62), Eq. (59), and Eq. (53) give

$$\tilde{\Psi} \approx \frac{\exp \left( \pm \frac{i}{\hbar} \int p_x dx \right)}{\sqrt{p_x}}. \qquad (63)$$

This result is identical to the Klein-Gordon WKB approximation because on the way I approximated FW Eq. (18) with Eq. (16).

## 10. Discussion of fermion and boson WKB approximations

When I made my WKB approximation in Section 7, from what I will show in Section 14 and Appendix B to be the correct form Eq. (18) of the FW equation, I saw in the non-exponential factor of Eq. (48) appearances of $p_i$ changed to appearances of $\alpha_i$ going from boson to fermion descriptions. I have noticed, that for tunneling through slowly varying potentials, Eq. (48), Eq. (50) and Eq. (63) give the implication:

$$\prod_{i=1}^{N} \left( \frac{2kq}{k^2 - q^2} \right)_i^2 \approx \frac{p_1}{p_N} \Rightarrow \prod_{i=1}^{N} \left( \frac{2\alpha_k \alpha_q}{\alpha_k^2 - \alpha_q^2} \right)_i^2 \approx \frac{\alpha_{p_1}^{(\pm)}}{\alpha_{p_N}^{(\pm)}}. \qquad (64)$$

Eq. (64) suggests that the WKB approximation for fermions, in the correct FW representation, will take the form

$$\tilde{\Psi} \approx \frac{\exp \left( \pm \frac{i}{\hbar} \int p_x dx \right)}{\sqrt{\alpha_{p_x}}} \qquad (65)$$

Going from Eq. (63) to Eq. (65) I swapped an appearance of $p_i$ to $\pm \alpha_i^{(\pm)}$, from boson to fermion descriptions, as suggested by Eq. (48), Eq. (50) and Eq. (63).

## 11. Connection formulae

It is suggested by Section 8 that for the WKB approximation to be appropriate

$$|p_x|^2 \gg |\hbar \frac{\partial p_x}{\partial x}|. \qquad (66)$$
in the connecting region is then given by the Airy function solutions to the Schrödinger equation. This approach leads to the two connected solutions for $\tilde{\Psi}$ between the classical and non-classical regions:

**Classical region $p^2 > 0$ and $x < x_1$**

$$A_1 \frac{\exp\left(+\frac{i}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}} + B_1 \frac{\exp\left(-\frac{i}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}}, \quad (68)$$

**Non-classical region $p^2 < 0$ and $x > x_1$**

$$A_2 \frac{\exp\left(-\frac{1}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}} + B_2 \frac{\exp\left(+\frac{1}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}}, \quad (69)$$

where

$$2A_2 = -i (A_1 + iB_1) \exp(+i\pi/4) \quad (70)$$

and

$$B_2 = +i (A_1 - iB_1) \exp(-i\pi/4). \quad (71)$$

In Eq. (68) and Eq. (69) I have used the asymptotic forms of the Airy functions and taken $p = 0$ to be at $x_1$.

The wavefunctions Eq. (68) and Eq. (69) are in the WKB form derived for relativistic fermions in Section 7 to Section 10, therefore the solutions in both regions can be extended to relativistic energies.

## 12. Extention to periodic structures

Consider the periodic potential

$$V(x) = \sum_{G} V_{G} \exp \left( + \frac{iGx}{\hbar} \right). \quad (72)$$

Eq. (72) is a Fourier expansion of the potential, where $G = 2\pi n/a$ and $n$ is positive and negative integers. The wave-function follows from Bloch's theorem for periodic crystals

$$\tilde{\Psi}^{(\pm)} = \exp\left(+\frac{ikx}{\hbar}\right) \sum_{g} A_g \exp\left(-\frac{igx}{\hbar}\right), \quad (73)$$

where $g = 2\pi m/a$ and $m$ is positive and negative integers. Substituting Eq. (72) and Eq. (73) into

$$[\mp \mathcal{L} - (E \pm mc^2 + V)] \tilde{\Psi}^{(\pm)} = 0 \quad (74)$$

and multiplying through by $\exp(i(G' - k)x/\hbar)$, where $G' = 2\pi n'/a$, then integrating $x$ over one period $a$, I find that the first, second, and third terms are only non-zero for $g = G'$ while the fourth term is only non-zero for $g - G = G'$. Therefore, with rearrangement

$$\sum_{g} (\delta_{gG'} W_{g} - V_{g-G'}) A_{g} = (E \pm mc^{2}) A_{G'}, \quad (75)$$

where

$$W_g = - (\alpha_k^{(\pm)} - \alpha_g^{(\pm)}) (ck - cg). \quad (76)$$

The eigenvalue problem described by Eq. (75) can be solved for varying $k$ to give the eigenvalues $(E \pm mc^2)$. The relationship between $k$ and $E$ defines the band structure.

For planar structures of period $D$ I have by Bloch's theorem

$$\tilde{\Psi}(x + D) = \exp\left(\frac{iKD}{\hbar}\right) \tilde{\Psi}(x), \quad (77)$$

which in terms of transfer matrices can be written as

$$M \tilde{\Psi} = \exp\left(\frac{iKD}{\hbar}\right) \tilde{\Psi} = \lambda \tilde{\Psi} \implies |M - \lambda I| = 0, \quad (78)$$

where

$$M = \begin{pmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{pmatrix}. \quad (79)$$

For fermions to propagate through the periodic structure $K$ must be real, therefore by Eqs. (78)

$$-1 \le \cos\left(\frac{KD}{\hbar}\right) = \frac{M_{11} + M_{22}}{2} \le +1. \quad (80)$$

Band-gaps arise in 1D periodic potentials if the condition in Eq. (80) is not met.

# 13. Calculation of transfer matrix parameters

## 13.1. Calculation of the angle of refraction

To match the Dirac spinors of Eq. (8) and Eq. (9) either side of a potential step at $x = 0$ it must be that

$$p_y^{(i)}|_{x=0} = p_y^{(t)}|_{x=0}, \quad (81)$$

therefore, since

$$\hat{p} = \hat{x}p \cos \phi + \hat{y}p \sin \phi, \quad (82)$$

I arrive at Snell's law for the relativistic fermion

$$\frac{\sin \phi^{(i)}}{\sin \phi^{(t)}} = \frac{p^{(t)}}{p^{(i)}}. \quad (83)$$

I do not analytically verify non-normal incident relativistic tunneling and transmission in this article.

## 13.2. Tunneling between layers of varying electron mass

In free space by special relativity

$$E^2 = m^2 c^4 \left(1 + \left(\frac{k}{c}\right)^2\right), \quad (84)$$

inside a layer by special relativity

$$(E + V)^2 = m^{*2} c^{*4} \left(1 + \left(\frac{q}{c^*}\right)^2\right), \quad (85)$$
Eq. (66) suggests that when $p_x \approx 0$ the WKB method fails. I need a connection formula between the regions where $p_x^2 \neq 0$.

I assume that when Eq. (66) fails

$$mc \gg |p_x|, \qquad (67)$$

so that I may treat the connecting region with the $c = \infty$ non-relativistic limit of the FW Eq. (18). The wavefunction in the connecting region is then given by the Airy function solutions to the Schrödinger equation. This approach leads to the two connected solutions for $\tilde{\Psi}$ between the classical and non-classical regions:

**Classical region $p^2 > 0$ and $x < x_1$**

$$A_1 \frac{\exp\left(+\frac{i}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}} + B_1 \frac{\exp\left(-\frac{i}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}}, \quad (68)$$

**Non-classical region $p^2 < 0$ and $x > x_1$**

$$A_2 \frac{\exp\left(-\frac{1}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}} + B_2 \frac{\exp\left(+\frac{1}{\hbar} \int_{x_1}^{x} |p| \, dx'\right)}{\sqrt{|p|}}, \quad (69)$$

where

$$2A_2 = -i (A_1 + iB_1) \exp(+i\pi/4) \quad (70)$$

and

$$B_2 = +i (A_1 - iB_1) \exp(-i\pi/4). \quad (71)$$

In Eq. (68) and Eq. (69) I have used the asymptotic forms of the Airy functions and taken $p = 0$ to be at $x_1$.

The wavefunctions Eq. (68) and Eq. (69) are in the WKB form derived for relativistic fermions in Section 7 to Section 10, therefore the solutions in both regions can be extended to relativistic energies.

## 12. Extension to periodic structures

Consider the periodic potential

$$V(x) = \sum_{G} V_{G} \exp \left( + \frac{iGx}{\hbar} \right). \quad (72)$$

Eq. (72) is a Fourier expansion of the potential, where $G = 2\pi n/a$ and $n$ is positive and negative integers. The wavefunction follows from Bloch's theorem for periodic crystals

$$\tilde{\Psi}^{(\pm)} = \exp\left(+\frac{ikx}{\hbar}\right) \sum_{g} A_g \exp\left(-\frac{igx}{\hbar}\right), \quad (73)$$

where $g = 2\pi m/a$ and $m$ is positive and negative integers. Substituting Eq. (72) and Eq. (73) into

$$[\mp\mathcal{L} - (E \pm mc^2 + V)] \tilde{\Psi}^{(\pm)} = 0 \quad (74)$$

and multiplying through by $\exp(i(G' - k)x/\hbar)$, where $G' = 2\pi n'/a$, then integrating $x$ over one period $a$, I find that the first, second, and third terms are only non-zero for $g = G'$ while the fourth term is only non-zero for $g - G = G'$. Therefore, with rearrangement

$$\sum_{g} (\delta_{gG'} W_{g} - V_{g-G'}) A_{g} = (E \pm mc^{2}) A_{G'}, \quad (75)$$

where

$$W_g = -c(k-g) \alpha_{k-g}^{(\pm)}. \quad (76)$$

The eigenvalue problem described by Eq. (75) can be solved for varying $k$ to give the eigenvalues $(E \pm mc^2)$. The relationship between $k$ and $E$ defines the band structure.

For planar structures of period $D$ I have by Bloch's theorem

$$\tilde{\Psi}(x+D) = \exp\left(\frac{iKD}{\hbar}\right) \tilde{\Psi}(x), \quad (77)$$

which in terms of transfer matrices can be written as

$$M \tilde{\Psi} = \exp\left(\frac{iKD}{\hbar}\right) \tilde{\Psi} = \lambda \tilde{\Psi} \implies |M - \lambda I| = 0, \quad (78)$$

where

$$M = \begin{pmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{pmatrix}. \quad (79)$$

For fermions to propagate through the periodic structure, $K$ must be real. Consequently, by Eqs. (78)

$$-1 \le \cos\left(\frac{KD}{\hbar}\right) = \frac{M_{11} + M_{22}}{2} \le +1. \quad (80)$$

Band gaps arise in 1D periodic potentials if the condition in Eq. (80) is not satisfied.

# 13. Calculation of transfer matrix parameters

## 13.1. Calculation of the angle of refraction

To match the Dirac spinors of Eq. (8) and Eq. (9) either side of a potential step at $x = 0$ it must be that

$$p_y^{(i)}|_{x=0} = p_y^{(t)}|_{x=0}, \quad (81)$$

therefore, since

$$\hat{p} = \hat{x}p \cos \phi + \hat{y}p \sin \phi, \quad (82)$$

I arrive at Snell's law for the relativistic fermion

$$\frac{\sin \phi^{(i)}}{\sin \phi^{(t)}} = \frac{p^{(t)}}{p^{(i)}}. \quad (83)$$

I do not analytically verify non-normal incident relativistic tunneling and transmission in this article.
where $m^*$ is the effective mass of the electron in the layer and $c^*$ is the corresponding effective speed of light. Comparing Eq. (84) and Eq. (85), I find

$$ \frac{q^2}{k^2} = \left( \frac{mc}{m^*c^*} \right)^2 + \frac{mc^2}{m^*c^{*2}} \frac{2g(k/c)V}{m^*k^2} + \left( \frac{V}{m^*c^*k} \right)^2 + \frac{(mc^2 - m^*c^{*2})(mc^2 + m^*c^{*2})}{(m^*c^*)^2} \frac{1}{k^2}. \quad (86) $$

In the last term of Eq. (86) I used the difference of two squares to improve numerical accuracy. When $|k| < c$ the function $g(k/c)$ is given by the Taylor expansion

$$ g(k/c) = \sqrt{1 + \left(\frac{k}{c}\right)^2} = 1 + \frac{(k/c)^2}{2} - \frac{(k/c)^4}{8} + \dots \quad (87) $$

and when $|k| > c$ the function $g(k/c)$ is given by the Taylor expansion

$$ g(k/c) = \left(\frac{k}{c}\right) \left( 1 + \frac{(c/k)^2}{2} - \frac{(c/k)^4}{8} + \dots \right). \quad (88) $$

Computationally I add smallest terms together first for improved numerical accuracy. In terms of $g(k/c)$ and the x-component of electron momentum $m^*q_x$ I calculate the parameter $\alpha_q^{(\pm)}$ to be

$$ \alpha_q^{(\pm)} = - \frac{mc^2 g(k/c) + V \pm m^*c^{*2}}{c^*m^*q_x}. \quad (89) $$

### 13.3. Tunneling into a layer where electron mass vanishes

In free space the fermion momentum dispersion is given by Eq. (84) whereas inside the layer electron mass has vanished and therefore

$$ E + V = qc^*, \quad (90) $$

where $q$ is the fermion momentum in the layer, and $c^*$ is the speed of light in the layer. Comparing Eq. (84) and Eq. (90), I find

$$ qc^* = mc^2g(k/c) + V \quad (91) $$

and

$$ \alpha_q = -\frac{q}{q_x}. \quad (92) $$

## 14. Results

This section is restricted to quantum tunneling. For the structures I investigate in this section I introduce a length scale $a$ given in terms of electron mass $m$, the speed of light $c$, and $\hbar$

$$ a = 5 \frac{\hbar}{mc}. \quad (93) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1:** (Color online) (a) Transmission through the potential barrier Eq. (94) as a function of incident momentum $p$, for fermions (red) and bosons (orange). (b) Convergence of the fermion transmission with increasing number of relativistic correction Taylor expansion terms 0, 5, 10, 20, 40, 80 as labelled.</figcaption></figure>

In Appendix B I investigate reflection of electrons from a sharp step in potential. In Appendix B I find exact analytic agreement for $R$ between the FW and Dirac spinor representations.

The striking feature of all my results is that I see exact analytic agreement between $T$ and $R$ calculated in the FW and Dirac spinor representations for all momentum $p$, even when $|p| > mc$. This is a significant finding since when $|p| > mc$ the Taylor expansion of $\mathcal{L}$ is a divergent series. I think that conservation of energy is restricting my evaluation of $\bar{\mathcal{L}}$ on $\tilde{\Psi}^{(\pm)}$ to the particular values required to meet the physical boundary conditions for the Dirac Hamiltonian.

I model a resonant-tunneling diode to compare my WKB approximations of Section 9 and Section 10 as the final example in this section.
## 13.2. Tunneling between layers of varying electron mass

In free space by special relativity

$$E^2 = m^2 c^4 \left( 1 + \left( \frac{k}{c} \right)^2 \right), \qquad (84)$$

inside a layer by special relativity

$$(E + V)^2 = m^{*2} c^{*4} \left( 1 + \left( \frac{q}{c^*} \right)^2 \right), \qquad (85)$$

where $m^*$ is the effective mass of the electron in the layer and $c^*$ is the corresponding effective speed of light. Comparing Eq. (84) and Eq. (85), I find

$$\frac{q^2}{k^2} = \left(\frac{mc}{m^*c^*}\right)^2 + \frac{mc^2}{m^*c^{*2}} \frac{2g(k/c)V}{m^*k^2} + \left(\frac{V}{m^*c^*k}\right)^2 + \frac{\left(m^2c^2 - m^*c^{*2}\right)\left(m^2c^2 + m^*c^{*2}\right)}{(m^*c^*)^2} \frac{1}{k^2}. \quad (86)$$

In the last term of Eq. (86) I used the difference of two squares to improve numerical accuracy. When $|k| < c$ the function $g(k/c)$ is given by the Taylor expansion

$$g(k/c) = \sqrt{1 + \left(\frac{k}{c}\right)^2} = 1 + \frac{(k/c)^2}{2} - \frac{(k/c)^4}{8} + \dots \quad (87)$$

and when $|k| > c$ the function $g(k/c)$ is given by the Taylor expansion

$$g(k/c) = \left(\frac{k}{c}\right) \left( 1 + \frac{(c/k)^2}{2} - \frac{(c/k)^4}{8} + \dots \right). \quad (88)$$

Computationally, I add the smallest terms together first to improve numerical accuracy. In terms of $g(k/c)$ and the x-component of electron momentum $m^*q_x$ I calculate the parameter $\alpha_q^{(\pm)}$ to be

$$\alpha_q^{(\pm)} = - \frac{mc^2 g(k/c) + V \pm m^* c^{*2}}{c^* m^* q_x}. \qquad (89)$$

## 13.3. Tunneling into a layer where electron mass vanishes

In free space the fermion momentum dispersion is given by Eq. (84) whereas inside the layer electron mass has vanished and therefore

$$E + V = qc^*, \qquad (90)$$

where $q$ is the fermion momentum in the layer, and $c^*$ is the speed of light in the layer. Comparing Eq. (84) and Eq. (90), I find

$$qc^* = mc^2g(k/c) + V \qquad (91)$$

and

$$\alpha_q = -\frac{q}{q_x}. \qquad (92)$$

<figure><img src="image_4.png" /><figcaption>Figure 1: (Color online) (a) Transmission through the potential barrier Eq. (94) as a function of incident momentum $p$, for fermions (red) and bosons (orange). (b) Convergence of the fermion transmission with increasing number of relativistic correction Taylor expansion terms 0, 5, 10, 20, 40, 80 as labelled.</figcaption></figure>

## 14. Results

This section is restricted to quantum tunneling. For the structures I investigate in this section I introduce a length scale $a$ given in terms of electron mass $m$, the speed of light $c$, and $\hbar$

$$a = 5 \frac{\hbar}{mc}. \qquad (93)$$

In Appendix B I investigate reflection of electrons from a sharp step in potential. In Appendix B I find exact analytic agreement for $R$ between the FW and Dirac spinor representations.

The striking feature of all my results is that I see exact analytic agreement between $T$ and $R$ calculated in the FW and Dirac spinor representations for all momentum $p$, even
<figure><img src="image_4.png" /><figcaption>**Figure 2:** (Color online) Repeat of Fig. 1 except that now electron mass vanishes inside the barrier.</figcaption></figure>

For my numerical computations I use Python 3.8.

## 14.1. Tunneling through a barrier

The system under consideration is a potential barrier described by

$$V(x) = \begin{cases} 0 & \text{for } |x| > a, \\ 3mc^2 & \text{for } |x| < a. \end{cases} \qquad (94)$$

Since $V$ is negative, the momentum of the electron inside the barrier becomes imaginary when the following condition is met:

$$-mc^2 < E + V < +mc^2, \qquad (95)$$

then the electron is quantum tunneling.

The analytic transmission for Eq. (94) is derived in Appendix D as:

$$T(k) = \frac{4\kappa^2}{4\kappa^2 + (1 - \kappa^2)^2 \sin^2\left(\frac{2mqa}{\hbar}\right)}, \qquad (96)$$

where I take $\kappa$ to be $\kappa_f$ in the FW representation, and $\kappa_d$ in the Dirac spinor representation. For the derivation of Appendix D in the FW representation

$$\kappa_f^{(\pm)} = \frac{\alpha_k^{(\pm)}}{\alpha_q^{(\pm)}} = \frac{q_x}{k_x} \frac{E \pm mc^2}{E \pm mc^2 + V}. \qquad (97)$$

It was shown by the authors of Refs. [14, 15] by deriving $T$ in the Dirac spinors representation of quantum tunneling

$$\kappa_d^{(\pm)} = \frac{q_x}{k_x} \frac{E \pm mc^2}{E \pm mc^2 + V}, \qquad (98)$$

which is in exact analytic agreement with $\kappa_f^{(\pm)}$, as expected.

The numerical results demonstrating relativistic tunneling of electrons are shown in Fig. 1. For comparison I also show the transmission of tunneling bosons with identical mass to the electron.

Fig. 1(b) shows the absolute error calculated as the difference between $T$ evaluated using the Taylor expansions $g(k/c)$ in $\kappa$ or the Python 3.8 `cmath` module square root function.

The region of slow convergence around $|k/c| = 1$ in Fig. 1(b) is an artifact of the point in $g(k/c)$ about which I have taken my Taylor expansions.

## 14.2. Tunneling through a barrier with vanishing electron mass inside the barrier

In this subsection I almost repeat the demonstration of Section 14.1 except that now inside the barrier the electron mass is strictly vanishing.

Since $m = 0$ when $|x| < a$ I must re-evaluate $\kappa_d$ for Eq. (96):

$$\kappa_d^{(\pm)} = \frac{q_x}{p_x} \frac{E \pm mc^2}{E + V}, \qquad (99)$$

which is in exact analytic agreement with $\kappa_f^{(\pm)}$, as expected.

The numerical results demonstrating tunneling of relativistic electrons are shown in Fig. 2. For comparison I also show the transmission of tunneling bosons with identical mass to the electron.

In Fig. 2(b) the absolute error is calculated as the difference between $T$ evaluated using the Taylor expansions $g(k/c)$ in $\kappa$ or the Python 3.8 `cmath` module square root function.

The region of slow convergence around $|k/c| = 1$ in Fig. 2(b) is an artifact of the point in $g(k/c)$ about which I have taken my Taylor expansions.
<figure><img src="image_5.png" /><figcaption>**Figure 2:** (Color online) Repeat of Fig. 1, except that now electron mass vanishes inside the barrier.</figcaption></figure>

when $|p| > mc$. This finding is significant, as the Taylor expansion of $\mathcal{L}$ diverges when $|p| > mc$. I think that conservation of energy is restricting my evaluation of $\tilde{\Psi}(\pm)$ to the particular values required to meet the physical boundary conditions for the Dirac Hamiltonian.

I model a resonant-tunneling diode to compare my WKB approximations of Section 9 and Section 10 as the final example in this section.

For my numerical computations I use Python 3.8.

## 14.1. Tunneling through a barrier

The system under consideration is a potential barrier described by

$$V(x) = - \begin{cases} 0 & \text{for } |x| > a, \\ 3mc^2 & \text{for } |x| < a. \end{cases} \qquad (94)$$

Since $V$ is negative, the momentum of the electron inside the barrier becomes imaginary when the following condition is

met:

$$-mc^2 < E + V < +mc^2, \qquad (95)$$

then the electron is quantum tunneling.

The analytic transmission for Eq. (94) is derived in Appendix D as:

$$T(k) = \frac{4\kappa^2}{4\kappa^2 + (1 - \kappa^2)^2 \sin^2\left(\frac{2mqa}{\hbar}\right)}, \qquad (96)$$

where I take $\kappa$ to be $\kappa_f$ in the FW representation, and $\kappa_d$ in the Dirac spinor representation. For the derivation of Appendix D in the FW representation

$$\kappa_f^{(\pm)} = \frac{\alpha_k^{(\pm)}}{\alpha_q^{(\pm)}} = \frac{q_x}{k_x} \frac{E \pm mc^2}{E \pm mc^2 + V}. \qquad (97)$$

It was shown by the authors of Refs. [16, 17] by deriving $T$ in the Dirac spinor representation of quantum tunneling

$$\kappa_d^{(\pm)} = \frac{q_x}{k_x} \frac{E \pm mc^2}{E \pm mc^2 + V}, \qquad (98)$$

which is in exact analytic agreement with $\kappa_f^{(\pm)}$, as expected.

The numerical results demonstrating relativistic tunneling of electrons are shown in Fig. 1. For comparison I also show the transmission of tunneling bosons with identical mass to the electron.

Fig. 1(b) shows the absolute error calculated as the difference between $T$ evaluated using the Taylor expansions $g(k/c)$ in $\kappa$ or the Python 3.8 `cmath` module square root function.

The region of slow convergence around $|k/c| = 1$ in Fig. 1(b) is an artifact of the point in $g(k/c)$ about which I have taken my Taylor expansions.

## 14.2. Tunneling through a barrier with vanishing electron mass inside the barrier

In this subsection I almost repeat the demonstration of Section 14.1 except that now inside the barrier the electron mass is strictly vanishing.

Since $m = 0$ when $|x| < a$, I must re-evaluate $\kappa_d$ for Eq. (96):

$$\kappa_d^{(\pm)} = \frac{q_x}{p_x} \frac{E \pm mc^2}{E + V}, \qquad (99)$$

which is in exact analytic agreement with $\kappa_f^{(\pm)}$, as expected.

The numerical results demonstrating tunneling of relativistic electrons are shown in Fig. 2. For comparison I also show the transmission of tunneling bosons with identical mass to the electron.

In Fig. 2(b) the absolute error is calculated as the difference between $T$ evaluated using the Taylor expansions $g(k/c)$ in $\kappa$ or the Python 3.8 `cmath` module square root function.

The region of slow convergence around $|k/c| = 1$ in Fig. 2(b) is an artifact of the point in $g(k/c)$ about which I have taken my Taylor expansions.
<figure><img src="image_5.png" /><figcaption>**Figure 3:** (color online) (a) Potential profile for a representation of a resonant-tunneling diode under a range of bias potentials $\Delta V/mc^2$ 0.00, 0.25, 0.50, 0.75, 1.00 as labelled. (b) Transmission as a function of bias potential $\Delta V$ for a series of left hand side electron incident energies $E^{(i)}/mc^2$ 0.400, 0.425, 0.450, 0.475 as labelled, with $\Gamma_i$ taken to be $\pm\alpha_i^{(\pm)}$. (c) Difference in transmission, as a function of bias potential $\Delta V$, between using $\pm\alpha_i^{(\pm)}$ or $p_i$ for $\Gamma_i$.</figcaption></figure>

### 14.3. Transmission through a resonant-tunneling diode using the WKB approximation

I calculate the transmission through a resonant diode under a range of bias potentials $\Delta V$. I show a series of potentials describing increasing $\Delta V$ in Fig. 3(a). The transmission as a function of $\Delta V$ is given in Fig. 3(b) for a series of incident energies. I calculate transmission with the two different WKB approximations, derived in Section 9 and Section 10, for comparison in Fig. 3(c).

The resonant-tunneling diode that I consider includes a planar cavity with a smoothly varying potential. For the wavefunction to traverse the cavity the employment of WKB transfer matrices developed in Section 9 and Section 10 are required:

$$ \left( \frac{\exp\left(-\frac{i}{\hbar} \int_{x_0}^{x_1} pdx\right)}{\sqrt{\Gamma_0/\Gamma_1}} \quad 0 \quad \left(\begin{array}{c} \tilde{a}_1 \\ \tilde{b}_1 \end{array}\right) = \left(\begin{array}{c} \tilde{a}_0 \\ \tilde{b}_0 \end{array}\right), \quad (100) $$

$\Gamma_i$ is $p_i$ or $\pm\alpha_i^{(\pm)}$ depending on whether I use FW Eq. (16) or Eq. (18) inside the cavity. The subscript 0 denotes evaluation at $x_0 = -a/2$ the left hand side of the smoothly varying region for which the transfer matrix is applied, 1 is at $x_1 = +a/2$ the right hand side of the smoothly varying region. I define $\tilde{a}_0$, $\tilde{a}_1$, $\tilde{b}_0$, and $\tilde{b}_1$ in the same way as Section 5.

Since I have established in Section 14.1 and Section 14.2 that the Taylor series for $g(k/c)$ gives convergence to the square root function of the Python 3.8 `cmath` module, I dispense with the Taylor series and use the `cmath` square root function to calculate the results in this subsection.

I see in Fig. 3(b) and (c) that as the energy of the incident electrons moves to higher values the first resonance of the diode moves towards lower values of $\Delta V$, as expected. I observe in Fig. 3(b) negative transmission as expected according to the Klein paradox [14, 16].

In Fig. 3(b) and (c) I see no qualitative difference between my two WKB approximations, showing that exponential decay is the dominant factor in the cavity for this example.

### 15. Conclusion

In this work I have developed Dirac fermion transfer matrix methods, for all energies, in the FW representations. I have also developed WKB approximations, to all orders, and discussed the validity of the approximations. I have introduced connection formulae between regions of real and imaginary fermion momentum.

I have shown the extension of my method to 2D periodic structures for band engineering [1, 2, 3].

I have alerted the reader to the limits of my methods by revealing a GF paradox that occurs when transforming between Dirac and FW representations.

I have verified my methods for massive and effectively massless fermions of all momentum by analytic comparison with the Dirac spinor representation of fermion tunneling. I have shown that the FW representation is an exact description of Dirac fermions and is not restricted to semi-relativistic energies as first envisioned [11].
<figure><img src="image_6.png" /><figcaption>**Figure 3:** (color online) (a) Potential profile for a representation of a resonant-tunneling diode under a range of bias potentials $\Delta V/mc^2$ 0.00, 0.25, 0.50, 0.75, 1.00 as labelled. (b) Transmission as a function of bias potential $\Delta V$ for a series of left hand side electron incident energies $E^{(i)}/mc^2$ 0.400, 0.425, 0.450, 0.475 as labelled, with $\Gamma_i$ taken to be $\pm\alpha_i^{(\pm)}$. (c) Difference in transmission, as a function of bias potential $\Delta V$, between using $\pm\alpha_i^{(\pm)}$ or $p_i$ for $\Gamma_i$.</figcaption></figure>

### 14.3. Transmission through a resonant-tunneling diode using the WKB approximation

I calculate the transmission through a resonant diode under a range of bias potentials $\Delta V$. I show a series of potentials describing increasing $\Delta V$ in Fig. 3(a). The transmission as a function of $\Delta V$ is given in Fig. 3(b) for a series of incident energies. I calculate transmission with the two different WKB approximations, derived in Section 9 and Section 10, for comparison in Fig. 3(c).

The resonant-tunneling diode that I consider includes a planar cavity with a smoothly varying potential. For the wavefunction to traverse the cavity the employment of WKB transfer matrices developed in Section 9 and Section 10 are required:

$$ \left( \frac{\exp\left(-\frac{i}{\hbar} \int_{x_0}^{x_1} p dx\right)}{\sqrt{\Gamma_0/\Gamma_1}} \quad 0 \quad \left( \frac{\exp\left(+\frac{i}{\hbar} \int_{x_0}^{x_1} p dx\right)}{\sqrt{\Gamma_1/\Gamma_0}} \right) \left( \begin{array}{c} \tilde{a}_1 \\ \tilde{b}_1 \end{array} \right) \right) = \left( \begin{array}{c} \tilde{a}_0 \\ \tilde{b}_0 \end{array} \right), \quad (100) $$

$\Gamma_i$ is $p_i$ or $\pm\alpha_i^{(\pm)}$ depending on whether I use FW Eq. (16) or Eq. (18) inside the cavity. The subscript 0 denotes evaluation at $x_0 = -a/2$ the left hand side of the smoothly varying region for which the transfer matrix is applied, 1 is at $x_1 = +a/2$ the right hand side of the smoothly varying region. I define $\tilde{a}_0$, $\tilde{a}_1$, $\tilde{b}_0$, and $\tilde{b}_1$ in the same way as Section 5.

Since I have established in Section 14.1 and Section 14.2 that the Taylor series for $g(k/c)$ gives convergence to the square root function of the Python 3.8 `cmath` module, I dispense with the Taylor series and use the `cmath` square root function to calculate the results in this subsection.

I see in Fig. 3(b) and (c) that as the energy of the incident electrons moves to higher values the first resonance of the diode moves towards lower values of $\Delta V$, as expected. I observe in Fig. 3(b) negative transmission as expected according to the Klein paradox [16, 18].

In Fig. 3(b) and (c) I see no qualitative difference between my two WKB approximations, showing that exponential decay is the dominant factor in the cavity for this example.

### 15. Conclusion

In this work I have developed Dirac fermion transfer matrix methods, for all energies, in the FW representations. I have also developed WKB approximations, to all orders, and discussed the validity of the approximations. I have introduced connection formulae between regions of real and imaginary fermion momentum.

I have demonstrated the applicability of my method to 2D periodic structures for band-gap engineering [1, 2, 3].

I have alerted the reader to the limits of my methods by revealing a GF paradox that occurs when transforming between Dirac and FW representations. In spite of the paradox, I have rigorously derived the 1D, 2D and 3D free space FW GFs. I note that since the FW Eq. (18) is a linear differential equation, the Dyson equations for the FW GFs are available for all types of Born approximation.

I have verified my methods and derivations for massive and effectively massless fermions of all momentum by analytic comparison with the Dirac spinor representation of fermion tunneling. I have shown that the FW representation is an exact description of Dirac fermions and is not restricted to semi-relativistic energies as first envisioned [11].
## 16. Acknowledgements

This research was not funded by any research council grants, if you want to help support this research please make a donation: https://gofund.me/34294631. Thank you in advance for whatever you donate.

Future developments will involve the inclusion of RSE perturbation theory [17, 18, 19, 20, 21]. I plan to extend the derivation of the normalization which I made in Ref. [18] to Dirac fermions described by the FW Eq. (18). I acknowledge that as a sub-sequence to my rigorous normalization derivations in Ref. [18], E. A. Muljarov incorporated his zero frequency mode discovery into the generalized normalization. I acknowledge the first of the RSE waveguide articles, which I reference, was predominantly the work of E. A. Muljarov.

## 17. Declaration of competing interests

I declare that I have no competing financial interests or personal relationships that could have influenced the work reported in this article.

### A. Discussion of commutation

In this appendix I will provide a discussion of

$$ \mathcal{L} = -\bar{\mathcal{L}}c\hat{p}_x = -c\hat{p}_x\bar{\mathcal{L}}. \quad (101) $$

Consider the form of the operator $\mathcal{L}$

$$ \mathcal{L} = mc^2 \left[ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 - \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 + \frac{1}{16} \left[ \frac{\hat{p}}{mc} \right]^6 - \dots \right], \quad (102) $$

where

$$ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 = - \frac{\hbar^2}{2m^2c^2} \left[ \frac{d^2}{dx^2} + \frac{d^2}{dy^2} \right], \quad (103) $$

$$ \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 = \frac{\hbar^4}{8m^4c^4} \left[ \frac{d^4}{dx^4} + 2 \frac{d^4}{dy^2 dx^2} + \frac{d^4}{dy^4} \right], \quad (104) $$

and the series of polynomial differential operators continues up to $n^{th}$ order. It is clear that I can write

$$ \frac{d^{n+m}}{dx^n dy^m} = \frac{d}{dx} \frac{d^{n+m-1}}{dx^{n-1} dy^{m-1}} = \frac{d}{dx^{n-1} dy^m} \frac{d}{dx}. \quad (105) $$

However it is also the case that

$$ \frac{d^n}{dy^n} f(x, y) = \frac{d}{dx} \frac{dx}{dy} \frac{d^{n-1} f}{dy^{n-1}} = \frac{d^{n-1}}{dy^{n-1}} \frac{dx}{dy} \frac{df}{dx} \quad (106) $$

when every term in $f$ is a function of not just $y$ but also $x$. I see from Eq. (105) and Eq. (106) that I can extract $d/dx$ from either side of the operator $\mathcal{L}$, in these cases.

When $f(x, y) = f(y)$, then in general

$$ \frac{d^n}{dy^n} f(y) \neq 0, \quad (107) $$

but

$$ \frac{d}{dx} f(y) = 0, \quad (108) $$

therefore in general

$$ 0 \neq \frac{d^n}{dy^n} f(y) \neq \frac{d^{n-1}}{dy^{n-1}} \frac{dx}{dy} \frac{d}{dx} f(y) = 0. \quad (109) $$

In the case described by Eq. (109) I cannot extract $d/dx$ from either side of the operator $d^n/dy^n$.

Eq. (101) is true when I assume that $\mathcal{L}$ acts on a function of not just $y$ but also $x$.

### B. Reflection at a step potential in Foldy-Wouthuysen and Dirac spinor representations for comparison

Consider a fermion incident on a step at $x = 0$. The wavefunction for $x < 0$ is

$$ \left( \frac{E \pm mc^2}{-c p_x} \right) \exp \left( \frac{i p_x x}{\hbar} \right) + B \left( \frac{E \pm mc^2}{+c p_x} \right) \exp \left( -\frac{i p_x x}{\hbar} \right), \quad (110) $$

and for $x > 0$

$$ F \left( \frac{E + V \pm mc^{*2}}{-c^* q_x} \right) \exp \left( \frac{i q_x x}{\hbar} \right). \quad (111) $$

By continuity of wavefunction at $x = 0$ I have

$$ (1 - B) c p_x = F c^* q_x, \quad (112) $$

$$ (1 + B)(E \pm mc^2) = F(E + V \pm mc^{*2}). \quad (113) $$

Dividing Eq. (112) by Eq. (113) gives

$$ \frac{1 - B}{1 + B} = \frac{q_x c^*}{p_x c} \frac{E \pm mc^2}{E + V \pm mc^{*2}} = \kappa_d^{(\pm)}. \quad (114) $$

Now I turn to the FW representation. The corresponding wavefunction for $x < 0$ is

$$ \exp \left( \frac{i p_x x}{\hbar} \right) + B \exp \left( -\frac{i p_x x}{\hbar} \right), \quad (115) $$

and for $x > 0$

$$ F \exp \left( \frac{i q_x x}{\hbar} \right). \quad (116) $$

Applying my boundary condition at $x = 0$, continuity of wavefunction and continuity of wavefunction with respect to $\bar{\mathcal{L}}$, I obtain

$$ (1 - B) \alpha_0^{(\pm)} = F \alpha_1^{(\pm)}, \quad (117) $$

$$ 1 + B = F. \quad (118) $$
## 16. Acknowledgements

This research was not funded by any research council grants. If you would like to support this research, please make a donation: https://gofund.me/34294631. Contributions to support this research are greatly appreciated.

Future developments will involve the inclusion of RSE perturbation theory [12, 13, 19, 20, 21]. I plan to extend the derivation of the normalization which I made in Ref. [20] to Dirac fermions described by the FW Eq. (18). I acknowledge that as a sub-sequence to my rigorous normalization derivations in Ref. [20], E. A. Muljarov incorporated his zero frequency mode discovery into the generalized normalization. I acknowledge the first of the RSE waveguide articles Ref. [21], was predominantly the work of E. A. Muljarov.

## 17. Declaration of competing interests

I declare that I have no competing financial interests or personal relationships that could have influenced the work reported in this article.

### A. Discussion of commutation

In this appendix I will provide a discussion of

$$ \mathcal{L} = -\hat{\mathcal{L}}c\hat{p}_x = -c\hat{p}_x\hat{\mathcal{L}}. \quad (101) $$

The operator $\mathcal{L}$ takes the form

$$ \mathcal{L} = mc^2 \left[ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 - \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 + \frac{1}{16} \left[ \frac{\hat{p}}{mc} \right]^6 - \dots \right], \quad (102) $$

where

$$ \frac{1}{2} \left[ \frac{\hat{p}}{mc} \right]^2 = - \frac{\hbar^2}{2m^2c^2} \left[ \frac{d^2}{dx^2} + \frac{d^2}{dy^2} \right], \quad (103) $$

$$ \frac{1}{8} \left[ \frac{\hat{p}}{mc} \right]^4 = \frac{\hbar^4}{8m^4c^4} \left[ \frac{d^4}{dx^4} + 2 \frac{d^4}{dy^2 dx^2} + \frac{d^4}{dy^4} \right], \quad (104) $$

and the series of polynomial differential operators continues up to $n^{th}$ order. It is clear that I can write

$$ \frac{d^{n+m}}{dx^n dy^m} = \frac{d}{dx} \frac{d^{n+m-1}}{dx^{n-1} dy^{m-1}} = \frac{d^{n+m-1}}{dx^{n-1} dy^m} \frac{d}{dx}. \quad (105) $$

However, it is also the case that

$$ \frac{d^n}{dy^n} f(x, y) = \frac{d}{dx} \frac{dx}{dy} \frac{d^{n-1} f}{dy^{n-1}} = \frac{d^{n-1}}{dy^{n-1}} \frac{dx}{dy} \frac{df}{dx} \quad (106) $$

when every term in $f$ is a function of not just $y$ but also $x$. I see from Eq. (105) and Eq. (106) that I can extract $d/dx$ from either side of the operator $\mathcal{L}$, in these cases.

When $f(x, y) = f(y)$, then in general

$$ \frac{d^n}{dy^n} f(y) \neq 0, \quad (107) $$

but

$$ \frac{d}{dx} f(y) = 0, \quad (108) $$

therefore in general

$$ 0 \neq \frac{d^n}{dy^n} f(y) \neq \frac{d^{n-1}}{dy^{n-1}} \frac{dx}{dy} \frac{d}{dx} f(y) = 0. \quad (109) $$

In the case described by Eq. (109) I cannot extract $d/dx$ from either side of the operator $d^n/dy^n$.

Eq. (101) is true when I assume that $\mathcal{L}$ acts on a function of not just $y$ but also $x$.

### B. Reflection at a step potential in Foldy-Wouthuysen and Dirac spinor representations for comparison

Consider a fermion incident on a step at $x = 0$. The wavefunction for $x < 0$ is

$$ \left( \frac{E \pm mc^2}{-c p_x} \right) \exp \left( \frac{i p_x x}{\hbar} \right) + B \left( \frac{E \pm mc^2}{+c p_x} \right) \exp \left( -\frac{i p_x x}{\hbar} \right), \quad (110) $$

and for $x > 0$

$$ F \left( \frac{E + V \pm mc^{*2}}{-c^* q_x} \right) \exp \left( \frac{i q_x x}{\hbar} \right). \quad (111) $$

Using wavefunction continuity at $x = 0$, I find:

$$ (1 - B) c p_x = F c^* q_x, \quad (112) $$

$$ (1 + B)(E \pm mc^2) = F(E + V \pm mc^{*2}). \quad (113) $$

Dividing Eq. (112) by Eq. (113) gives

$$ \frac{1 - B}{1 + B} = \frac{q_x c^*}{p_x c} \frac{E \pm mc^2}{E + V \pm mc^{*2}} = \kappa_d^{(\pm)}. \quad (114) $$

Now I turn to the FW representation. The corresponding wavefunction for $x < 0$ is

$$ \exp \left( \frac{i p_x x}{\hbar} \right) + B \exp \left( -\frac{i p_x x}{\hbar} \right), \quad (115) $$

and for $x > 0$

$$ F \exp \left( \frac{i q_x x}{\hbar} \right). \quad (116) $$

Applying my boundary condition at $x = 0$, continuity of wavefunction and continuity of wavefunction with respect to $\hat{\mathcal{L}}$, I obtain

$$ (1 - B) \alpha_0^{(\pm)} = F \alpha_1^{(\pm)}, \quad (117) $$

$$ 1 + B = F. \quad (118) $$
Dividing Eq. (117) by Eq. (118) gives

$$ \frac{1-B}{1+B} = \frac{\alpha_1^{(\pm)}}{\alpha_0^{(\pm)}} = \frac{p_x c}{q_x c^*} \frac{E+V \pm mc^{*2}}{E \pm mc^2} = \kappa_f^{(\pm)}. \quad (119) $$

In Section 3 I derived $R = |B|^2$, therefore

$$ R = \left| \frac{1 - \kappa}{1 + \kappa} \right|^2, \qquad (120) $$

where I use $\kappa_d^{(\pm)}$ for the Dirac equation and $\kappa_f^{(\pm)}$ for the FW equation. Substituting either $\kappa_f^{(\pm)}$ or $\kappa_d^{(\pm)}$ into Eq. (120) gives analytically identical T and R between the Dirac and FW representation, as expected. I also calculated $\kappa_d^{(\pm)}$ in the case of m vanishing for $x > 0$ to find exact analytic agreement in T and R between the Dirac and FW representations, as expected.

$\kappa_d^{(\pm)}$ corresponds to $\alpha^{(\pm)}$, and is the mathematical consequence of there being two forms for the FW wavefunction Eq. (18). The FW wavefunction describes the electron $\tilde{\Psi}^{(-)}$ or positron $\tilde{\Psi}^{(+)}$ [11].

I also calculate, for useful comparison, the transmission of bosons with the following Klein-Gordon result

$$ \kappa_{kg} = \frac{q_x c^*}{p_x c}. \qquad (121) $$

### C. Examination of a massless fermion tunneling through a magnetic delta barrier in the Foldy-Wouthuysen representation

Consider a massless fermion tunneling through a magnetic delta barrier $B = B_0\delta(x)\hat{z}$, where $B = \nabla \times A$. The barrier is described by $A = 0$ except for $A = A_y\hat{y} \neq 0$ when $x > 0$. The Dirac spinor wavefunction is a solution of

$$ [\alpha c \cdot (\hat{p} - A)] \Psi = E \Psi. \qquad (122) $$

Using the approach of Section 2, Eq. (122) can be transformed to

$$ [\mp c \sqrt{(\hat{p} - A)^2}] \tilde{\Psi}^{(\pm)} = E \tilde{\Psi}^{(\pm)}. \qquad (123) $$

Taylor expanding Eq. (123) around $\hat{p} = 0$ gives

$$ [\mp M - (E \pm cA_y)] \tilde{\Psi}^{(\pm)} = 0, \qquad (124) $$

where $M$ is the corresponding linear differential operator. Using the approach of Section 4 I calculate, that at the barrier, the wavefunction is continuous and also continuous with respect to $\tilde{M}$, where I define

$$ \tilde{M} \tilde{\Psi}^{(\pm)} = \pm \frac{E \pm c A_y}{c p_x} \tilde{\Psi}^{(\pm)}. \qquad (125) $$

I now return to the Dirac spinor representation. In the region $A = A_y\hat{y} \neq 0$, the general solution for $\hat{p} = \pm q_x\hat{x}$ is given by [22, 23, 24]

$$ \Psi = \begin{pmatrix} +E \\ \pm c q_x - i c A_y \end{pmatrix} \exp \left( \pm \frac{i q_x x}{\hbar} \right), \qquad (126) $$

where $c q_x = \sqrt{E^2 - c^2 A_y^2}$. The reflection is calculated by matching wavefunctions at the step in vector potential, similar to Appendix C for a step in electrostatic potential.

Making use of my boundary conditions I find the reflection coefficient R takes the form of Eq. (120), with

$$ \kappa_d = \frac{c^* q_x - i A_y}{c p_x} \qquad (127) $$

for the Dirac spinor representation and

$$ \kappa_f = \frac{c p_x}{c^* q_x} \frac{E \pm c^* A_y}{E} \qquad (128) $$

for the FW representation. $\hat{p} = p_x\hat{x} = (E/c)\hat{x}$ when $A = 0$.

The conflict between Eq. (127) and Eq. (128) shows that R calculated in the FW representation does not reproduce R calculated in the Dirac spinor representation, FW wavefunctions cannot be used for the calculation of tunneling through magnetic barriers.

### D. Analytic Green's functions in 1-dimension

Consider

$$ [\mp \mathcal{L} - (V + E \pm mc^2)] G^{(\pm)}(x, x') = \delta(x - x'). \qquad (129) $$

Let $\tilde{\Psi}_L$ and $\tilde{\Psi}_R$ be the solutions of Eq. (18) which separately satisfy the boundary conditions at $x = -a$ and $x = +a$. The GF of Eq. (129) for the case of a 1D barrier can be written as

$$ G^{(\pm)}(x, x') = \frac{\tilde{\Psi}_{L}^{(\pm)}(x_c) \tilde{\Psi}_{R}^{(\pm)}(x_c)}{\tilde{W}^{(\pm)}}, \qquad (130) $$

where $x_c = \min(x, x')$ and $x_c = \max(x, x')$ and

$$ \tilde{W}^{(\pm)} = \mp \tilde{\Psi}_{L}^{(\pm)} i c \hbar \tilde{\mathcal{L}} \tilde{\Psi}_{R}^{(\pm)} \pm \tilde{\Psi}_{R}^{(\pm)} i c \hbar \tilde{\mathcal{L}} \tilde{\Psi}_{L}^{(\pm)} \qquad (131) $$

is the Wronskian, which does not depend on x.

To prove Eq. (130) and Eq. (131), note that by Eq. (129)

$$ \mp \lim_{\epsilon \to 0} \int_{x'-\epsilon}^{x'+\epsilon} \mathcal{L} G^{(\pm)}(x, x') dx = 1, \qquad (132) $$

which implies

$$ \mp i c \hbar \tilde{\mathcal{L}} G^{(\pm)}(x_+, x') \pm i c \hbar \tilde{\mathcal{L}} G^{(\pm)}(x_-, x') = 1, \qquad (133) $$

then, applying $\tilde{\mathcal{L}}$ to Eq. (130) written out explicitly

$$ \tilde{\mathcal{L}} G(x, x') = \frac{1}{\tilde{W}} \begin{cases} \tilde{\mathcal{L}} \tilde{\Psi}_L(x) \tilde{\Psi}_R(x') & \text{for } x < x', \\ \tilde{\Psi}_L(x') \tilde{\mathcal{L}} \tilde{\Psi}_R(x) & \text{for } x > x', \end{cases} \qquad (134) $$

and substituting into Eq. (133) gives Eq. (131).
Dividing Eq. (117) by Eq. (118) gives

$$ \frac{1-B}{1+B} = \frac{\alpha_1^{(\pm)}}{\alpha_0^{(\pm)}} = \frac{p_x c}{q_x c^*} \frac{E+V \pm mc^{*2}}{E \pm mc^2} = \kappa_f^{(\pm)}. \quad (119) $$

In Section 3 I derived $R = |B|^2$, therefore

$$ R = \left| \frac{1 - \kappa}{1 + \kappa} \right|^2, \quad (120) $$

where I use $\kappa_d^{(\pm)}$ for the Dirac equation and $\kappa_f^{(\pm)}$ for the FW equation. Substituting either $\kappa_f^{(\pm)}$ or $\kappa_d^{(\pm)}$ into Eq. (120) gives analytically identical T and R between the Dirac and FW representation, as expected. I also calculated $\kappa_d^{(\pm)}$ in the case of m vanishing for $x > 0$ and found exact analytic agreement in T and R between the Dirac and FW representations, as expected.

I calculate, for useful comparison, the transmission of bosons with the following Klein-Gordon result

$$ \kappa_{kg} = \frac{q_x c^*}{p_x c}. \quad (121) $$

### C. Examination of a massless fermion tunneling through a magnetic delta barrier in the Foldy-Wouthuysen representation

Consider a massless fermion tunneling through a magnetic delta barrier $B = B_0\delta(x)\hat{z}$, where $B = \nabla \times A$. The barrier is described by $A = 0$ except for $A = A_y\hat{y} \neq 0$ when $x > 0$. The Dirac spinor wavefunction is a solution of

$$ [c\alpha \cdot (\hat{p} - A)] \Psi = E\Psi. \quad (122) $$

Using the approach of Section 2, Eq. (122) can be transformed to

$$ [\mp c \sqrt{(\hat{p} - A)^2}] \tilde{\Psi}^{(\pm)} = E \tilde{\Psi}^{(\pm)}. \quad (123) $$

Taylor expanding Eq. (123) around $\hat{p} = 0$ gives

$$ [\mp M - (E \pm cA_y)] \tilde{\Psi}^{(\pm)} = 0, \quad (124) $$

where $M$ is the corresponding linear differential operator. Using the approach of Section 4 I calculate, that at the barrier, the wavefunction is continuous and also continuous with respect to $\tilde{M}$, where I define

$$ \tilde{M} \tilde{\Psi}^{(\pm)} = \pm \frac{E \pm c A_y}{c p_x} \tilde{\Psi}^{(\pm)}. \quad (125) $$

I now return to the Dirac spinor representation. In the region $A = A_y\hat{y} \neq 0$, the general solution for $\hat{p} = \pm q_x \hat{x}$ is given by [22, 23]

$$ \Psi = \begin{pmatrix} +E \\ \pm c q_x - i c A_y \end{pmatrix} \exp \left( \pm \frac{i q_x x}{\hbar} \right), \quad (126) $$

where $c q_x = \sqrt{E^2 - c^2 A_y^2}$. The reflection is calculated by matching wavefunctions at the step in vector potential, similar to Appendix C for a step in electrostatic potential.

Making use of my boundary conditions I find the reflection coefficient R takes the form of Eq. (120), with

$$ \kappa_d = \frac{c^* q_x - i A_y}{c p_x} \quad (127) $$

for the Dirac spinor representation and

$$ \kappa_f = \frac{c p_x}{c^* q_x} \frac{E \pm c^* A_y}{E} \quad (128) $$

for the FW representation. $\hat{p} = p_x \hat{x} = (E/c)\hat{x}$ when $A = 0$.

The conflict between Eq. (127) and Eq. (128) indicates R calculated in the FW representation does not reproduce R calculated in the Dirac spinor representation. This suggests that FW wavefunctions are unsuitable for calculating tunneling through magnetic barriers. However, this conflict might be resolved by employing the form of the FW transformation in Ref. [11] that derives the Pauli equation from the Dirac equation in the non-relativistic limit.

### D. Analytic Green's functions in one dimension

Consider

$$ [\mp \mathcal{L} - (V + E \pm mc^2)] G^{(\pm)}(x, x') = \delta(x - x'). \quad (129) $$

Let $\tilde{\Psi}_L$ and $\tilde{\Psi}_R$ be the solutions of Eq. (18) which separately satisfy the boundary conditions at $x = -a$ and $x = +a$. The GF of Eq. (129) for the case of a 1D barrier can be written as

$$ G^{(\pm)}(x, x') = \frac{\tilde{\Psi}_{L}^{(\pm)}(x_{<}) \tilde{\Psi}_{R}^{(\pm)}(x_{>})}{W^{(\pm)}}, \quad (130) $$

where $x_{<} = \min(x, x')$ and $x_{>} = \max(x, x')$ and

$$ W^{(\pm)} = \mp \tilde{\Psi}_{L}^{(\pm)} i c \hbar \tilde{\mathcal{L}} \tilde{\Psi}_{R}^{(\pm)} \pm \tilde{\Psi}_{R}^{(\pm)} i c \hbar \tilde{\mathcal{L}} \tilde{\Psi}_{L}^{(\pm)} \quad (131) $$

is the Wronskian, which does not depend on x.

To prove Eq. (130) and Eq. (131), note that by Eq. (129)

$$ \mp \lim_{\epsilon \to 0} \int_{x'-\epsilon}^{x'+\epsilon} \mathcal{L} G^{(\pm)}(x, x') dx = 1, \quad (132) $$

which implies

$$ \mp i c \hbar \tilde{\mathcal{L}} G^{(\pm)}(x_+, x') \pm i c \hbar \tilde{\mathcal{L}} G^{(\pm)}(x_-, x') = 1, \quad (133) $$

then, applying $\tilde{\mathcal{L}}$ to Eq. (130) written out explicitly

$$ \tilde{\mathcal{L}} G(x, x') = \frac{1}{W} \begin{cases} \tilde{\mathcal{L}} \tilde{\Psi}_L(x) \tilde{\Psi}_R(x') & \text{for } x < x', \\ \tilde{\Psi}_L(x') \tilde{\mathcal{L}} \tilde{\Psi}_R(x) & \text{for } x > x', \end{cases} \quad (134) $$
I now calculate the GF for a barrier of width $2a$.

$$ \tilde{\Psi}_1 = \begin{cases} \exp\left( + \frac{imkx}{\hbar} \right) & x > a, \\ C \exp\left( + \frac{im^*qx}{\hbar} \right) + D \exp\left( - \frac{im^*qx}{\hbar} \right) & |x| < a, \\ A \exp\left( + \frac{imkx}{\hbar} \right) + B \exp\left( - \frac{imkx}{\hbar} \right) & x < -a. \end{cases} \quad (135) $$

Continuity of $\Psi_1(x)$ and continuity of $\Psi_1(x)$ with respect to $\tilde{\mathcal{L}}$ gives

$$ C = \frac{\alpha_q + \alpha_k}{2\alpha_q} \exp\left( + \frac{imka}{\hbar} \right) \exp\left( - \frac{im^*qa}{\hbar} \right), \quad (136) $$

$$ D = \frac{\alpha_q - \alpha_k}{2\alpha_q} \exp\left( + \frac{imka}{\hbar} \right) \exp\left( + \frac{im^*qa}{\hbar} \right). \quad (137) $$

$\tilde{\Psi}_2$ the left hand solution is given by $\tilde{\Psi}_2(+x) = \tilde{\Psi}_1(-x)$. Therefore the Wronskian is given by

$$ \bar{W} = 2ic\hbar\alpha_q (C^2 - D^2). \quad (138) $$

I evaluate Eq. (130) to give $G(+a, -a)$:

$$ \frac{i\alpha_q}{ic\hbar(\alpha_q^2 + \alpha_k^2) \sin\left(\frac{2m^*qa}{\hbar}\right) - 2c\hbar\alpha_q\alpha_k \cos\left(\frac{2m^*qa}{\hbar}\right)}. \quad (139) $$

The transmission for a barrier embedded in a vacuum is given as the ratio of $G(+a, -a)$ and the free space GF $G^{fs}(+a, -a)$:

$$ T(k) = |2c\hbar\alpha_k G(+a, -a)|^2 = \left| \frac{G(+a, -a)}{G^{fs}(+a, -a)} \right|^2 . \quad (140) $$

The free space GF can be written as

$$ [\mp \mathcal{L} - (E \pm mc^2)] G^{(\pm)}(x, x') = \delta(x - x'), \quad (141) $$

therefore by symmetry, when $x' = 0$,

$$ G^{(\pm)}(x, 0) = A^{(\pm)} \begin{cases} \exp\left(-\frac{imkx}{\hbar}\right) & \text{for } x < 0, \\ \exp\left(+\frac{imkx}{\hbar}\right) & \text{for } x > 0. \end{cases} \quad (142) $$

Integrating Eq. (141) across $[0_-, 0_+]$ gives

$$ \mp ic\hbar A^{(\pm)} \bar{\mathcal{L}} \exp\left( + \frac{imk0_+}{\hbar} \right) \pm ic\hbar A^{(\pm)} \bar{\mathcal{L}} \exp\left( - \frac{imk0_-}{\hbar} \right) = 1 \quad (143) $$

and letting $0_+ \to 0$ and $0_- \to 0$, I find

$$ A^{(\pm)} = \frac{1}{2ic\hbar\alpha_k^{(\pm)}}. \quad (144) $$

## E. Green's function paradox

I reveal a GF paradox occurring when transforming between the FW and Dirac representations.

Let

$$ [H - (E + V)] \mathcal{G}(r, r') = \delta(r - r') \hat{1}, \quad (145) $$

therefore

$$ \begin{aligned} & e^{+iS} [H - (E + V)] e^{-iS} e^{+iS} [\mathcal{G}(r, r')] e^{-iS} \\ & = e^{+iS} \hat{1} e^{-iS} \delta(r - r'). \end{aligned} \quad (146) $$

However since

$$ e^{+iS} \hat{1} e^{-iS} \delta(r - r') = \hat{1} \delta(r - r'), \quad (147) $$

and since I have already seen in Section 2

$$ e^{+iS} [H - (E + V)] e^{-iS} = [\mp \mathcal{L} - (V + E \pm mc^2)] \hat{1}, \quad (148) $$

I have

$$ \begin{aligned} & e^{+iS} [\mathcal{G}(r, r')] e^{-iS} = G^{(\pm)}(r, r') \hat{1} \\ & \implies e^{-iS} G^{(\pm)}(r, r') e^{+iS} = \mathcal{G}(r, r') \hat{1}. \end{aligned} \quad (149) $$

Eqs. (149) imply that the valid GFs of the Dirac equation are diagonal. Now consider the identity

$$ [H - (E + V)][H + (E + V)] K(r, r') = \hat{1} \delta(r - r'), \quad (150) $$

so that

$$ [H + (E + V)] K(r, r') = \mathcal{G}(r, r'). \quad (151) $$

Multiplying out the brackets in Eq. (150) I obtain

$$ -c^2 \hbar^2 \left[ \nabla^2 + \left(\frac{p}{\hbar}\right)^2 \right] K(r, r') = \delta(r - r'). \quad (152) $$

I see from Eq. (152) that in 1D free space ($V = 0$):

$$ c^2 \hbar^2 K(x, x') = - \frac{\hbar \exp\left(\frac{ip |x - x'|}{\hbar}\right)}{2ip}, \quad (153) $$

in 2D free space:

$$ c^2 \hbar^2 K(\rho, \rho') = \frac{i}{4} H_0^{(1)} \left( \frac{p |\rho - \rho'|}{\hbar} \right), \quad (154) $$

and in 3D free space:

$$ c^2 \hbar^2 K(r, r') = \frac{\exp\left(\frac{ip |r - r'|}{\hbar}\right)}{4\pi |r - r'|}. \quad (155) $$
and substituting into Eq. (133) gives Eq. (131).

I now calculate the GF for a rectangular barrier of width $2a$.

$$ \tilde{\Psi}_1 = \begin{cases} \exp\left(\frac{imkx}{\hbar}\right) & x > +a, \\ C \exp\left(\frac{im^*qx}{\hbar}\right) + D \exp\left(-\frac{im^*qx}{\hbar}\right) & |x| < a, \\ A \exp\left(\frac{imkx}{\hbar}\right) + B \exp\left(-\frac{imkx}{\hbar}\right) & x < -a. \end{cases} \quad (135) $$

Continuity of $\tilde{\Psi}_1(x)$ and continuity of $\tilde{\Psi}_1(x)$ with respect to $\mathcal{L}$ gives

$$ C = \frac{\alpha_q + \alpha_k}{2\alpha_q} \exp\left(\frac{imka}{\hbar}\right) \exp\left(-\frac{im^*qa}{\hbar}\right), \quad (136) $$

$$ D = \frac{\alpha_q - \alpha_k}{2\alpha_q} \exp\left(\frac{imka}{\hbar}\right) \exp\left(\frac{im^*qa}{\hbar}\right). \quad (137) $$

$\tilde{\Psi}_2$ the left hand solution is given by $\tilde{\Psi}_2(+x) = \tilde{\Psi}_1(-x)$. Therefore I evaluate Eq. (130) to give $G(+a, -a)$:

$$ \frac{i\alpha_q}{ic\hbar \left(\alpha_q^2 + \alpha_k^2\right) \sin\left(\frac{2m^*qa}{\hbar}\right) - 2c\hbar\alpha_q\alpha_k \cos\left(\frac{2m^*qa}{\hbar}\right)}. \quad (138) $$

In Eq. (138), the Wronskian is given by

$$ W = 2ic\hbar\alpha_q (C^2 - D^2). \quad (139) $$

The transmission for a barrier embedded in a vacuum is given as the ratio of $G(+a, -a)$ and the free space GF $G^{fs}(+a, -a)$:

$$ T(k) = |2c\hbar\alpha_k G(+a, -a)|^2 = \left| \frac{G(+a, -a)}{G^{fs}(+a, -a)} \right|^2. \quad (140) $$

The free space GF can be written as

$$ [\mp \mathcal{L} - (E \pm mc^2)] G^{(\pm)}(x, x') = \delta(x - x'), \quad (141) $$

therefore by symmetry, when $x' = 0$,

$$ G^{(\pm)}(x, 0) = A^{(\pm)} \begin{cases} \exp\left(-\frac{imkx}{\hbar}\right) & \text{for } x < 0, \\ \exp\left(\frac{imkx}{\hbar}\right) & \text{for } x > 0. \end{cases} \quad (142) $$

Integrating Eq. (141) across $[0_-, 0_+]$ gives

$$ \mp ic\hbar A^{(\pm)} \mathcal{L} \exp\left(\frac{imk0_+}{\hbar}\right) \pm ic\hbar A^{(\pm)} \mathcal{L} \exp\left(-\frac{imk0_-}{\hbar}\right) = 1 \quad (143) $$

and letting $0_+ \to 0$ and $0_- \to 0$, I find

$$ A^{(\pm)} = \frac{1}{2ic\hbar\alpha_k^{(\pm)}}. \quad (144) $$

## E. Green's function paradox

I reveal a GF paradox occurring when transforming between the FW and Dirac representations.

Let

$$ [H - (E + V)] G(r, r') = \delta(r - r') \hat{1}, \quad (145) $$

therefore

$$ \begin{aligned} & e^{+iS} [H - (E + V)] e^{-iS} e^{+iS} [G(r, r')] e^{-iS} \\ & \qquad = e^{+iS} \hat{1} e^{-iS} \delta(r - r'). \end{aligned} \quad (146) $$

However since

$$ e^{+iS} \hat{1} e^{-iS} \delta(r - r') = \hat{1} \delta(r - r'), \quad (147) $$

and since I have already seen in Section 2

$$ e^{+iS} [H - (E + V)] e^{-iS} = [\mp \mathcal{L} - (V + E \pm mc^2)] \hat{1}, \quad (148) $$

I have

$$ \begin{aligned} & e^{+iS} [G(r, r')] e^{-iS} = G^{(\pm)}(r, r') \hat{1} \\ & \qquad \Rightarrow e^{-iS} G^{(\pm)}(r, r') e^{+iS} = G(r, r') \hat{1}. \end{aligned} \quad (149) $$

Eqs. (149) imply that the valid GFs of the Dirac equation are diagonal. Now consider the identity

$$ [H - (E + V)][H + (E + V)] K(r, r') = \hat{1} \delta(r - r'), \quad (150) $$

so that

$$ [H + (E + V)] K(r, r') = G(r, r'). \quad (151) $$

Multiplying out the brackets in Eq. (150) I obtain

$$ -c^2 \hbar^2 \left[ \nabla^2 + \left(\frac{p}{\hbar}\right)^2 \right] K(r, r') = \delta(r - r'). \quad (152) $$

I see from Eq. (152) that in 1D free space ($V = 0$):

$$ c^2 \hbar^2 K(x, x') = - \frac{\hbar \exp\left(\frac{ip |x - x'|}{\hbar}\right)}{2ip}, \quad (153) $$

in 2D free space:

$$ c^2 \hbar^2 K(\rho, \rho') = \frac{i}{4} H_0^{(1)} \left( \frac{p |\rho - \rho'|}{\hbar} \right), \quad (154) $$

and in 3D free space:

$$ c^2 \hbar^2 K(r, r') = \frac{\exp\left(\frac{ip |r - r'|}{\hbar}\right)}{4\pi |r - r'|}. \quad (155) $$
Substituting the identity Eq. (151) into Eq. (149)

$$e^{+iS} [H + (E + V)] e^{-iS} K(r, r') = G^{(\pm)}(r, r') \hat{1} \quad (156)$$

and so

$$\left[ \sqrt{\hat{\rho}^2 c^2 + m^2 c^4} + (E + V) \right] K(r, r') = G^{(\pm)}(r, r') \hat{1}. \quad (157)$$

By inspection of Eq. (157) and the energy momentum relation $(E + V) = \sqrt{p^2 c^2 + m^2 c^4}$,

$$G^{(\pm)}(r, r') = 2(E + V) K(r, r') \quad (158)$$

and substituting Eq. (158) into Eq. (149) I also have

$$G(r, r') = 2(E + V) K(r, r'). \quad (159)$$

Eq. (159) is consistent with the relativistic correction to the 1st Born approximation in 3D.

Let us test Eq. (158) and Eq. (159) by making an analytic comparison between the GF for 1D free space Eq. (142) and the GF composed of Eq. (159) and Eq. (153):

$$2EK = -E \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{i\hbar pc^2} \neq \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{2ic\hbar\alpha_k} \quad (160)$$

and so I have an apparent paradox.

Let us attempt to resolve the paradox. I have shown in Appendix D that there are two possible forms of the FW GF:

$$G^{(\pm)}(x, x') = \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{2ic\hbar\alpha_k^{(\pm)}} \quad (161)$$

for an electron delta source (−) or a positron source (+). If I add my two expressions $G^{(\pm)}(x, x')$ I find

$$G^{(-)}(x, x') + G^{(+)}(x, x') = 2EK(x, x'). \quad (162)$$

However if I try to prove that in general

$$G(r, r') = G^{(-)}(r, r') + G^{(+)}(r, r'), \quad (163)$$

I arrive at a contradiction. The Klein-Gordon GF has the continuity conditions of the Klein-Gordon equation, $(E + V)$ can be a discontinuous function, therefore $2(E + V)K$ can be discontinuous at boundaries. Therefore Eq. (159) implies $G$ may be discontinuous. $G$ has the same boundary conditions as $\Psi$, Eq. (159) may violate the boundary conditions for $G$.

Let us now discuss a possible cause of the paradox. To evaluate the FW transformations defined to be

$$e^{\pm iS} = \cos \theta \pm \beta\alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta, \quad (164)$$

I am required to know the momentum $\hat{p}$ also at the delta source. To make this calculation of $\hat{p}$ at the delta source consider Eq. (142) from Appendix D,

$$G(x, 0) = \frac{1}{2ic\hbar\alpha_k} \begin{cases} \exp\left(-\frac{imkx}{\hbar}\right) & \text{for } x < 0, \\ \exp\left(+\frac{imkx}{\hbar}\right) & \text{for } x > 0. \end{cases} \quad (165)$$

I see from Eq. (165) that

$$p_x(x = 0) G(x, 0) = \lim_{0_{\pm} \to 0} \hat{p}G(0_{\pm}, 0) = \pm mkG(0, 0). \quad (166)$$

Eq. (166) indicates that when $x = x'$, the momentum $\hat{p}$ is undefined. Therefore at the delta source the FW transformation is undefined. In 1D it might be impossible to evaluate $e^{\pm iS}$ at the delta source.

## References

[1] M. Barbier, P. Vasilopoulos, F. Peeters, Extra dirac points in the energy spectrum for superlattices on single-layer graphene, Physical Review B 81 (2010) 075438. doi:https://doi.org/10.1103/PhysRevB.81.075438.

[2] M. Barbier, F. Peeters, P. Vasilopoulos, J. M. Pereira, Dirac and klein-gordon particles in one-dimensional periodic potentials, Phys. Rev. B 77 (2008) 115446. doi:https://doi.org/10.1103/PhysRevB.77.115446.

[3] D. Arovas, L. Brey, H. Fertig, E.-A. Kim, K. Ziegler, Dirac spectrum in piecewise constant one-dimensional (1d) potentials, New Journal of Physics 12 (2010) 123020. doi:https://doi.org/10.1088/1367-2630/12/12/123020.

[4] F. Guinea, M. Katsnelson, A. Geim, Energy gaps and a zero-field quantum hall effect in graphene by strain engineering, Nature Physics 6 (2010) 30–33. doi:https://doi.org/10.1038/nphys1420.

[5] K. S. Novoselov, A. Mishchenko, A. Carvalho, A. H. C. Neto, 2d materials and van der waals heterostructures, Science 353 (2016). doi:https://doi.org/10.1126/science.aac9439.

[6] C. Wan, X. Gu, F. Dang, T. Itoh, Y. Wang, H. Sasaki, M. Kondo, K. Koga, K. Yabuki, G. J. Snyder, et al., Flexible n-type thermoelectric materials by organic intercalation of layered transition metal dichalcogenide tis 2, Nature materials 14 (2015) 622–627. doi:https://doi.org/10.1038/nmat4251.

[7] Y. Song, H.-C. Wu, Y. Guo, Negative differential resistances in graphene double barrier resonant tunneling diodes, Applied Physics Letters 102 (2013) 093118. doi:https://doi.org/10.1063/1.4794952.

[8] E. Sonin, Effect of klein tunneling on conductance and shot noise in ballistic graphene, Physical Review B 79 (2009) 195438. doi:https://doi.org/10.1103/PhysRevB.79.195438.

[9] T. Tudorovskiy, K. Reijnders, M. I. Katsnelson, Chiral tunneling in single-layer and bilayer graphene, Physica Scripta 2012 (2012) 014010. doi:https://doi.org/10.1088/0031-8949/2012/T146/014010.

[10] V. Zalipaev, C. Linton, M. Croitoru, A. Vagov, Resonant tunneling and localized states in a graphene monolayer with a mass gap, Physical Review B 91 (2015) 085405. doi:https://doi.org/10.1103/PhysRevB.91.085405.

[11] L. L. Foldy, S. A. Wouthuysen, On the dirac theory of spin 1/2 particles and its non-relativistic limit, Physical Review 78 (1950) 29. doi:https://doi.org/10.1103/PhysRev.78.29.

[12] D. Y. K. Ko, J. Sambles, Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals, JOSA A 5 (1988) 1863–1866. doi:https://doi.org/10.1364/JOSAA.5.001863.

[13] S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, T. Ishihara, Quasiguided modes and optical properties of photonic crystal slabs, Phys. Rev. B 66 (2002) 045102. URL: https://doi.org/10.1103/PhysRevB.66.045102.
Substituting the identity Eq. (151) into Eq. (149)

$$e^{+iS} [H + (E + V)] e^{-iS} K(r, r') = G^{(\pm)}(r, r') \hat{1} \quad (156)$$

and so

$$[\sqrt{\hat{p}^2 c^2 + m^2 c^4} + (E + V)] K(r, r') = G^{(\pm)}(r, r') \hat{1}. \quad (157)$$

By inspection of Eq. (157) and the energy momentum relation $(E + V) = \sqrt{p^2 c^2 + m^2 c^4}$,

$$G^{(\pm)}(r, r') = 2(E + V) K(r, r') \quad (158)$$

and substituting Eq. (158) into Eq. (149) I also have

$$\mathcal{G}(r, r') = 2(E + V) K(r, r'). \quad (159)$$

Let us test Eq. (158) and Eq. (159) by making an analytic comparison between the GF for 1D free space Eq. (142) and the GF composed of Eq. (158) and Eq. (153):

$$2EK = -E \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{i\hbar pc^2} \neq \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{2ic\hbar\alpha_k} \quad (160)$$

and so I have an apparent paradox.

Let us attempt to resolve the paradox. I have shown in Appendix D that there are two possible forms of the FW GF:

$$G^{(\pm)}(x, x') = \frac{\exp\left(\frac{ip|x-x'|}{\hbar}\right)}{2ic\hbar\alpha_k^{(\pm)}} \quad (161)$$

for an electron delta source (−) or a positron source (+). If I add these two expressions $G^{(\pm)}(x, x')$ I find

$$G^{(-)}(x, x') + G^{(+)}(x, x') = 2EK(x, x'). \quad (162)$$

However if I try to prove that in general

$$\mathcal{G}(r, r') = G^{(-)}(r, r') + G^{(+)}(r, r'), \quad (163)$$

I arrive at a contradiction. The Klein-Gordon GF has the continuity conditions of the Klein-Gordon equation, $(E + V)$ can be a discontinuous function, therefore $2(E + V)K$ can be discontinuous at boundaries. Therefore Eq. (159) implies $\mathcal{G}$ may be discontinuous. $\mathcal{G}$ has the same boundary conditions as $\Psi$, Eq. (159) may violate the boundary conditions for $\mathcal{G}$.

Let us now discuss a possible cause of the paradox. To evaluate the FW transformations defined to be

$$e^{\pm iS} = \cos \theta \pm \beta\alpha \cdot \frac{\hat{p}}{|\hat{p}|} \sin \theta, \quad (164)$$

I am required to know the momentum $\hat{p}$ also at the delta source. To make this calculation of $\hat{p}$ at the delta source consider Eq. (142) from Appendix D,

$$G(x, 0) = \frac{1}{2ic\hbar\alpha_k} \begin{cases} \exp\left(-\frac{imkx}{\hbar}\right) & \text{for } x < 0, \\ \exp\left(+\frac{imkx}{\hbar}\right) & \text{for } x > 0. \end{cases} \quad (165)$$

I see from Eq. (165) that

$$p_x(x = 0) G(x, 0) = \lim_{0_{\pm} \to 0} \hat{p} G(0_{\pm}, 0) = \pm mkG(0, 0). \quad (166)$$

Eq. (166) indicates that when $x = x'$, the momentum $\hat{p}$ is undefined. Therefore at the delta source the FW transformation is undefined. For 1D space it might be impossible to evaluate $e^{\pm iS}$ at the delta source.

## F. Correct free space Foldy-Wouthuysen Green's functions in two and three dimensions

I begin with the 3D FW GF in free space

$$G_{3D}^{(\pm)}(r, r') = - \frac{p}{c\hbar^2 \alpha_k^{(\pm)}} \frac{1}{4\pi |r-r'|} e^{ip|r-r'|/\hbar}, \quad (167)$$

which at the end of this Appendix F, I will show to be consistent with the proven 1D FW GF, Eq. (142) and Eq. (144). Eq. (167) is the solution of

$$[\mp \mathcal{L} - (E \pm mc^2)] G_{3D}^{(\pm)}(r, r') = \delta_{3D}(r - r'). \quad (168)$$

Eq. (167) can be thought of as an outgoing exponential plane FW wavefunction spreading out over a sphere centered on the source, as we move further from the source. In Eq. (167) I make a choice of amplitude to obtain agreement with the correct 1D FW GF.

To derive the 1D free space FW GF, I integrate Eq. (167) in the plane $(y, z)$ and find

$$G_{1D}^{(\pm)}(x, x') = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} G_{3D}^{(\pm)}(r, r') dy dz, \quad (169)$$

where $r = (x, y, z)$. Eq. (169) is a consequence of

$$\delta_{1D}(x - x') = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \delta_{3D}(r - r') dy dz. \quad (170)$$

Let $\rho = \sqrt{y^2 + z^2}$ in the cylindrical coordinates $(\rho, \theta, x)$, and write $\xi = \sqrt{(x - x')^2 + \rho^2}$, so that Eq. (169) becomes

$$G_{1D}^{(\pm)}(x, x') = - \frac{p}{c\hbar^2 \alpha_k^{(\pm)}} \int_0^{2\pi} \int_0^{\infty} \frac{1}{4\pi\xi} e^{ip\xi/\hbar} \rho d\rho d\theta. \quad (171)$$

Integrating Eq. (171) with respect to $\theta$ and changing variable to $\xi$ gives, using the chain rule

$$d\xi = \frac{\rho}{\sqrt{(x-x')^2 + \rho^2}} d\rho = \frac{\rho}{\xi} d\rho, \quad (172)$$

noting in the integration limits that I have $\xi = |x - x'|$ as a minimum when $\rho = 0$:

$$G_{1D}^{(\pm)}(x, x') = - \frac{p}{c\hbar^2 \alpha_k^{(\pm)}} \int_{|x-x'|}^{\infty} \frac{1}{2} e^{ip\xi/\hbar} d\xi. \quad (173)$$
link.aps.org/doi/10.1103/PhysRevB.66.045102. doi:10.1103/PhysRevB.66.045102.

[14] A. Calogeracos, N. Dombey, History and physics of the klein paradox, Contemporary physics 40 (1999) 313–321. doi:https://doi.org/10.1080/001075199181387.

[15] H. Dosch, V. Muller, J. Jensen, Kleins paradox, Physica Norvegica 5 (1971) 151.

[16] O. Klein, Die reflexion von elektronen an einem potentialsprung nach der relativistischen dynamik von dirac, Zeitschrift für Physik 53 (1929) 157–165. doi:https://doi.org/10.1007/BF01339716.

[17] E. A. Muljarov, W. Langbein, R. Zimmermann, Brillouin-wigner perturbation theory in open electromagnetic systems, EPL (Europhysics Letters) 92 (2011) 50010. doi:https://doi.org/10.1209/0295-5075/92/50010.

[18] M. B. Doost, W. Langbein, E. A. Muljarov, Resonant-state expansion applied to three-dimensional open optical systems, Physical Review A 90 (2014) 013834. doi:https://doi.org/10.1103/PhysRevA.90.013834.

[19] L. J. Armitage, M. B. Doost, W. Langbein, E. A. Muljarov, Resonant-state expansion applied to planar waveguides, Physical Review A 89 (2014) 053832. doi:https://doi.org/10.1103/PhysRevA.89.053832.

[20] M. B. Doost, Resonant-state-expansion born approximation with a correct eigen-mode normalisation, Journal of Optics 18 (2016) 085607. doi:https://doi.org/10.1088/2040-8978/18/8/085607.

[21] M. B. Doost, Resonant-state-expansion born approximation for waveguides with dispersion, Physical Review A 93 (2016) 023835. doi:https://doi.org/10.1103/PhysRevA.93.023835.

[22] M. R. Masir, P. Vasilopoulos, A. Matulis, F. Peeters, Direction-dependent tunneling through nanostructured magnetic barriers in graphene, Physical Review B 77 (2008) 235443. doi:https://doi.org/10.1103/PhysRevB.77.235443.

[23] M. R. Masir, P. Vasilopoulos, F. Peeters, Wavevector filtering through single-layer and bilayer graphene with magnetic barrier structures, Applied Physics Letters 93 (2008) 242103. doi:https://doi.org/10.1063/1.3049600.

[24] M. R. Masir, P. Vasilopoulos, F. Peeters, Magnetic kronig–penney model for dirac electrons in single-layer graphene, New Journal of Physics 11 (2009) 095009. doi:https://doi.org/10.1088/1367-2630/11/9/095009.
Eq. (173) can be integrated while dropping the oscillating upper integration limit at $\xi = \infty$, assuming the contributions at $\xi = \infty$ cancel out leading to a zero net contribution,

$$ \int_{|x-x'|}^{\infty} \frac{1}{2} e^{ip\xi/\hbar} d\xi = \left[ \frac{\hbar}{2ip} e^{ip\xi/\hbar} \right]_{|x-x'|}^{\infty} = -\frac{\hbar}{2ip} e^{ip|x-x'|/\hbar}. \quad (174) $$

Hence I arrive at the required result, Eq. (142) combined with Eq. (144):

$$ G_{1D}^{(\pm)}(x, x') = \frac{p}{c\hbar^2\alpha_k^{(\pm)}} \frac{\hbar}{2ip} e^{ip|x-x'|/\hbar} = \frac{e^{ip|x-x'|/\hbar}}{2ic\hbar\alpha_k^{(\pm)}}. \quad (175) $$

The integration in Eq. (174) can be seen clearly when I take $p$ to be complex with even only a vanishingly small imaginary part so that the upper integration limit is an exponentially decaying function, rather than just an oscillating function.

An interesting comparison can be made between the FW GF Eqs. (167) and Eqs. (175) and their Klein-Gordon counterparts, Eqs. (155) and Eqs. (153), which were derived in Appendix E. The derivation in this Appendix F can be repeated, as a simple checking exercise, for Eqs. (155) and Eqs. (153). These comparisons and derivations, repeated in 1D and 2D space, suggest that the 2D confined free space FW GF is given by

$$ G_{2D}^{(\pm)}(r, r') = -\frac{p}{c\hbar^2\alpha_k^{(\pm)}} \frac{i}{4} H_0^{(1)}(p|r-r'|/\hbar). \quad (176) $$

Eq. (167) and Eq. (176) can be simplified with

$$ \frac{p}{c\hbar^2\alpha_k^{(\pm)}} = -\frac{E \mp mc^2}{c^2\hbar^2}. \quad (177) $$

I note that since the FW Eq. (18) is a linear differential equation, the Dyson equations for the FW GFs are available for all types of Born approximation [12, 13]. There might be cases where the decoupling of electron and positron states by the FW approach could simplify or eliminate certain divergences in the Dyson series, potentially addressing aspects of the renormalization problem.

## References

[1] M. Barbier, P. Vasilopoulos, F. Peeters, Extra dirac points in the energy spectrum for superlattices on single-layer graphene, Physical Review B 81 (2010) 075438. doi:https://doi.org/10.1103/PhysRevB.81.075438.

[2] M. Barbier, F. Peeters, P. Vasilopoulos, J. M. Pereira, Dirac and klein-gordon particles in one-dimensional periodic potentials, Phys. Rev. B 77 (2008) 115446. doi:https://doi.org/10.1103/PhysRevB.77.115446.

[3] D. Arovas, L. Brey, H. Fertig, E.-A. Kim, K. Ziegler, Dirac spectrum in piecewise constant one-dimensional (1d) potentials, New Journal of Physics 12 (2010) 123020. doi:https://doi.org/10.1088/1367-2630/12/12/123020.

[4] F. Guinea, M. Katsnelson, A. Geim, Energy gaps and a zero-field quantum hall effect in graphene by strain engineering, Nature Physics 6 (2010) 30–33. doi:https://doi.org/10.1038/nphys1420.

[5] K. S. Novoselov, A. Mishchenko, A. Carvalho, A. H. C. Neto, 2d materials and van der waals heterostructures, Science 353 (2016). doi:https://doi.org/10.1126/science.aac9439.

[6] C. Wan, X. Gu, F. Dang, T. Itoh, Y. Wang, H. Sasaki, M. Kondo, K. Koga, K. Yabuki, G. J. Snyder, et al., Flexible n-type thermoelectric materials by organic intercalation of layered transition metal dichalcogenide tis 2, Nature materials 14 (2015) 622–627. doi:https://doi.org/10.1038/nmat4251.

[7] Y. Song, H.-C. Wu, Y. Guo, Negative differential resistances in graphene double barrier resonant tunneling diodes, Applied Physics Letters 102 (2013) 093118. doi:https://doi.org/10.1063/1.4794952.

[8] E. Sonin, Effect of klein tunneling on conductance and shot noise in ballistic graphene, Physical Review B 79 (2009) 195438. doi:https://doi.org/10.1103/PhysRevB.79.195438.

[9] T. Tudorovskiy, K. Reijnders, M. I. Katsnelson, Chiral tunneling in single-layer and bilayer graphene, Physica Scripta 2012 (2012) 014010. doi:https://doi.org/10.1088/0031-8949/2012/T146/014010.

[10] V. Zalipaev, C. Linton, M. Croitoru, A. Vagov, Resonant tunneling and localized states in a graphene monolayer with a mass gap, Physical Review B 91 (2015) 085405. doi:https://doi.org/10.1103/PhysRevB.91.085405.

[11] L. L. Foldy, S. A. Wouthuysen, On the dirac theory of spin 1/2 particles and its non-relativistic limit, Physical Review 78 (1950) 29. doi:https://doi.org/10.1103/PhysRev.78.29.

[12] M. B. Doost, Resonant-state-expansion born approximation with a correct eigen-mode normalisation, Journal of Optics 18 (2016) 085607. doi:https://doi.org/10.1088/2040-8978/18/8/085607.

[13] M. B. Doost, Resonant-state-expansion born approximation for waveguides with dispersion, Physical Review A 93 (2016) 023835. doi:https://doi.org/10.1103/PhysRevA.93.023835.

[14] D. Y. K. Ko, J. Sambles, Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals, JOSA A 5 (1988) 1863–1866. doi:https://doi.org/10.1364/JOSAA.5.001863.

[15] S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, T. Ishihara, Quasiguided modes and optical properties of photonic crystal slabs, Phys. Rev. B 66 (2002) 045102. doi:https://doi.org/10.1103/PhysRevB.66.045102.

[16] A. Calogeracos, N. Dombey, History and physics of the klein paradox, Contemporary physics 40 (1999) 313–321. doi:https://doi.org/10.1080/001075199181387.

[17] H. Dosch, V. Muller, J. Jensen, Kleins paradox, Physica Norvegica 5 (1971) 151.

[18] O. Klein, Die reflexion von elektronen an einem potentialsprung nach der relativistischen dynamik von dirac, Zeitschrift für Physik 53 (1929) 157–165. doi:https://doi.org/10.1007/BF01339716.

[19] E. A. Muljarov, W. Langbein, R. Zimmermann, Brillouin-wigner perturbation theory in open electromagnetic systems, EPL (Europhysics Letters) 92 (2011) 50010. doi:https://doi.org/10.1209/0295-5075/92/50010.

[20] M. B. Doost, W. Langbein, E. A. Muljarov, Resonant-state expansion applied to three-dimensional open optical systems, Physical Review A 90 (2014) 013834. doi:https://doi.org/10.1103/PhysRevA.90.013834.

[21] L. J. Armitage, M. B. Doost, W. Langbein, E. A. Muljarov, Resonant-state expansion applied to planar waveguides, Physical Review A 89 (2014) 053832. doi:https://doi.org/10.1103/PhysRevA.89.053832.

[22] M. R. Masir, P. Vasilopoulos, A. Matulis, F. Peeters, Direction-dependent tunneling through nanostructured magnetic barriers in graphene, Physical Review B 77 (2008) 235443. doi:https://doi.org/10.1103/PhysRevB.77.235443.

[23] M. R. Masir, P. Vasilopoulos, F. Peeters, Wavevector filtering through single-layer and bilayer graphene with magnetic barrier structures, Applied Physics Letters 93 (2008) 242103. doi:https://doi.org/10.1063/1.3049600.