Optical Dirac equation

The properties of an electron in relativistic quantum theory are represented by operators, just as in the more familiar non-relativistic theory. The position and momentum operators, ${\bf x}$ and ${\bf p}$ satisfy the canonical commutation relation

$$\begin{eqnarray}&&[{{x}_{j}},{{p}_{k}}]=i\hbar {{\delta }_{jk}}\end{eqnarray} \tag{ 17 }$$

and the Dirac equation (12) leads directly to the Hamiltonian operator

$$\begin{eqnarray}&&H={\boldsymbol{ \alpha }} \cdot {\bf p}+\beta m.\end{eqnarray} \tag{ 18 }$$

This, together with the Heisenberg equation of motion, provides the evolution equations for the operators representing mechanical properties of interest. The position operator, for example, satisfies the equation

$$\begin{eqnarray}&&{\bf \dot{x}}=\frac{i}{\hbar }[H,{\bf x}]={\boldsymbol{ \alpha }} ,\end{eqnarray} \tag{ 19 }$$

which we recognize as the velocity operator. The momentum operator commutes with the Hamiltonian and so, unlike the velocity operator, it is a constant of the motion

$$\begin{eqnarray}&&{\bf \dot{p}}=\frac{i}{\hbar }[H,{\bf p}]=0.\end{eqnarray} \tag{ 20 }$$

This demonstrates the conservation of momentum for a free electron.

It is with the angular momentum that the surprises begin, as the orbital angular momentum for a free electron is not conserved. We can define an orbital angular momentum operator precisely as in non-relativistic quantum theory

$$\begin{eqnarray}&&{\bf L}={\bf x}\times {\bf p}.\end{eqnarray} \tag{ 21 }$$

This operator does not commute with the Dirac Hamiltonian and so satisfies the non-trivial equation of motion

$$\begin{eqnarray}&&{\bf \dot{L}}=\frac{i}{\hbar }[H,{\bf L}]={\boldsymbol{ \alpha }} \times {\bf p}.\end{eqnarray} \tag{ 22 }$$

In contrast with the non-relativistic theory, the orbital angular momentum of a free electron, or of one confined a by a rotationally symmetric potential, is not a constant of the motion. Diracʼs resolution of this paradoxical situation was to add to the orbital angular momentum a spin angular momentum

$$\begin{eqnarray}{\boldsymbol{ \Sigma }} =\frac{\hbar }{2}\left( \begin{array}{ccccccccccccccc} {\boldsymbol{ \sigma }} & 0 \\ 0 & {\boldsymbol{ \sigma }} \\ \end{array} \right).\end{eqnarray} \tag{ 23 }$$

Like the orbital angular momentum, this is not a constant of the motion, but satisfies the equation of motion

$$\begin{eqnarray}&&{\boldsymbol{\dot{ \Sigma }}} =-{\boldsymbol{ \alpha }} \times {\bf p}.\end{eqnarray} \tag{ 24 }$$

In relativistic quantum theory the spin and orbital angular momenta are naturally coupled together and it is only the total angular momentum, given by the sum of these, that is conserved. The total angular momentum operator

$$\begin{eqnarray}&&{\bf J}={\bf L}+{\boldsymbol{ \Sigma }} \end{eqnarray} \tag{ 25 }$$

is a constant of the motion. It is this quantity, rather than the orbital or spin angular momenta, that appears in the relativistic expressions for atomic spectra [50, 51, 54]. We should note that the separation implicit in (25) is not, itself, Lorentz covariant.

It is apparent that the velocity operator for a Dirac electron is very different from that in non-relativistic quantum theory. This feature, together with the non-conservation of spin and orbital angular momentum and a desire to understand the non-relativistic limit motivated the discovery of a unitary transformation that diagonalizes the Dirac Hamiltonian [55]. The required unitary operator is generated by an Hermitian and time-independent operator, S, and has the form [51, 55]

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm I}S}}={\rm exp} \left( \theta \beta {\boldsymbol{ \alpha }} \cdot {\bf p} \right)={\rm cos} (p\theta )+\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p}{\rm sin} (p\theta ),\end{eqnarray} \tag{ 26 }$$

where $p=|{\bf p}|$ is the magnitude of the momentum and θ is a parameter to be determined. It is not difficult to find that the required value of θ is

$$\begin{eqnarray}&&\theta =\frac{1}{2p}{{{\rm tan} }^{-1}}\left( \frac{p}{m} \right).\end{eqnarray} \tag{ 27 }$$

The Foldy–Wouthuysen transformation produces a new state and Dirac equation of the form

$$\begin{eqnarray}\begin{array}{rcl} \psi ^{\prime} & = & {{{\rm e}}^{{\rm I}S}}\psi \\ i\hbar \frac{\partial }{\partial t}\psi ^{\prime} & = & H^{\prime} \psi ^{\prime} , \\ \end{array}\end{eqnarray} \tag{ 28 }$$

where $H^{\prime}$ is the transformed (diagonal) Hamiltonian

$$\begin{eqnarray}&&H^{\prime} ={{{\rm e}}^{{\rm I}S}}H{{{\rm e}}^{-{\rm I}S}}=\beta \sqrt{{{m}^{2}}+{{p}^{2}}}.\end{eqnarray} \tag{ 29 }$$

