Quantum mechanical spin is regularly described as a phenomenon that arises naturally only in relativistic physics. Here for example is a quote from *What Remains to Be Discovered*, a book written by a former editor of *Nature*:

Pauli was the first to guess at an explanation: spin, he argued, exists because of the way the time dimension must be added to the three dimensions of ordinary space to provide a relativistic description of objects. Spin is nature’s way of signaling the correctness of Einstein’s theory of relativity. [pp.72-3]

This sort of perspective is unfortunately not entirely correct, as explained over fifty years ago by Jean-Marc Lévy-Leblond. Spin can be shown to arise in a natural way in the completely nonrelativistic context of the Schrödinger equation. The presentation of this surprising insight is based on this paper, with a slightly different derivation.

The method used by Lévy-Leblond is to apply the traditional derivation of the Dirac equation to the Schrödinger equation. Dirac had the idea to derive first-order wave equations by factoring the second-order Klein-Gordon equation:

$[\frac{1}{{c}^{2}}\frac{{\partial}^{2}}{\partial {t}^{2}}-\sum _{k}\frac{{\partial}^{2}}{\partial {x}_{k}^{2}}+(\frac{mc}{\hslash}{)}^{2}]\psi =0$

The sum here runs over all spatial variables in the system to which the equation is applied. In practice that is three dimensions, but it is worth keeping in mind that the following procedure can be carried out for other choices using Clifford algebras.

Physicists generally set $\hslash =c=1$ in this equation to simplify algebra. Explicit constants are retained here to facilitate comparison with the nonrelativistic derivation to come.

To get the Dirac equation, form the linear combinations with constant coefficients

$\begin{array}{l}{\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}+i\frac{mc}{\hslash}\\ {\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}-i\frac{mc}{\hslash}\end{array}$

where repeated spatial indices represent sums over those indices. Multiply these combinations:

$\begin{array}{l}({\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}+i\frac{mc}{\hslash})({\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{l}\frac{\partial}{\partial {x}_{l}}-i\frac{mc}{\hslash})\\ \phantom{\rule{3em}{0ex}}=({\gamma}^{0}{)}^{2}\frac{1}{{c}^{2}}\frac{{\partial}^{2}}{\partial {t}^{2}}+({\gamma}^{k}{)}^{2}\frac{{\partial}^{2}}{\partial {x}_{k}^{2}}+(\frac{mc}{\hslash}{)}^{2}\\ \phantom{\rule{6em}{0ex}}+({\gamma}^{0}{\gamma}^{k}+{\gamma}^{k}{\gamma}^{0})\frac{1}{c}\frac{{\partial}^{2}}{\partial t\partial {x}_{k}}+({\gamma}^{k}{\gamma}^{l}+{\gamma}^{l}{\gamma}^{k})\frac{{\partial}^{2}}{\partial {x}_{k}\partial {x}_{l}}\end{array}$

The multiplicative order of the constant coefficients with different indices is deliberate: Dirac realized that he could factor the Klein-Gordon equation if the coefficients are noncommuting matrices with the following properties:

$\begin{array}{l}({\gamma}^{0}{)}^{2}=1\\ ({\gamma}^{k}{)}^{2}=-1\end{array}\phantom{\rule{3em}{0ex}}\begin{array}{l}{\gamma}^{0}{\gamma}^{k}+{\gamma}^{k}{\gamma}^{0}=0\\ \phantom{\rule{.3em}{0ex}}{\gamma}^{k}{\gamma}^{l}+{\gamma}^{l}{\gamma}^{k}=0\phantom{\rule{.5em}{0ex}},\phantom{\rule{.5em}{0ex}}k\ne l\end{array}$

With these choices, the Dirac equation is either of

$\begin{array}{l}[{\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}+i\frac{mc}{\hslash}]\psi =0\\ [{\gamma}^{0}\frac{1}{c}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}-i\frac{mc}{\hslash}]\psi =0\end{array}$