The operator β has eigenvalues $\pm 1$, with the positive eigenvalue corresponding to the upper two components of ψ and the negative to the lower two. Clearly, the transformation has succeeded in separating the positive energy states from those with negative energy. The eigenvalues of $H^{\prime}$, moreover, are $\pm \sqrt{{{m}^{2}}+{{p}^{2}}}$ in agreement with the relationship between energy, mass and momentum, ${{E}^{2}}={{m}^{2}}+{{p}^{2}}$, familiar from relativistic kinematics.

Our transformed Dirac equation has the form

$$\begin{eqnarray}&&i\hbar \frac{\partial }{\partial t}\psi ^{\prime} =\beta \sqrt{{{m}^{2}}+{{p}^{2}}}\psi ^{\prime} .\end{eqnarray} \tag{ 30 }$$

If we consider only low energies in the non-relativistic regime then this, on making a binomial expansion, reduces to

$$\begin{eqnarray}&&i\hbar \frac{\partial }{\partial t}\psi ^{\prime} \approx \left( m-\frac{{{\hbar }^{2}}}{2\;m}{{\nabla }^{2}} \right)\beta \psi ^{\prime} .\end{eqnarray} \tag{ 31 }$$

If we further restrict our attention only to the two-component positive energy part of the state and shift our zero of energy, so as to remove the constant m term, then we recover the familiar Schrödinger equation for a free particle.

The operators representing the mechanical and other properties of the electron are also transformed by the Foldy–Wouthuysen transformation. In particular, the position operator ${\bf x}$ is transformed to [55]

$$\begin{eqnarray}\begin{array}{rcl} {\bf x}^{\prime} & = & {{{\rm e}}^{{\rm I}S}}{\bf x}{{{\rm e}}^{-{\rm I}S}} \\ {} & = & {\bf x}-\frac{i\hbar \beta {\boldsymbol{ \alpha }} }{2{{E}_{p}}}-\frac{i{\bf p}\hbar \beta {\boldsymbol{ \alpha }} \cdot {\bf p}-2p{\boldsymbol{ \Sigma }} \times {\bf p}}{2{{E}_{p}}({{E}_{p}}+m)p}, \\ \end{array}\end{eqnarray} \tag{ 32 }$$

where E_p is the operator $\sqrt{{{m}^{2}}+{{p}^{2}}}$. The momentum operator commutes with S and so is unchanged by the transformation. The orbital and spin angular momenta are also changed by the Foldy–Wouthuysen transformation, with the orbital angular momentum operator becoming

$$\begin{eqnarray}&&{\bf L}^{\prime} ={\bf x}^{\prime} \times {\bf p}.\end{eqnarray} \tag{ 33 }$$

The total angular momentum commutes with S and so is unchanged by the transformation. We can readily use this fact to obtain the transformed spin angular momentum operator. It is important to realize that using the transformed Hamiltonian, states and observables is entirely equivalent to working with the untransformed quantities and the original Dirac equation.

In the non-relativistic limit the position and momentum operators are ${\bf x}$ and ${\bf p}$, as they are in the Dirac theory. With this non-relativistic limit in mind, it is reasonable to ask what role these operators play in the transformed system. We note that in the transformed picture, the original position operator satisfies the equation of motion

$$\begin{eqnarray}&&{\bf \dot{x}}=\frac{i}{\hbar }[H^{\prime} ,{\bf x}]=\frac{{\bf p}(\beta m+{\boldsymbol{ \alpha }} \cdot {\bf p})}{{{m}^{2}}+{{p}^{2}}}.\end{eqnarray} \tag{ 34 }$$

For a positive energy state this gives the velocity operator we might have expected: ${\bf p}/{{E}_{p}}$. For a negative energy state the velocity is the same apart from an overall minus sign. The, Zitterbewegung or rapid oscillations, associated with an instantaneous velocity equal to that of light, has disappeared. For this reason, Foldy and Wouthuysen refer to the operator ${\bf x}$ in their transformed system as the, mean-position operator. Working with this mean-position operator, ${\bf x}$, and the transformed Dirac equation (30) gives the familiar Schrödinger evolution in the non-relativistic limit. The other original operators ${\bf L}$ and Σ are similarly the, mean-orbital and mean-spin angular momenta. It is important to note that it is these mean quantities that correspond to the those familiar from non-relativistic quantum theory [55].

The mean quantities can also be used with the original Dirac equation by simply inverting the Foldy–Wouthuysen transformation. This means, in particlar, that the mean position, orbital and spin angular momentum operators appropriate for use with the Dirac equation are, respectively,

$$\begin{eqnarray}\begin{array}{rcl} {{\bf \bar{x}}} & = & {{{\rm e}}^{-{\rm I}S}}{\bf x}{{{\rm e}}^{{\rm I}S}}, \\ {{\bf \bar{L}}} & = & {{{\rm e}}^{-{\rm I}S}}{\bf L}{{{\rm e}}^{{\rm I}S}}, \\ {{\boldsymbol{\bar{ \Sigma }}} } & = & {{{\rm e}}^{-{\rm I}S}}{\boldsymbol{ \Sigma }} {{{\rm e}}^{{\rm I}S}}, \\ \end{array}\end{eqnarray} \tag{ 35 }$$

echoing the quantities introduced by Pryce [56, 57]. The mean orbital and spin angular momenta, unlike their unaveraged counterparts, commute with the Hamiltonian and so are true constants of the motion. In the Foldy–Wouthuysen picture this means that

$$\begin{eqnarray}&&[H^{\prime} ,{\bf L}]=0=[H^{\prime} ,{\boldsymbol{ \Sigma }} ]\end{eqnarray} \tag{ 36 }$$

or, for the original Dirac Hamiltonian

$$\begin{eqnarray}&&[H,{\bf \bar{L}}]=0=[H,{\boldsymbol{\bar{ \Sigma }}} ].\end{eqnarray} \tag{ 37 }$$

The derivation of Poyntingʼs theorem as the continuity equation for the optical Dirac equation led us to identify the density of energy and the flux of energy, or the momentum density, as

$$\begin{eqnarray}\begin{array}{rcl} w & = & {{\psi }^{\dagger }}\psi \\ {\bf S} & = & {{\psi }^{\dagger }}{\boldsymbol{ \alpha }} \psi . \\ \end{array}\end{eqnarray} \tag{ 42 }$$

To these we can add the density of the optical angular momentum, which we write as

$$\begin{eqnarray}&&{\bf j}={{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi ={\bf r}\times ({\bf E}\times {\bf B}).\end{eqnarray} \tag{ 43 }$$

We trust that no confusion will arise if we stick with established convention in which both the optical angular-momentum density and the electron probability current are represented by ${\bf j}$. Each of these mechanical quantities is globally conserved, of course, and we should confirm this fact using the properties of our optical Dirac equation. The global conservation of energy follows directly from the local conservation law (41)

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}\int {{{\rm d}}^{3}}x{{\psi }^{\dagger }}\psi & = & -\int {{{\rm d}}^{3}}x\nabla \cdot \left( {{\psi }^{\dagger }}{\boldsymbol{ \alpha }} \psi \right) \\ {} & = & -\int {{\psi }^{\dagger }}{\boldsymbol{ \alpha }} \psi \cdot {\rm d}{\bf s} \\ {} & = & 0, \\ \end{array}\end{eqnarray} \tag{ 44 }$$

where we have used Gaussʼs theorem and the physical requirement that the fields tend to zero at infinity.

For the remaining properties we can proceed in the same way, using the Dirac equation for ψ and for ${{\psi }^{\dagger }}$. Alternatively, we can exploit further the analogy between the (classical) optical Dirac equation and the (quantum) Dirac equation for the electron, by introducing the 'Hamiltonian' operator

$$\begin{eqnarray}&&H={\boldsymbol{ \alpha }} \cdot {\bf p}=-i\hbar {\boldsymbol{ \alpha }} \cdot \nabla \end{eqnarray} \tag{ 45 }$$

and use Heisenbergʼs equation of motion. Let us emphasize that this is an entirely equivalent procedure to working with the optical Dirac equation itself and that, despite the appearance of ℏ, the description is essentially classical. The optical momentum, for example is, expressed in terms of the operator α and the momentum, density at position ${\bf r}$ therefore corresponds to the operator ${\boldsymbol{ \alpha }} \delta ({\bf r}-{\bf x})$, where ${\bf x}$ is the position operator, as before [66–68]. Both approaches lead to the same result. For the electromagnetic momentum density we find

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial }{\partial t}\left( {{\psi }^{\dagger }}{{\alpha }_{i}}\psi \right)=-\left( {{\nabla }_{j}}{{\psi }^{\dagger }} \right){{\alpha }_{j}}{{\alpha }_{i}}\psi -{{\psi }^{\dagger }}{{\alpha }_{i}}{{\alpha }_{j}}\left( {{\nabla }_{j}}\psi \right), \\ \Rightarrow \frac{\partial }{\partial t}{{S}_{i}}+{{\nabla }_{j}}{{T}_{ji}}=0, \\ \end{array}\end{eqnarray} \tag{ 46 }$$

where we have introduced the summation convention in which a repeated index implies a summation over the three cartesian coordinates and T_ji is the familiar momentum flux density [36]

$$\begin{eqnarray}&&{{T}_{ji}}=\frac{1}{2}{{\delta }_{ij}}({{E}^{2}}+{{B}^{2}})-{{E}_{i}}{{E}_{j}}-{{B}_{i}}{{B}_{j}}.\end{eqnarray} \tag{ 47 }$$

In deriving the momentum continuity equation (46) we use the commutation properties of the Kemmer matrices, given in appendix A and the transverse nature of ${\bf E}$ and ${\bf B}$ which mean that ${\bf p}\cdot {\bf E}=0={\bf p}\cdot {\bf B}$. If we follow the same procedure for the operator ${\bf x}\times {\boldsymbol{ \alpha }}$ then we find the familiar local conservation for the electromagnetic angular momentum

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}\left( {{\psi }^{\dagger }}{{({\bf x}\times {\boldsymbol{ \alpha }} )}_{i}}\psi \right) & = & -\left( {{\nabla }_{j}}{{\psi }^{\dagger }} \right){{\alpha }_{j}}{{\psi }^{\dagger }}{{({\bf x}\times {\boldsymbol{ \alpha }} )}_{i}}\psi -{{\psi }^{\dagger }}{{({\bf x}\times {\boldsymbol{ \alpha }} )}_{i}}{{\alpha }_{j}}\left( {{\nabla }_{j}}\psi \right), \\ \Rightarrow \frac{\partial }{\partial t}{{j}_{i}}+{{\nabla }_{j}}{{M}_{ji}} & = & 0, \\ \end{array}\end{eqnarray} \tag{ 48 }$$