While the first choice is conventional, either sign on the mass leads to equivalent physics.

There are several explicit representations of the gamma matrices, depending on the emphasis of the particular solution to the Dirac equation. For the standard solution in one temporal and three spatial dimensions, 4×4 matrices are necessary to accommodate all four variables. In the Dirac basis these matrices are

${\gamma}^{0}=\left(\begin{array}{cc}I& 0\\ 0& -I\end{array}\right)\phantom{\rule{5em}{0ex}}{\gamma}^{k}=\left(\begin{array}{cc}0& {\sigma}_{k}\\ -{\sigma}_{k}& 0\end{array}\right)$

where the submatrices are the 2×2 Pauli matrices

${\sigma}_{1}=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\phantom{\rule{3em}{0ex}}{\sigma}_{2}=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right)\phantom{\rule{3em}{0ex}}{\sigma}_{3}=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)$

supplemented by a 2×2 identity matrix. The Dirac basis is convenient for describing low-energy systems, and thus useful for comparison with the Schrödinger equation. Another basis is that of Weyl,

${\gamma}^{0}=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\phantom{\rule{5em}{0ex}}{\gamma}^{k}=\left(\begin{array}{cc}0& -{\sigma}_{k}\\ {\sigma}_{k}& 0\end{array}\right)$

which is more suited to describing high-energy systems. The choice of sign for ${\gamma}^{k}$ is related to the form of a fifth gamma matrix not needed here.

The point of factoring the Klein-Gordon equation is to allow new details to emerge. This can be illustrated simply by comparing the second-order differential equation

$(\frac{{d}^{2}}{d{x}^{2}}+1)f=0$

to the corresponding factored equation:

$(\frac{d}{dx}+i)(\frac{d}{dx}-i)f=0$

The solutions to the first equation are $sinx$ and $cosx$ . The factored form of the equation picks out the constituent functions making up the circular functions: the first factor has only the solution ${e}^{-ix}$ , while the second has only the solution ${e}^{ix}$ . These two functions are related by a reversal of sign of the independent variable that picks out odd and even parts. Further, the factors act as projection operators on the constituent function space, selecting only the single function that satisfies each factor separately.

Similarly, the two forms of the Dirac equation act as projection operators on the space of solutions to the Klein-Gordon equation, bringing out detail that washes out in forming the second-order equation. And just as for the simple factoring example, there will be a relation connecting solutions to the two forms of the Dirac equation.

Nontrivial details are expected to arise for nonrelativistic systems from a parallel factoring process.

For physical context, first consider applying equations above to plane wave states. For a mass without spin, one can use a scalar with a Lorentz-invariant exponent,

$\psi (x,p)=exp[-\frac{i}{\hslash}px]=exp\left[\frac{i}{\hslash}(\mathbf{p}\xb7\mathbf{x}-Et)\right]$

where the negative sign on energy ensures the correct relationship to nonrelativistic kinetic energy in the Schrödinger equation. Acting upon this function with the Klein-Gordon equation leads to

${E}^{2}={\mathbf{p}}^{2}{c}^{2}+{m}^{2}{c}^{4}$

which is simply the definition of the relativistic energy of a free particle. For small momenta, one has the nonrelativistic approximation

$E=m{c}^{2}\sqrt{1+\frac{{\mathbf{p}}^{2}}{{m}^{2}{c}^{2}}}\approx m{c}^{2}+\frac{{\mathbf{p}}^{2}}{2m}$

For the four-dimensional Dirac equation, plane waves must include two-dimensional bispinors

$\psi (x,p)=\left[\begin{array}{c}u\left(p\right)\\ v\left(p\right)\end{array}\right]exp\left[\frac{i}{\hslash}(\mathbf{p}\xb7\mathbf{x}-Et)\right]$

that are functions of energy and momentum. Applying the first form of the Dirac equation to this gives