where ${\bf j}$ is given by (43) and M_ji is the angular-momentum flux density [36, 69]

$$\begin{eqnarray}&&{{M}_{ji}}={{\varepsilon }_{ilk}}{{x}_{l}}{{T}_{jk}}.\end{eqnarray} \tag{ 49 }$$

The expressions obtained in the preceding subsection are exact (within the Maxwell theory of the free field) and contain both rapidly and slowly oscillating contributions. For a monochromatic field of angular frequency ω, for example, the densities and flux densities of the energy, momentum and angular momentum include constant terms and terms oscillating at frequency $2\omega$. In the majority of situations in optics, this frequency is too high to be observed in the interaction with matter and it is appropriate to remove the rapidly oscillating contributions so as to arrive at slowly-varying or mean properties. This process is known by different terms in different parts of optics: we may think of it as averaging over many optical cycles, the slowly-varying amplitude approximation or, within quantum optics, as making the rotating-wave approximation in the interaction with a detector. The situation is reminiscent and indeed precisely analogous to the, Zitterbewegung apparent in the physical quantities calculated from the Dirac equation for the electron. It is natural to construct the appropriate cycle-averaged optical properties by adopting the same procedure as for the electron: we apply a Foldy–Wouthuysen transformation to diagonalize the Hamiltonian for our optical Dirac equation and then seek suitable operators from which to obtain mean mechanical properties.

The optical Foldy–Wouthuysen transformation we seek is the natural analogue of that for the electron [47] with the required unitary operator having the form

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm I}S}}={\rm exp} \left( \theta \beta {\boldsymbol{ \alpha }} \cdot {\bf p} \right),\end{eqnarray} \tag{ 50 }$$

where θ is a function, to be determined, of the, operator $p=|{\bf p}|$. We can expand the unitary operator (50) in a form analogous to that given for the electron transformation in (26), although the form is complicated, somewhat, by the algebra of the Kemmer, as opposed to Dirac matrices. We find the form

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm I}S}}=\frac{{{\mathcal{P}}^{2}}}{{{p}^{2}}}+\left( 1-\frac{{{\mathcal{P}}^{2}}}{{{p}^{2}}} \right){\rm cos} (p\theta )+\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p}{\rm sin} (p\theta ),\end{eqnarray} \tag{ 51 }$$

where we define the operator

$$\begin{eqnarray}&&{{\mathcal{P}}^{2}}={{p}^{2}}-{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}},\end{eqnarray} \tag{ 52 }$$

the properties of which are given in appendix A.

The desired diagonalized Hamiltonian is obtained by setting $\theta =\pi /(4p)$, as shown in appendix B

$$\begin{eqnarray}\begin{array}{rcl} H^{\prime} & = & {{{\rm e}}^{{\rm I}S}}H{{{\rm e}}^{-{\rm I}S}} \\ {} & = & \beta \left( p-\frac{{{\mathcal{P}}^{2}}}{p} \right). \\ \end{array}\end{eqnarray} \tag{ 53 }$$

It follows that the eigenvalues of the Hamiltonian are $\pm p$, corresponding to the wavenumber or angular frequency $\omega =p/\hbar$, as ${{\mathcal{P}}^{2}}\psi =0$.

To complete the transformation we need to find the transformed spinor

$$\begin{eqnarray}\begin{array}{rcl} \psi ^{\prime} & = & {{{\rm e}}^{{\rm I}S}}\psi ={\rm exp} \left( \frac{\pi }{4p}\beta {\boldsymbol{ \alpha }} \cdot {\bf p} \right)\psi \\ {} & = & \frac{1}{\sqrt{2}}\left( 1+\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p} \right)\psi . \\ \end{array}\end{eqnarray} \tag{ 54 }$$

We can write this in terms of our electric and magnetic fields if we introduce complex electric and magnetic fields corresponding to the positive, ${{{\bf E}}^{(+)}},{{{\bf B}}^{(+)}}$, and negative, ${{{\bf E}}^{(-)}},{{{\bf B}}^{(-)}}$, frequency parts of the field

$$\begin{eqnarray}\begin{array}{rcl} {\bf E} & = & {{{\bf E}}^{(+)}}+{{{\bf E}}^{(-)}} \\ {\bf B} & = & {{{\bf B}}^{(+)}}+{{{\bf B}}^{(-)}}, \\ \end{array}\end{eqnarray} \tag{ 55 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \hbar \nabla \times {{{\bf E}}^{(\pm )}} & = & \pm ip{{{\bf B}}^{(\pm )}} \\ \hbar \nabla \times {{{\bf B}}^{(\pm )}} & = & \mp ip{{{\bf E}}^{(\pm )}}. \\ \end{array}\end{eqnarray} \tag{ 56 }$$

These positive and negative frequency parts of the field are those that arise naturally in the theory of optical coherence and are associated, in the quantum theory of light, with the photon annihilation and creation operators respectively [70–72]. The transformed spinor naturally separates these components with positive frequency electric components in the top three entries and negative frequency components in the lower three

$$\begin{eqnarray}&&\psi ^{\prime} =\left( \begin{array}{ccccccccccccccc} {{{\bf E}}^{(+)}} \\ i{{{\bf B}}^{(-)}} \\ \end{array} \right).\end{eqnarray} \tag{ 57 }$$

It is useful to note that the action of ${\boldsymbol{ \alpha }} \times {\bf p}/p$ on our transformed spinor gives the minus complex conjugate

$$\begin{eqnarray}&&\frac{{\boldsymbol{ \alpha }} \cdot {\bf p}}{p}\psi ^{\prime} =-\psi {{^{\prime} }^{*}},\end{eqnarray} \tag{ 58 }$$

and similary

$$\begin{eqnarray}&&\psi {{^{\prime} }^{\dagger }}\frac{{\boldsymbol{ \alpha }} \cdot {\bf p}}{p}=-\psi {{^{\prime} }^{*\dagger }}=-\psi {{^{\prime} }^{{\rm T}}}.\end{eqnarray} \tag{ 59 }$$

We emphasize that the combination of this transformed spinor and the transformed Hamiltonian, $H^{\prime}$, is exact and fully equivalent to our untransformed quantities.

We conclude this subsection by noting that the paraxial approximation, a mainstay of laser physics [73, 74], is readily derived from our transformed optical Dirac equation. It follows from (53) and the transverse nature of the fields that our transformed Dirac equation has the form

$$\begin{eqnarray}&&i\hbar \frac{\partial }{\partial t}\psi ^{\prime} =\beta p\psi ^{\prime} .\end{eqnarray} \tag{ 60 }$$

There is no mass term, as in the corresponding equation for the electron (30), but we can arrive at a similar Schrödinger-like equation if one component of the momentum, which we take to be the z-component, is far larger than the others, so that ${{p}_{z}}\approx p=E$. In this case we write

$$\begin{eqnarray}&&p=\sqrt{p_{z}^{2}+p_{x}^{2}+p_{y}^{2}}\approx {{p}_{z}}+\frac{p_{x}^{2}+p_{y}^{2}}{2p}.\end{eqnarray} \tag{ 61 }$$

If we make the ansatz

$$\begin{eqnarray}&&\psi ^{\prime} ={{{\rm e}}^{{\rm I}p(\beta z-t)/\hbar }}\psi ^{\prime\prime} \end{eqnarray} \tag{ 62 }$$

and multiply the equation by β, then our transformed Dirac equation becomes

$$\begin{eqnarray}&&i\left( \frac{\partial }{\partial z}+\beta \frac{\partial }{\partial t} \right)\psi ^{\prime\prime} =-\frac{\hbar }{2p}\left( \frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}} \right)\psi ^{\prime\prime} ,\end{eqnarray} \tag{ 63 }$$

so that cartesian components of the positive and negative frequency parts of both the electric and magnetic fields obey, in this approximation, a suitable paraxial wave equation [73, 74].

We can use our transformed spinor, $\psi ^{\prime}$, to arrive at mean mechanical properties by following the procedure Foldy and Wouthuysen applied to the electron Dirac equation [55]. We transform the operators corresponding to the mechanical properties of interest and drop terms that are of odd order in the α, as these introduce cross terms between the positive and negative frequency parts of the spinor. This serves to remove the rapidly oscillating contributions in the same way as the corresponding procedure for the electron removes the, Zitterbewegung.

Let us begin by considering the energy density. We have transformed the spinor and our form for the energy density must be unchanged if we also transform the operator the expectation value of which is the energy density. It is straightforward to see that the operator for the energy density at position ${\bf r}$ must be $\delta ({\bf r}-{\bf x})$, where ${\bf x}$ is the position operator, as

$$\begin{eqnarray}&&w({\bf r})=\int {{{\rm d}}^{3}}x{{\psi }^{\dagger }}\delta ({\bf r}-{\bf x})\psi .\end{eqnarray} \tag{ 64 }$$

Naturally, we obtain exactly the same result if we work with our transformed spinor $\psi ^{\prime}$ and the transformed energy density operator

$$\begin{eqnarray}\begin{array}{rcl} w({\bf r}) & = & \int {{{\rm d}}^{3}}x\psi {{^{\prime} }^{\dagger }}{{{\rm e}}^{{\rm I}S}}\delta ({\bf r}-{\bf x}){{{\rm e}}^{-{\rm I}S}}\psi ^{\prime} \\ {} & = & \frac{1}{2}\int {{{\rm d}}^{3}}x\psi {{^{\prime} }^{\dagger }}\left( 1+\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p} \right)\delta ({\bf r}-{\bf x})\left( 1-\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p} \right)\psi ^{\prime} , \\ \end{array}\end{eqnarray} \tag{ 65 }$$

where we have used the fact that both $\psi ^{\prime}$ and $\psi {{^{\prime} }^{\dagger }}$ contain only transverse fields so that ${{\mathcal{P}}^{2}}\psi ^{\prime} =0$ and $\psi {{^{\prime} }^{\dagger }}{{\mathcal{P}}^{2}}=0$. We obtain the mean or cycle-averaged energy density by keeping only those terms that are of even order in α

$$\begin{eqnarray}&&\begin{array}{l} \bar{w}({\bf r})=\frac{1}{2}\int {{{\rm d}}^{3}}x\psi {{^{\prime} }^{\dagger }}\delta ({\bf r}-{\bf x})\psi ^{\prime} -\frac{1}{2}\int {{{\rm d}}^{3}}x\psi {{^{\prime} }^{\dagger }}\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p}\delta ({\bf r}-{\bf x})\frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p}\psi ^{\prime} \\ =\frac{1}{2}\left( \psi {{^{\prime} }^{\dagger }}({\bf r})\psi ^{\prime} ({\bf r})+\psi {{^{\prime} }^{*\dagger }}({\bf r})\psi {{^{\prime} }^{*}}({\bf r}) \right) \\ =\psi {{^{\prime} }^{\dagger }}({\bf r})\psi ^{\prime} ({\bf r}) \\ ={{{\bf E}}^{(-)}}\cdot {{{\bf E}}^{(+)}}+{{{\bf B}}^{(-)}}\cdot {{{\bf B}}^{(+)}}, \\ \end{array}\end{eqnarray} \tag{ 66 }$$