$({\gamma}^{0}E-c{\gamma}^{k}{p}_{k}-m{c}^{2})\left[\begin{array}{c}u\\ v\end{array}\right]=0$

which in the Dirac basis, suitable for low-energy phenomena, becomes

$\left[\begin{array}{cc}E-m{c}^{2}& -c\mathbf{\sigma}\xb7\mathbf{p}\\ c\mathbf{\sigma}\xb7\mathbf{p}& -E-m{c}^{2}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=0$

which is a matrix equation for the four components of the two bispinors. Expanding into component equations one has

$(E-m{c}^{2})u=c\mathbf{\sigma}\xb7\mathbf{p}v\phantom{\rule{5em}{0ex}}c\mathbf{\sigma}\xb7\mathbf{p}u=(E+m{c}^{2})v$

Now introduce an electromagnetic field into this relativistic system with a minimal coupling in Gaussian units

$E\to E-e\phi \phantom{\rule{5em}{0ex}}\mathbf{p}\to \mathbf{p}-\frac{e}{c}\mathbf{A}$

so that the component equations become

$(E-m{c}^{2}-e\phi )u=c\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})v\phantom{\rule{3em}{0ex}}c\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})u=(E+m{c}^{2}-e\phi )v$

For small momentum and electrostatic potential, the rightmost coefficient can be approximated as

$E+m{c}^{2}-e\phi \approx 2m{c}^{2}+\frac{{\mathbf{p}}^{2}}{2m}-e\phi \approx 2m{c}^{2}$

and then eliminating the bispinor $v$ leads to the single nonrelativistic equation

$[\frac{1}{2m}\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})+e\phi ]u=(E-m{c}^{2})u$

where the factor on the right-hand side is energy apart from rest mass. Using the identity of the Pauli matrices for the product of two dot products,

$(\mathbf{\sigma}\xb7\mathbf{a})(\mathbf{\sigma}\xb7\mathbf{b})=\mathbf{a}\xb7\mathbf{b}+i\mathbf{\sigma}\xb7\mathbf{a}\times \mathbf{b}$

the corresponding product in square brackets becomes

$\begin{array}{l}\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})\\ \phantom{\rule{4em}{0ex}}=(\mathbf{p}-\frac{e}{c}\mathbf{A}{)}^{2}+i\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})\times (\mathbf{p}-\frac{e}{c}\mathbf{A})\\ \phantom{\rule{4em}{0ex}}=(\mathbf{p}-\frac{e}{c}\mathbf{A}{)}^{2}-\frac{ie}{c}\mathbf{\sigma}\xb7[\mathbf{p}\times \mathbf{A}+\mathbf{A}\times \mathbf{p}]\end{array}$

If the momentum in the second term is treated as a quantum mechanical operator, and remembering that cross products are antisymmetric, then one can write

$\begin{array}{l}[\mathbf{p}\times \mathbf{A}+\mathbf{A}\times \mathbf{p}]u=\left[\right(\frac{\hslash}{i}\nabla \times \mathbf{A})-(\mathbf{A}\times \frac{\hslash}{i}\nabla )+(\mathbf{A}\times \frac{\hslash}{i}\nabla \left)\right]u\\ \phantom{[\mathbf{p}\times \mathbf{A}+\mathbf{A}\times \mathbf{p}]u}=\left[\right(\frac{\hslash}{i}\nabla \times \mathbf{A}\left)\right]u=\frac{\hslash}{i}\mathbf{B}u\end{array}$

where the last step uses the definition of the magnetic field in terms of the vector potential. The single nonrelativistic equation is then

$[\frac{1}{2m}(\mathbf{p}-\frac{e}{c}\mathbf{A}{)}^{2}+e\phi -\frac{e\hslash}{2mc}\mathbf{\sigma}\xb7\mathbf{B}]u=(E-m{c}^{2})u$