which is the required average. In reaching this result we have used the anti-commutation relation (A.1) and the complex-conjugation properties (58) and (59). This simple expression for the cycle-averaged energy density illustrates the benefit of working with the transformed spinor.

We can continue in the same way, by keeping terms of only even order in α to find the averaged forms of our other mechanical properties. For the momentum density we find

$$\begin{eqnarray}\begin{array}{rcl} {{\bf \bar{S}}} & = & \frac{1}{2}\int {{{\rm d}}^{3}}x\psi {{^{\prime} }^{\dagger }}\left[ \frac{\beta {\boldsymbol{ \alpha }} \cdot {\bf p}}{p},{\boldsymbol{ \alpha }} \delta ({\bf r}-{\bf x}) \right]\psi ^{\prime} \\ {} & = & -\frac{1}{2}\int {{{\rm d}}^{3}}x\left[ \psi {{^{\prime} }^{{\rm T}}}\beta {\boldsymbol{ \alpha }} \delta ({\bf r}-{\bf x})\psi ^{\prime} +\psi {{^{\prime} }^{\dagger }}\delta ({\bf r}-{\bf x}){\boldsymbol{ \alpha }} \beta \psi {{^{\prime} }^{*}} \right] \\ {} & = & {{{\bf E}}^{(+)}}\times {{{\bf B}}^{(-)}}+{{{\bf E}}^{(-)}}\times {{{\bf B}}^{(+)}}, \\ \end{array}\end{eqnarray} \tag{ 67 }$$

which is the required cycle-averaged form of Poyntingʼs vector. In reaching this expression, we have made use of the readily proved identity

$$\begin{eqnarray}&&{\bf a}\cdot {\boldsymbol{ \alpha }} \psi ^{\prime} =\left( \begin{array}{ccccccccccccccc} -{\bf a}\times {{{\bf B}}^{(-)}} \\ i{\bf a}\times {{{\bf E}}^{(+)}} \\ \end{array} \right).\end{eqnarray} \tag{ 68 }$$

The averaged momentum density (66), although evidently correct, is not in its most useful form. We can obtain a form closer to that used in optics [75] by evaluating the operator product $({\boldsymbol{ \alpha }} \times {\bf p}){\boldsymbol{ \alpha }}$ using (A.9). After a little work we find

$$\begin{eqnarray}&&\begin{array}{r} {\bf \bar{S}}=\frac{1}{2}\left[ E_{i}^{(-)}\overleftarrow{\frac{{\bf p}}{p}}E_{i}^{(+)}-B_{i}^{(+)}\overrightarrow{\frac{{\bf p}}{p}}B_{i}^{(-)}+E_{i}^{(-)}\overrightarrow{\frac{{\bf p}}{p}}E_{i}^{(+)}-B_{i}^{(+)}\overleftarrow{\frac{{\bf p}}{p}}B_{i}^{(-)} \right. \\ \qquad \,-{{{\bf E}}^{(-)}}\left( \overleftarrow{\frac{{\bf p}}{p}}\cdot {{{\bf E}}^{(+)}} \right)+{{{\bf B}}^{(+)}}\left( \overleftarrow{\frac{{\bf p}}{p}}\cdot {{{\bf B}}^{(-)}} \right) \\ \quad \ \ \;\left. -\left( {{{\bf E}}^{(-)}}\cdot \overrightarrow{\frac{{\bf p}}{p}} \right){{{\bf E}}^{(+)}}+\left( {{{\bf B}}^{(+)}}\cdot \overrightarrow{\frac{{\bf p}}{p}} \right){{{\bf B}}^{(-)}} \right]. \\ \end{array}\end{eqnarray} \tag{ 69 }$$

We recall that ${\bf \vec{p}}=-i\hbar \nabla$ and ${\bf \overset{\scriptscriptstyle\leftarrow}{p}}=i\hbar \nabla$ and so write this in the simpler form

$$\begin{eqnarray}&&\begin{array}{r} {\bf \bar{S}}={\rm Im}\left[ E_{i}^{(-)}\left( \frac{\nabla }{\omega }E_{i}^{(+)} \right)+B_{i}^{(-)}\left( \frac{\nabla }{\omega }B_{i}^{(+)} \right) \right] \\ \qquad +{{\nabla }_{i}}{\rm Im}\left[ E_{i}^{(+)}\left( \frac{{{{\bf E}}^{(-)}}}{\omega } \right)+B_{i}^{(+)}\left( \frac{{{{\bf B}}^{(-)}}}{\omega } \right) \right]. \\ \end{array}\end{eqnarray} \tag{ 70 }$$

We emphasize that here ω is an, operator and that no restriction to monochromatic or near-monochromatic fields is necessary. The action of the operator $1/\omega$ simply divides each temporal Fourier component of the field by its associated angular frequency. The second line in (70) is a divergence and therefore does not contribute to the total momentum density. It should be noted that this operator introduces some spatial averring in that, for example, ${{{\bf B}}^{(-)}}({\bf r})/\omega$ is not simply related to ${\bf B}$ at position ${\bf r}$ alone. The density of a mechanical property is only defined up to such a divergence and so we are at liberty to neglect this term. If this statement needs justification, we recall that conserved quantities, like momentum, may be obtained from symmetries by means of Noetherʼs theorem and that the local resulting local conservation law only, defines the density up to a divergence [76–79]² . We note, however, that retaining such divergence terms allows us to express the momentum density in a variety of suggestive forms [80], although, because of the non-uniqueness described above, it is difficult to assign to these any particular physical significance. If we drop the divergence term then we are left with the simple form

$$\begin{eqnarray}\begin{array}{rcl} {{\bf \bar{S}}} & = & {\rm Im}\left[ E_{i}^{(-)}\frac{\nabla }{\omega }E_{i}^{(+)}+B_{i}^{(-)}\frac{\nabla }{\omega }B_{i}^{(+)} \right] \\ {} & = & {\rm Im}\left( \psi {{^{\prime} }^{\dagger }}\frac{\nabla }{\omega }\psi ^{\prime} \right), \\ \end{array}\end{eqnarray} \tag{ 71 }$$

which is reminiscent of the momentum density in non-relativistic quantum theory [81], although the presence of the frequency operator, ω is an important difference.

The problem of separating optical angular momentum into spin and orbital parts has an interesting history. It has been suggested, on theoretical grounds, that such a separation is not physically meaningful [86–91], and such a conclusion is, at first sight, natural as these quantities are not separately conserved for the Dirac electron [48]. This position is somewhat at odds, however, with experimental work which suggested, at least in some circumstances, that precisely such a separation is meaningful [92, 93]. The resolution of this apparent contradiction is that the separation of the total angular momentum into spin and orbital parts is indeed a physical one, but that neither the spin nor the orbital components is, by itself, an angular momentum [94, 95]. Both the spin and orbital parts are the generators of rotations that are constrained by the requirement that the rotated fields must be transverse [38, 39].

We have seen that the total angular momentum density is expressed, in the formalism of the optical Dirac equation in terms of the alpha matrices in the form