The first two terms in square brackets are the Hamiltonian of a charged particle in a external field. The third term is the spin-orbit coupling of a charged particle with spin of one half, *i.e.*, an electron.

This nonrelativistic equation is known as the Pauli equation. It appeared a year before Dirac’s derivation of the full relativistic equation, but was considered ad hoc in its construction. The Dirac equation is notable for producing the correct value of the spin term automatically.

The Schrödinger equation for a freely moving mass is

$[i\hslash \frac{\partial}{\partial t}+\frac{{\hslash}^{2}}{2m}\sum _{k}\frac{{\partial}^{2}}{\partial {x}_{k}^{2}}]\psi =0$

Multiplying by 2*m* and applying to plane waves gives

$(2mE-{\mathbf{p}}^{2})\psi =0$

The factorization will be done in terms of these physical variables rather than derivatives, simply because the notation is more compact.

The relevant section of Lévy-Leblond’s paper, “III b) Linearization of the Schrödinger Equation”, is unfortunately not quite as clear as it could be, and contains several obvious typographical errors. There is also an arbitrariness in the selection of coefficient matrices as part of the process.

Rather than attacking the factorization of the Schrödinger equation directly, consider reverse engineering the nonrelativistic approximation to the Dirac equation. This provides a nonarbitrary choice for coefficient matrices, and although different from the choice in the paper, will produce the expected physics. This is related to the freedom of basis available in the gamma matrices.

Here again is the Dirac equation applied to plane waves in the Dirac basis, before any approximations:

$\left[\begin{array}{cc}{E}_{\mathrm{rel}}-m{c}^{2}& -c\mathbf{\sigma}\xb7\mathbf{p}\\ c\mathbf{\sigma}\xb7\mathbf{p}& -{E}_{\mathrm{rel}}-m{c}^{2}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=0$

To convert this to the nonrelativistic case, replace the diagonal terms with nonrelativistic quantities and omit the speed of light as nonphysical:

$\left[\begin{array}{cc}E& -\mathbf{\sigma}\xb7\mathbf{p}\\ \mathbf{\sigma}\xb7\mathbf{p}& -2m\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=0$

Expand this into component equations

$Eu=\mathbf{\sigma}\xb7\mathbf{p}v\phantom{\rule{5em}{0ex}}\mathbf{\sigma}\xb7\mathbf{p}u=2mv$

and introduce an electromagnetic field using the same minimal coupling:

$(E-e\phi )u=\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})v\phantom{\rule{3em}{0ex}}\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})u=2mv$

These can be combined without approximation into a single equation

$[\frac{1}{2m}\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})\mathbf{\sigma}\xb7(\mathbf{p}-\frac{e}{c}\mathbf{A})+e\phi ]u=Eu$

The appearance of the speed of light here is a result of how Gaussian units are defined, not some relativistic effect. The product of dot products is treated as above,

$[\frac{1}{2m}(\mathbf{p}-\frac{e}{c}\mathbf{A}{)}^{2}+e\phi -\frac{e\hslash}{2mc}\mathbf{\sigma}\xb7\mathbf{B}]u=Eu$

resulting in the same nonrelativistic Pauli equation with automatic spin-orbit coupling. This is the crux of the appearance of spin in a natural way mentioned at the outset.

To complete the factoring process, rewrite the nonrelativistic equation in terms of gamma matrices in the Dirac basis:

$[{\gamma}^{0}\frac{E+2m}{2}-{\gamma}^{k}{p}_{k}+\frac{E-2m}{2}]\left[\begin{array}{c}u\\ v\end{array}\right]=0$

This is the first factor, representing the Schrödinger equation in a four-dimensional form. Following Lévy-Leblond, there must be some other factor that when applied on the left produces the usual Schrödinger equation. Write this second factor with arbitrary coefficients and multiply the two using properties of the gamma matrices:

$\begin{array}{l}(A{\gamma}^{0}+B{\gamma}^{k}{p}_{k}+C)({\gamma}^{0}\frac{E+2m}{2}-{\gamma}^{k}{p}_{k}+\frac{E-2m}{2})\\ \phantom{\rule{3em}{0ex}}=A\frac{E+2m}{2}+B{\mathbf{p}}^{2}+C\frac{E-2m}{2}\\ \phantom{\rule{6em}{0ex}}-{\gamma}^{0}{\gamma}^{k}{p}_{k}(A+B\frac{E+2m}{2})+{\gamma}^{0}(A\frac{E-2m}{2}+C\frac{E+2m}{2})\\ \phantom{\rule{6em}{0ex}}+{\gamma}^{k}{p}_{k}(B\frac{E-2m}{2}-C)\end{array}$

Comparing with the Schrödinger equation applied to plane waves immediately gives

$A=\frac{E+2m}{2}\phantom{\rule{4em}{0ex}}B=-1\phantom{\rule{4em}{0ex}}C=-\frac{E-2m}{2}$

since these choices make the last three terms identically zero, and the first and third terms sum to $2mE$ . Given the number of equations involved, this simple result implies quite a bit of symmetry.

Restoring derivatives, the Schrödinger equation in four-dimensional form is one of either

$\begin{array}{l}\left[{\gamma}^{0}\right(\frac{i\hslash}{2}\frac{\partial}{\partial t}+m)+i\hslash {\gamma}^{k}\frac{\partial}{\partial {x}_{k}}+\frac{i\hslash}{2}\frac{\partial}{\partial t}-m]\psi =0\\ \left[{\gamma}^{0}\right(\frac{i\hslash}{2}\frac{\partial}{\partial t}+m)+i\hslash {\gamma}^{k}\frac{\partial}{\partial {x}_{k}}-\frac{i\hslash}{2}\frac{\partial}{\partial t}+m]\psi =0\end{array}$

or simplifying a bit

$\begin{array}{l}[\frac{{\gamma}^{0}+1}{2}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}-i\frac{m}{\hslash}({\gamma}^{0}-1)]\psi =0\\ [\frac{{\gamma}^{0}-1}{2}\frac{\partial}{\partial t}+{\gamma}^{k}\frac{\partial}{\partial {x}_{k}}-i\frac{m}{\hslash}({\gamma}^{0}+1)]\psi =0\end{array}$

The two forms are related in the same way as the two forms of the Dirac equation, with a change in sign of mass but essentially the same physics.

As a final note, Lévy-Leblond points out that one can include an arbitrary invertible matrix between the two factors of the squared Schrödinger equation. This is equivalent to multiplying the first form of the equation by this same matrix. The first and third coefficients in that equation are

$\frac{{\gamma}^{0}+1}{2}=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\phantom{\rule{5em}{0ex}}{\gamma}^{0}-1=\left(\begin{array}{cc}0& 0\\ 0& -2\end{array}\right)$

Reading past typographical errors, the invertible matrix used by Lévy-Leblond is $\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)$ , which is a two-dimensional representation of the imaginary unit. Applying this matrix to the coefficient matrices gives

$\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)\phantom{\rule{3em}{0ex}}\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{cc}0& 0\\ 0& -2\end{array}\right)=\left(\begin{array}{cc}0& 2\\ 0& 0\end{array}\right)$

which are the coefficients *A* and *C*, apart from mass, chosen by Lévy-Leblond. The second coefficient follows similarly,

$\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right){\gamma}^{k}=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\left(\begin{array}{cc}0& {\sigma}_{k}\\ -{\sigma}_{k}& 0\end{array}\right)=\left(\begin{array}{cc}{\sigma}_{k}& 0\\ 0& {\sigma}_{k}\end{array}\right)$

providing a connection between the somewhat arbitrarily chosen coefficients in the paper and those employed in this presentation.

*Uploaded 2023.05.31*
analyticphysics.com