$$\begin{eqnarray}&&{\bf j}={{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi ={\bf r}\times ({\bf E}\times {\bf B}).\end{eqnarray} \tag{ 72 }$$

The spin and orbital parts of this total angular momentum are usually extracted by writing ${\bf B}$ in terms of the vector potential [94, 95] or both ${\bf B}$ in terms of the vector potential and ${\bf E}$ in terms of a second vector potential [38, 39] and using integration by parts. Here we make this division directly form the properties of the optical Dirac equation.

At first sight it is far from obvious how to break up ${{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi$ into constant spin and orbital parts. We achieve this by making use of a property of the Kemmer matrices (A.11)

$$\begin{eqnarray}\begin{array}{rcl} {{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}} & = & {{p}^{2}}-{{\mathcal{P}}^{2}} \\ \Rightarrow \frac{{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}}}{{{p}^{2}}} & = & {\rm I}-\frac{{{\mathcal{P}}^{2}}}{{{p}^{2}}}. \\ \end{array}\end{eqnarray} \tag{ 73 }$$

We recall that ${{\mathcal{P}}^{2}}\psi =0$, by virtue of the transverse nature of the ${\bf E}$ and ${\bf B}$ fields, and hence for any physically allowed fields we have

$$\begin{eqnarray}&&\frac{{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}}}{{{p}^{2}}}\psi =\psi ,\qquad {{\psi }^{\dagger }}\frac{{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}}}{{{p}^{2}}}={{\psi }^{\dagger }}.\end{eqnarray} \tag{ 74 }$$

It follows that we are free to rewrite the total angular momentum density in the form

$$\begin{eqnarray}\begin{array}{rcl} {{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi & = & {{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\frac{{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}}}{{{p}^{2}}}\psi \\ {} & = & {{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} ){\boldsymbol{ \alpha }} \cdot {\bf p}\frac{{\boldsymbol{ \alpha }} \cdot {\bf p}}{{{p}^{2}}}\psi . \\ \end{array}\end{eqnarray} \tag{ 75 }$$

We proceed by evaluating the effect on ψ of the operator ${\boldsymbol{ \alpha }} \cdot {\bf p}/{{p}^{2}}$ In order to do so we introduce the transverse (divergenceless) and, gauge-invariant parts of two vector potentials, ${{{\bf A}}^{\bot }}$ and ${{{\bf C}}^{\bot }}$, related to our (transverse) electric and magnetic fields by [9, 96]

$$\begin{eqnarray}\begin{array}{rcl} {\bf E} & = & -\frac{\partial }{\partial t}{{{\bf A}}^{\bot }}=-\nabla \times {{{\bf C}}^{\bot }} \\ {\bf B} & = & -\frac{\partial }{\partial t}{{{\bf C}}^{\bot }}=\nabla \times {{{\bf A}}^{\bot }}. \\ \end{array}\end{eqnarray} \tag{ 76 }$$

Hence with the aid of elementary vector identities we find

$$\begin{eqnarray}&&\phi =\frac{{\boldsymbol{ \alpha }} \cdot {\bf p}}{{{p}^{2}}}\psi =\frac{1}{\hbar \sqrt{2}}\left( \begin{array}{ccccccccccccccc} i{{{\bf A}}^{\bot }} \\ -{{{\bf C}}^{\bot }} \\ \end{array} \right).\end{eqnarray} \tag{ 77 }$$

We can complete our calculation of the total angular momentum density by acting with $({\bf r}\times {\boldsymbol{ \alpha }} ){\boldsymbol{ \alpha }} \cdot {\bf p}$ to the, left on ${{\psi }^{\dagger }}$. The result is

$$\begin{eqnarray}&&{{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi ={{\psi }^{\dagger }}({\bf r}\times {\bf p})\phi +{{\psi }^{\dagger }}{\boldsymbol{ \Sigma }} \phi .\end{eqnarray} \tag{ 78 }$$

It is natural to identify these with the orbital, l, and spin, s parts of the angular momentum density and if we write these in terms of the fields and the potentials we find

$$\begin{eqnarray}\begin{array}{rcl} {\bf l} & = & {{\psi }^{\dagger }}({\bf r}\times {\bf p})\phi =\frac{1}{2}\left( {{E}_{i}}{\bf r}\times \nabla A_{i}^{\bot }+{{B}_{i}}{\bf r}\times \nabla C_{i}^{\bot } \right) \\ {\bf s} & = & {{\psi }^{\dagger }}{\boldsymbol{ \Sigma }} \phi =\frac{1}{2}\left( {\bf E}\times {{{\bf A}}^{\bot }}+{\bf B}\times {{{\bf C}}^{\bot }} \right), \\ \end{array}\end{eqnarray} \tag{ 79 }$$

which we recognize as the required Heaviside–Larmor symmetric forms [38, 39, 78]. Had we chosen to act with $({\bf r}\times {\boldsymbol{ \alpha }} ){\boldsymbol{ \alpha }} \cdot {\bf p}$ to the right we would have recovered ${\bf r}\times ({\bf E}\times {\bf B})$ and it follows that this expression and ${\bf l}+{\bf s}$, given in (79), are, identical within the formalism of the optical Dirac equation. We note that, as for the electron, this separation into spin and orbital parts is not Lorentz-covariant.

We do not wish to give the impression that the spin and orbital angular momentum densities given in (79) are unique. Had we written

$$\begin{eqnarray*}&&{{\psi }^{\dagger }}({\bf r}\times {\boldsymbol{ \alpha }} )\psi ={{\psi }^{\dagger }}\frac{{{({\boldsymbol{ \alpha }} \cdot {\bf p})}^{2}}}{{{p}^{2}}}({\bf r}\times {\boldsymbol{ \alpha }} )\psi ,\end{eqnarray*}$$

for example, then we would have obtained expressions for the densities differing only by a surface term, with the total volume-integrated spin and orbital angular momenta unchanged. The possibility of defining different densities for locally conserved quantities associated with the same total properties is familiar, of course, from Noetherʼs theorem [76–78].

In particle physics, the helicity of a single particle is defined as the component of its spin along its momentum [52, 97–100], thus it is natural to define this as the expectation value of $({\boldsymbol{ \Sigma }} \cdot {\bf p})/p$. If we apply this idea to our optical Dirac equation then we are led to a helicity density in the form

$$\begin{eqnarray}&&{{h}_{1}}={{\psi }^{\dagger }}\frac{{\boldsymbol{ \Sigma }} \cdot {\bf p}}{p}\psi ,\end{eqnarray} \tag{ 80 }$$

It is certainly, possible to define the helicity in this way and, indeed, doing so has some practical advantages [101]. There are two clear reasons, however, why this does not suit our present purposes. The first is that it gives a value that is proportional to ℏ which should not appear in any physical quantity in our essentially classical analysis. The second, more serious, problem is that this quantity does not have the dimensions of an angular momentum density. We can remove both these features if we take our inspiration from (78) to give

$$\begin{eqnarray}&&{{h}_{2}}={{\psi }^{\dagger }}\frac{{\boldsymbol{ \Sigma }} \cdot {\bf p}}{p}\phi ,\end{eqnarray} \tag{ 81 }$$

but this quantity is a purely imaginary one.

We can recover the required form of the helicity density within from optical Dirac equation by defining the helicity operator to be $({\boldsymbol{ \Sigma }} \cdot {\bf p})/{{p}^{2}}$, so that the helicity density is

$$\begin{eqnarray}\begin{array}{rcl} h & = & {{\psi }^{\dagger }}\frac{{\boldsymbol{ \Sigma }} \cdot {\bf p}}{{{p}^{2}}}\psi \\ {} & = & \frac{1}{2}\left( {{{\bf A}}^{\bot }}\cdot {\bf B}-{{{\bf C}}^{\bot }}\cdot {\bf E} \right), \\ \end{array}\end{eqnarray} \tag{ 82 }$$

which is the familiar expression [40–45, 78]. The total helicity obtained by integrating this helicity density over all space is the conserved quantity associated with the Heaviside–Larmor symmetry and, in the quantum theory, is ℏ multiplied by the difference between the number of left and right circularly polarized photons.

The difference in form for the helicity operator in the optical Dirac when compared with the analogous quantity for the electron should not come as a surprise, as this has been true for each of the mechanical properties we have considered and has its origins in the fact that ${{\psi }^{\dagger }}\psi$ has the dimensions of a probability density for the electron but an energy density for the electromagnetic field [16].