format: text

text:
# Untitled Document

##

Wavelength-Dependent Effects in Maxwell Optics

Sameen Ahmed KHAN

khan@fis.unam.mx, rohelakhan@hotmail.com

http://www.pd.infn.it/$``$khan/

http://www.imsc.ernet.in/$``$jagan/khan-cv.html

Centro de Ciencias Físicas,

Universidad Nacional Autónoma de México (UNAM)

Apartado Postal 48-3, Cuernavaca 62251, Morelos, MÉXICO

## Abstract

We present a new formalism for light beam optics starting with an exact eight-dimensional matrix representation of the Maxwell equations. The Foldy-Wouthuysen iterative diagonalization technique is employed to obtain a Hamiltonian description for a system with varying refractive index. Besides, reproducing all the traditional quasiparaxial terms, this method leads to additional contributions, which are dependent on the wavelength, in the optical Hamiltonian. This alternate prescription to obtain the aberration expansion is applied to the axially symmetric graded index fiber. This results in the wavelength-dependent modifications of the paraxial behaviour and the aberration coefficients. Furthermore it predicts a wavelength-dependent image rotation. In the low wavelength limit our formalism reproduces the Lie algebraic formalism of optics. The Foldy-Wouthuysen technique employed by us is ideally suited for the Lie algebraic approach to optics. The present study further strengthens the close analogy between the various prescription of light and charged-particle optics. All the associated machinery used in this formalism is described in the text and the accompanying appendices.

##

## 1 Introduction

The traditional scalar wave theory of optics (including aberrations to all orders) is based on the beam-optical Hamiltonian derived using the Fermat’s principle. This approach is purely geometrical and works adequately in the scalar regime. The other approach is based on the Helmholtz equation which is derived from the Maxwell equations. In this approach one takes the square-root of the Helmholtz operator followed by an expansion of the radical . This approach works to all orders and the resulting expansion is no different from the one obtained using the geometrical approach of the Fermat’s principle.

Another way of obtaining the aberration expansion is based on the algebraic similarities between the Helmholtz equation and the Klein-Gordon equation. Exploiting this algebraic similarity the Helmholtz equation is linearized in a procedure very similar to the one due to Feschbach-Villars, for linearizing the Klein-Gordon equation. This brings the Helmholtz equation to a Dirac-like form and then follows the procedure of the Foldy-Wouthuysen expansion used in the Dirac electron theory. This approach, which uses the algebraic machinery of quantum mechanics, was developed recently , providing an alternative to the traditional square-root procedure. This scalar formalism gives rise to wavelength-dependent contributions modifying the aberration coefficients . The algebraic machinery of this formalism is very similar to the one used in the quantum theory of charged-particle beam optics, based on the Dirac - and the Klein-Gordon  equations respectively. The detailed account for both of these is available in . A treatment of beam optics taking into account the anomalous magnetic moment is available in -.

As for the polarization: A systematic procedure for the passage from scalar to vector wave optics to handle paraxial beam propagation problems, completely taking into account the way in which the Maxwell equations couple the spatial variation and polarization of light waves, has been formulated by analysing the basic Poincaré invariance of the system, and this procedure has been successfully used to clarify several issues in Maxwell optics -.

In all the above approaches, the beam-optics and the polarization are studied separately, using very different machineries. The derivation of the Helmholtz equation from the Maxwell equations is an approximation as one neglects the spatial and temporal derivatives of the permittivity and permeability of the medium. Any prescription based on the Helmholtz equation is bound to be an approximation, irrespective of how good it may be in certain situations. It is very natural to look for a prescription based fully on the Maxwell equations. Such a prescription is sure to provide a deeper understanding of beam-optics and polarization in a unified manner. With this as the chief motivation we construct a formalism starting with the Maxwell equations in a matrix form: a single entity containing all the four Maxwell equations.

In our approach we require an exact matrix representation of the Maxwell equations in a medium taking into account the spatial and temporal variations of the permittivity and permeability. It is necessary and sufficient to use $`8\times 8`$ matrices for such an exact representation. This representation makes use of the Riemann-Silberstein vector, which is described in Appendix-A. The derivation of the required matrix representation, and how it differs from the numerous other ones is presented in Appendix-B.

The derived matrix representation of the Maxwell equations has a very close algebraic correspondence with the Dirac equation. This enables us to apply the machinery of the Foldy-Wouthuysen expansion used in the Dirac electron theory. The Foldy-Wouthuysen transformation technique is outlined in Appendix-C. General expressions for the Hamiltonians are derived without assuming any specific form for the refractive index. These Hamiltonians are shown to contain the extra wavelength-dependent contributions which arise very naturally in our approach. In Section-IV we apply the general formalism to the specific examples: A. Medium with Constant Refractive Index. This example is essentially for illustrating some of the details of the machinery used.

The other application, B. Axially Symmetric Graded Index Medium is used to demonstrate the power of the formalism. Two points are worth mentioning, Image Rotation: Our formalism gives rise to the image rotation (proportional to the wavelength) and we have derived an explicit relationship for the angle of the image rotation. The other pertains to the aberrations: In our formalism we get all the nine aberrations permitted by the axial symmetry. The traditional approaches give six aberrations. Our formalism modifies these six aberration coefficients by wavelength-dependent contributions and also gives rise to the remaining three permitted by the axial symmetry. The existence of the nine aberrations and image rotation are well-known in axially symmetric magnetic lenses, even when treated classically. The quantum treatment of the same system leads to the wavelength-dependent modifications . The alternate procedure for the Helmholtz optics in  gives the usual six aberrations (though modified by the wavelength-dependent contributions) and does not give any image rotation. These extra aberrations and the image rotation are the exclusive outcome of the fact that the formalism is based on a treatment starting with an exact matrix representation of the Maxwell equations.

The traditional beam-optics is completely obtained from our approach in the limit wavelength, $`\text{}\lambda 0`$, which we call as the traditional limit of our formalism. This is analogous to the classical limit obtained by taking $`\mathrm{}0`$ in the quantum prescriptions. The scheme of using the Foldy-Wouthuysen machinery in this formalism is very similar to the one used in the quantum theory of charged-particle beam optics -. There too one recovers the classical prescriptions (Lie algebraic formalism of charged-particle beam optics, to be precise) in the limit $`\text{}\lambda _00`$, where $`\text{}\lambda _0=\mathrm{}/p_0`$ is the de Broglie wavelength and $`p_0`$ is the design momentum of the system under study.

In this article we focus on the Hamiltonian description of the beam optics, as is customary in the traditional prescriptions of beam optics. This also enables us to relate our formalism with the traditional prescriptions. The studies on the evolution of the fields and the polarization are very much in progress. Some of the results in  have been obtained as the lowest order approximation of the more general framework developed here. These shall be presented elsewhere .

## 2 Traditional Prescriptions

Recalling, that in the traditional scalar wave theory for treating monochromatic quasiparaxial light beam propagating along the positive $`z`$-axis, the $`z`$-evolution of the optical wave function $`\psi (𝒓)`$ is taken to obey the Schrödinger-like equation

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\psi (𝒓)=\widehat{H}\psi (𝒓),`$ (1)
where the optical Hamiltonian $`\widehat{H}`$ is formally given by the radical

$`\widehat{H}=\left(n^2(𝒓)\widehat{𝒑}_{}^2\right)^{1/2},`$ (2)
and $`n(𝒓)=n(x,y,z)`$ is the varying refractive index. In beam optics the rays are assumed to propagate almost parallel to the optic-axis, chosen to be $`z`$-axis, here. That is, $`\left|\widehat{𝒑}_{}\right|1`$. The refractive index is the order of unity. For a medium with uniform refractive index, $`n(𝒓)=n_0`$ and the Taylor expansion of the radical is

$`\left(n^2(𝒓)\widehat{𝒑}_{}^2\right)^{1/2}`$ $`=`$ $`n_0\left\{1{\displaystyle \frac{1}{n_0^2}}\widehat{𝒑}_{}^2\right\}^{1/2}`$ (3)
$`=`$ $`n_0\{1{\displaystyle \frac{1}{2n_0^2}}\widehat{𝒑}_{}^2{\displaystyle \frac{1}{8n_0^4}}\widehat{𝒑}_{}^4{\displaystyle \frac{1}{16n_0^6}}\widehat{𝒑}_{}^6`$
$`{\displaystyle \frac{5}{128n_0^8}}\widehat{𝒑}_{}^8{\displaystyle \frac{7}{256n_0^{10}}}\widehat{𝒑}_{}^{10}\mathrm{}\}.`$
In the above expansion one retains terms to any desired degree of accuracy in powers of $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)`$. In general the refractive index is not a constant and varies. The variation of the refractive index $`n(𝒓)`$, is expressed as a Taylor expansion in the spatial variables $`x`$, $`y`$ with $`z`$-dependent coefficients. To get the beam optical Hamiltonian one makes the expansion of the radical as before, and retains terms to the desired order of accuracy in $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)`$ along with all the other terms (coming from the expansion of the refractive index $`n(𝒓)`$) in the phase-space components up to the same order. In this expansion procedure the problem is partitioned into paraxial behaviour $`+`$ aberrations, order-by-order.

In relativistic quantum mechanics too, one has the problem of understanding the behaviour in terms of nonrelativistic limit $`+`$ relativistic corrections, order-by-order. In the Dirac theory of the electron this is done most conveniently through the Foldy-Wouthuysen transformation . The beam optical Hamiltonian derived, starting with the exact matrix representation of the Maxwell equations has a very close algebraic resemblance with the Dirac case, accompanied by the analogous physical interpretations. The details of this correspondence and the Foldy-Wouthuysen transformation are given in Appendix-C.

## 3 The Beam-Optical Formalism

Matrix representations of the Maxwell equations are very well-known -. However, all these representations lack an exactness or/and are given in terms of a pair of matrix equations. A treatment expressing the Maxwell equations in a single matrix equation instead of a pair of matrix equations was obtained recently -. This representation contains all the four Maxwell equations in presence of sources taking into account the spatial and temporal variations of the permittivity $`ϵ(𝒓,t)`$ and the permeability $`\mu (𝒓,t)`$.

Maxwell equations  in an inhomogeneous medium with sources are

$`\mathbf{}𝑫(𝒓,t)=\rho ,`$
$`\mathbf{}\times 𝑯(𝒓,t){\displaystyle \frac{}{t}}𝑫(𝒓,t)=𝑱,`$
$`\mathbf{}\times 𝑬(𝒓,t)+{\displaystyle \frac{}{t}}𝑩(𝒓,t)=0,`$
$`\mathbf{}𝑩(𝒓,t)=0.`$ (4)
We assume the media to be linear, that is $`𝑫=ϵ(𝒓,t)𝑬`$, and $`𝑩=\mu (𝒓,t)𝑯`$, where $`ϵ`$ is the permittivity of the medium and $`\mu `$ is the permeability of the medium. The magnitude of the velocity of light in the medium is given by $`v(𝒓,t)=\left|𝒗(𝒓,t)\right|=1/\sqrt{ϵ(𝒓,t)\mu (𝒓,t)}`$. In vacuum we have, $`ϵ_0=8.85\times 10^{12}C^2/N.m^2`$ and $`\mu _0=4\pi \times 10^7N/A^2`$. Following the notation in  we use the Riemann-Silberstein vector given by

$`𝑭^\pm (𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\sqrt{ϵ(𝒓,t)}𝑬(𝒓,t)\pm \mathrm{i}{\displaystyle \frac{1}{\sqrt{\mu (𝒓,t)}}}𝑩(𝒓,t)\right).`$ (5)
We further define,

$`\mathrm{\Psi }^\pm (𝒓,t)=\left[\begin{array}{c}F_x^\pm \pm \mathrm{i}F_y^\pm \\ F_z^\pm \\ F_z^\pm \\ F_x^\pm \pm \mathrm{i}F_y^\pm \end{array}\right],W^\pm =\left({\displaystyle \frac{1}{\sqrt{2ϵ}}}\right)\left[\begin{array}{c}J_x\pm \mathrm{i}J_y\\ J_zv\rho \\ J_z+v\rho \\ J_x\pm \mathrm{i}J_y\end{array}\right],`$ (14)
where $`W^\pm `$ are the vectors for the sources. Following the notation in  the exact matrix representation of the Maxwell equations is

$`{\displaystyle \frac{}{t}}\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]{\displaystyle \frac{\dot{v}(𝒓,t)}{2v(𝒓,t)}}\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (23)
$`+{\displaystyle \frac{\dot{h}(𝒓,t)}{2h(𝒓,t)}}\left[\begin{array}{cc}\mathrm{𝟎}& \mathrm{i}\beta \alpha _y\\ \mathrm{i}\beta \alpha _y& \mathrm{𝟎}\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (28)
$`=v(𝒓,t)\left[\begin{array}{ccc}\left\{𝑴\mathbf{}+𝚺𝒖\right\}& & \mathrm{i}\beta \left(𝚺𝒘\right)\alpha _y\\ \mathrm{i}\beta \left(𝚺^{}𝒘\right)\alpha _y& & \left\{𝑴^{}\mathbf{}+𝚺^{}𝒖\right\}\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (33)
$`\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{c}W^+\\ W^{}\end{array}\right],`$ (38)
where $`()^{}`$ denotes complex-conjugation, $`\dot{v}=\frac{v}{t}`$ and $`\dot{h}=\frac{h}{t}`$. The various matrices are

$`M_x`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{𝟎}& 1\mathrm{l}\\ 1\mathrm{l}& \mathrm{𝟎}\end{array}\right],M_y=\left[\begin{array}{cc}\mathrm{𝟎}& \mathrm{i1l}\\ \mathrm{i1l}& \mathrm{𝟎}\end{array}\right],M_z=\beta =\left[\begin{array}{cc}1\mathrm{l}& \mathrm{𝟎}\\ \mathrm{𝟎}& 1\mathrm{l}\end{array}\right],`$ (45)
$`𝚺`$ $`=`$ $`\left[\begin{array}{cc}𝝈& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝝈\end{array}\right],𝜶=\left[\begin{array}{cc}\mathrm{𝟎}& 𝝈\\ 𝝈& \mathrm{𝟎}\end{array}\right],𝑰=\left[\begin{array}{cc}1\mathrm{l}& \mathrm{𝟎}\\ \mathrm{𝟎}& 1\mathrm{l}\end{array}\right],`$ (52)
and $`1\mathrm{l}`$ is the $`2\times 2`$ unit matrix. The triplet of the Pauli matrices, $`𝝈`$ are

$`𝝈`$ $`=`$ $`\left[\sigma _x=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\sigma _y=\left[\begin{array}{cc}0\hfill & \hfill \mathrm{i}\\ \mathrm{i}\hfill & \hfill 0\end{array}\right],\sigma _z=\left[\begin{array}{cc}1\hfill & \hfill 0\\ 0\hfill & \hfill 1\end{array}\right]\right],`$ (59)
and

$`𝒖(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2v(𝒓,t)}}\mathbf{}v(𝒓,t)={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}v(𝒓,t)\right\}={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}n(𝒓,t)\right\}`$
$`𝒘(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2h(𝒓,t)}}\mathbf{}h(𝒓,t)={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}h(𝒓,t)\right\}.`$ (60)
Lastly,

$`\mathrm{Velocity}\mathrm{Function}:v(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{ϵ(𝒓,t)\mu (𝒓,t)}}}`$
$`\mathrm{Resistance}\mathrm{Function}:h(𝒓,t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mu (𝒓,t)}{ϵ(𝒓,t)}}}.`$ (61)
As we shall soon see, it is advantageous to use the above derived functions instead of the permittivity, $`ϵ(𝒓,t)`$ and the permeability, $`\mu (𝒓,t)`$. The functions, $`v(𝒓,t)`$ and $`h(𝒓,t)`$ have the dimensions of velocity and resistance respectively.

Let us consider the case without any sources ($`W^\pm =0`$). We further assume,

$`\mathrm{\Psi }^\pm (𝒓,t)=\psi ^\pm \left(𝒓\right)e^{\mathrm{i}\omega t},\omega >0,`$ (62)
with $`\dot{v}(𝒓,t)=0`$ and $`\dot{h}(𝒓,t)=0`$. Then,

$`\left[\begin{array}{cc}M_z& \mathrm{𝟎}\\ \mathrm{𝟎}& M_z\end{array}\right]{\displaystyle \frac{}{z}}\left[\begin{array}{cc}\psi ^+& \\ \psi ^{}& \end{array}\right]`$ (67)
$`=\mathrm{i}{\displaystyle \frac{\omega }{v(𝒓)}}\left[\begin{array}{cc}\psi ^+& \\ \psi ^{}& \end{array}\right]`$ (70)
$`v(𝒓)\left[\begin{array}{ccc}\left\{𝑴_{}\mathbf{}_{}+𝚺𝒖\right\}& & \mathrm{i}\beta \left(𝚺𝒘\right)\alpha _y\\ \mathrm{i}\beta \left(𝚺^{}𝒘\right)\alpha _y& & \left\{𝑴_{}^{}\mathbf{}_{}+𝚺^{}𝒖\right\}\end{array}\right]\left[\begin{array}{cc}\psi ^+& \\ \psi ^{}& \end{array}\right].`$ (75)
(76)
At this stage we introduce the process of wavization, through the familiar Schrödinger replacement

$`\mathrm{i}\text{}\lambda \mathbf{}_{}\widehat{𝒑}_{},\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}p_z,`$ (77)
where $`\text{}\lambda =\lambda /2\pi `$ is the reduced wavelength, $`c=\text{}\lambda \omega `$ and $`n(𝒓)=c/v(𝒓)`$ is the refractive index of the medium. Noting, that $`(pqqp)=\mathrm{i}\text{}\lambda `$, which is very similar to the commutation relation, $`(pqqp)=\mathrm{i}\mathrm{}`$, in quantum mechanics. In our formalism, ‘$`\lambda `$’ plays the same role which is played by the Planck constant, ‘$`\mathrm{}`$’ in quantum mechanics. The traditional beam-optics is completely obtained from our formalism in the limit $`\text{}\lambda 0`$.

Noting, that $`M_z^1=M_z=\beta `$, we multiply both sides of equation (76) by

$`\left[\begin{array}{cc}M_z& \mathrm{𝟎}\\ \mathrm{𝟎}& M_z\end{array}\right]^1=\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]`$ (82)
and $`(\mathrm{i}\text{}\lambda )`$ , then, we obtain

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\left[\begin{array}{cc}\psi ^+(𝒓_{},z)& \\ \psi ^{}(𝒓_{},z)& \end{array}\right]`$ $`=`$ $`\widehat{\mathrm{H}}_g\left[\begin{array}{cc}\psi ^+(𝒓_{},z)& \\ \psi ^{}(𝒓_{},z)& \end{array}\right].`$ (87)
This is the basic optical equation, where

$`\widehat{\mathrm{H}}_g`$ $`=`$ $`n_0\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]+\widehat{}_g+\widehat{𝒪}_g`$ (90)
$`\widehat{}_g`$ $`=`$ $`\left(n\left(𝒓\right)n_0\right)\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]\beta _g`$ (95)
$`+\left[\begin{array}{ccc}\beta \left\{𝑴_{}𝒑_{}\mathrm{i}\text{}\lambda 𝚺𝒖\right\}& & \mathrm{𝟎}\\ \mathrm{𝟎}& & \beta \left\{𝑴_{}^{}𝒑_{}\mathrm{i}\text{}\lambda 𝚺^{}𝒖\right\}\end{array}\right]`$
$`\widehat{𝒪}_g`$ $`=`$ $`\left[\begin{array}{ccc}\mathrm{𝟎}& & \text{}\lambda \left(𝚺𝒘\right)\alpha _y\\ \text{}\lambda \left(𝚺^{}𝒘\right)\alpha _y& & \mathrm{𝟎}\end{array}\right],`$ (98)
where ‘$`g`$’ stands for grand, signifying the eight dimensions and

$`\beta _g=\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right].`$ (101)
The above optical Hamiltonian is exact (as exact as the Maxwell equations in a time-independent linear media). The approximations are made only at the time of doing specific calculations. Apart from the exactness, the optical Hamiltonian is in complete algebraic correspondence with the Dirac equation with appropriate physical interpretations. The relevant point is:

$`\beta _g\widehat{}_g=\widehat{}_g\beta _g,\beta _g\widehat{𝒪}_g=\widehat{𝒪}_g\beta _g.`$ (102)
We note that the upper component ($`\psi ^+`$) is coupled to the lower component ($`\psi ^{}`$) through the logarithmic divergence of the resistance function. If this coupling function, $`𝒘=0`$, or is approximated to be zero, then the eight dimensional equations for ($`\psi ^+`$) and ($`\psi ^{}`$) get completely decoupled, leading to two independent four dimensional equations. Each of these two equations is equivalent to the other. These are the leading equations for our studies of beam-optics and polarization. In the optics context any contribution from the gradient of the resistance function can be assumed to be negligible. With this reasonable assumption we can decouple the equations and reduce the problem from eight dimensions to four dimensions. In the following sections we shall present a formalism with the approximation $`𝒘0`$. After constructing the formalism in four dimensions we shall also address the question of dealing with the contributions coming from the gradient of the resistance function. This will require the application of the Foldy-Wouthuysen transformation technique in cascade as we shall see. This justifies the usage of the two derived laboratory functions in place of permittivity and permeability respectively.

We drop the ‘<sup>+</sup>’ throughout, then the beam-optical Hamiltonian is

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\psi \left(𝒓\right)`$ $`=`$ $`\widehat{\mathrm{H}}\psi \left(𝒓\right)`$
$`\widehat{\mathrm{H}}`$ $`=`$ $`n_0\beta +\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`\left(n\left(𝒓\right)n_0\right)\beta \mathrm{i}\text{}\lambda \beta 𝚺𝒖`$
$`\widehat{𝒪}`$ $`=`$ $`\mathrm{i}\left(M_yp_xM_xp_y\right)`$ (103)
$`=`$ $`\beta \left(𝑴_{}\widehat{𝒑}_{}\right).`$
If we were to neglect the derivatives of the permittivity and permeability, we would have missed the term, $`(\mathrm{i}\text{}\lambda \beta 𝚺𝒖)`$. This is an outcome of the exact treatment.

Proceeding with our analogy with the Dirac equation: this extra term is analogous to the anomalous magnetic/electric moment term coupled to the magnetic/electric field respectively in the Dirac equation. The term we dropped (while going from the eight dimensional exact to the four dimensional almost-exact) is analogous to the anomalous magnetic/electric moment term coupled to the electric/magnetic fields respectively. However it should be borne in mind that in our exact treatment, both the terms were derived from the Maxwell equations, where as in the Dirac theory the anomalous terms are added based on experimental results (some even predating the Dirac equation) and certain arguments of invariances. In our exact treatment of the Maxwell optics, these are the only two terms one gets, where as in the Dirac equation the scheme of invariances permits addition of any number of terms! The term, $`(\mathrm{i}\text{}\lambda \beta 𝚺𝒖)`$ is related to the polarization and we shall call it as the polarization term.

One of the other similarities worth noting, relates to the square of the optical Hamiltonian.

$`\widehat{\mathrm{H}}^2`$ $`=`$ $`\left\{n^2\left(𝒓\right)\widehat{𝒑}_{}^2\right\}\text{}\lambda ^2u^2+[𝑴_{}\widehat{𝒑}_{},n\left(𝒓\right)]`$ (104)
$`+2\mathrm{i}\text{}\lambda n(𝒓)𝚺𝒖+\mathrm{i}\text{}\lambda [𝑴_{}\widehat{𝒑}_{},𝚺𝒖]`$
$`=`$ $`\left\{n\left(𝒓\right)+\mathrm{i}\text{}\lambda 𝚺𝒖\right\}^2\widehat{𝒑}_{}^2`$
$`+[𝑴_{}\widehat{𝒑}_{},\left\{n\left(𝒓\right)+\mathrm{i}\text{}\lambda 𝚺𝒖\right\}],`$
where, $`[A,B]=(ABBA)`$ is the commutator. It is to be noted that the square of the Hamiltonian in our formalism differs from the square of the Hamiltonian in the square-root approaches  and the scalar approach in . This is essentially the same type of difference which exists in the Dirac case. There too, the square of the Dirac Hamiltonian gives rise to extra pieces (such as, $`\mathrm{}q𝚺𝑩`$, the Pauli term which couples the spin to the magnetic field) which is absent in the Schrödinger and the Klein-Gordon descriptions. It is this difference in the square of the Hamiltonians which give rise to the various extra wavelength-dependent contributions in our formalism. These differences persist even in the approximation when the polarization term is neglected.

The beam optical Hamiltonian derived in (103) has a very close algebraic correspondence with the Dirac equation, accompanied by the analogous physical interpretations. This enables us to employ the machinery of the Foldy-Wouthuysen transformation technique. The details are available in Appendix-C. To the leading order, that is to order, $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)`$ the beam-optical Hamiltonian in terms of $`\widehat{}`$ and $`\widehat{𝒪}`$ is formally given by

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}|\psi `$ $`=`$ $`\widehat{}^{(2)}|\psi ,`$
$`\widehat{}^{(2)}`$ $`=`$ $`n_0\beta +\widehat{}{\displaystyle \frac{1}{2n_0}}\beta \widehat{𝒪}^2.`$ (105)
Note that, $`\widehat{𝒪}^2=\widehat{𝒑}_{}^2`$ and $`\widehat{}=\left(n\left(𝒓\right)n_0\right)\beta \mathrm{i}\text{}\lambda \beta 𝚺𝒖`$. Since, we are primarily interested in the forward propagation, we drop the $`\beta `$ from the non-matrix parts of the Hamiltonian. The matrix terms are related to the polarization. The formal Hamiltonian in (105), expressed in terms of the phase-space variables is:

$`\widehat{}^{(2)}`$ $`=`$ $`\left\{n\left(𝒓\right){\displaystyle \frac{1}{2n_0}}\widehat{𝒑}_{}^2\right\}\mathrm{i}\text{}\lambda \beta 𝚺𝒖.`$ (106)
Note that one retains terms up to quadratic in the Taylor expansion of the refractive index $`n(𝒓)`$, to be consistent with the order of $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)`$. This is the paraxial Hamiltonian which also contains an extra matrix dependent term, which we call as the polarization term. Rest of it is similar to the one obtained in the traditional approaches.

To go beyond the paraxial approximation one goes a step further in the Foldy-Wouthuysen iterative procedure. Note that, $`\widehat{𝒪}`$ is the order of $`\widehat{𝒑}_{}`$. To order $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)^2`$, the beam-optical Hamiltonian in terms of $`\widehat{}`$ and $`\widehat{𝒪}`$ is formally given by

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}|\psi `$ $`=`$ $`\widehat{}^{(4)}|\psi ,`$
$`\widehat{}^{(4)}`$ $`=`$ $`n_0\beta +\widehat{}{\displaystyle \frac{1}{2n_0}}\beta \widehat{𝒪}^2`$ (107)
$`{\displaystyle \frac{1}{8n_0^2}}[\widehat{𝒪},\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\widehat{𝒪}\right)]`$
$`+{\displaystyle \frac{1}{8n_0^3}}\beta \left\{\widehat{𝒪}^4+\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\widehat{𝒪}\right)^2\right\}.`$
Note that $`\widehat{𝒪}^4=\widehat{𝒑}_{}^4`$, and $`\frac{}{z}\widehat{𝒪}=0`$. The formal Hamiltonian in (107) when expressed in terms of the phase-space variables is

$`\widehat{}^{(4)}`$ $`=`$ $`\left\{n(𝒓){\displaystyle \frac{1}{2n_0}}\widehat{𝒑}_{}^2{\displaystyle \frac{1}{8n_0^3}}\widehat{𝒑}_{}^4\right\}`$ (108)
$`{\displaystyle \frac{1}{8n_0^2}}\{[\widehat{𝒑}_{}^2,(n(𝒓)n_0)]_+`$
$`+2(p_x(n(𝒓)n_0)p_x+p_y(n(𝒓)n_0)p_y)\}`$
$`{\displaystyle \frac{\mathrm{i}}{8n_0^2}}\left\{[p_x,[p_y,\left(n(𝒓)n_0\right)]_+][p_y,[p_x,\left(n(𝒓)n_0\right)]_+]\right\}`$
$`+{\displaystyle \frac{1}{8n_0^3}}\left\{[p_x,\left(n(𝒓)n_0\right)]_+^2+[p_y,\left(n(𝒓)n_0\right)]_+^2\right\}`$
$`+{\displaystyle \frac{\mathrm{i}}{8n_0^3}}\left\{[[p_x,\left(n(𝒓)n_0\right)]_+,[p_y,\left(n(𝒓)n_0\right)]_+]\right\}`$
$`\mathrm{}`$
where $`[A,B]_+=(AB+BA)`$ and ‘$`\mathrm{}`$’ are the contributions arising from the presence of the polarization term. Any further simplification would require information about the refractive index $`n(𝒓)`$.

Note that, the paraxial Hamiltonian (106) and the leading order aberration Hamiltonian (108) differs from the ones derived in the traditional approaches. These differences arise by the presence of the wavelength-dependent contributions which occur in two guises. One set occurs totally independent of the polarization term in the basic Hamiltonian. This set is a multiple of the unit matrix or at most the matrix $`\beta `$. The other set involves the contributions coming from the polarization term in the starting optical Hamiltonian. This gives rise to both matrix contributions and the non-matrix contributions, as the squares of the polarization matrices is unity. We shall discuss the contributions of the polarization to the beam optics elsewhere. Here, it suffices to note existence of the the wavelength-dependent contributions in two distinguishable guises, which are not present in the traditional prescriptions.

### 3.1 When w $``$ 0

In the previous sections we assumed, $`𝒘=0`$, and this enabled us to develop a formalism using $`4\times 4`$ matrices via the Foldy-Wouthuysen machinery. The Foldy-Wouthuysen transformation enables us to eliminate the odd part in the $`4\times 4`$ matrices, to any desired order of accuracy. Here too we have the identical problem, but a step higher in dimensions. So, we need to apply the Foldy-Wouthuysen to reduce the strength of the odd part in eight dimensions. This will reduce the problem from eight to four dimensions.

We start with the grand beam optical equation in (87) and proceed with the Foldy-Wouthuysen transformations as before, but with each quantity in double the number of dimensions. Symbolically this means:

$`\widehat{\mathrm{H}}\widehat{\mathrm{H}}_g,\psi \psi _g=\left[\begin{array}{cc}\psi ^+& \\ \psi ^{}& \end{array}\right],`$ (111)
$`\widehat{}\widehat{}_g,\widehat{𝒪}\widehat{𝒪}_g`$
$`n_0n_g=n_0\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right].`$ (114)
The first Foldy-Wouthuysen iteration gives

$`\widehat{}_g^{(2)}`$ $`=`$ $`n_0\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]+\widehat{}_g{\displaystyle \frac{1}{2n_0}}\beta _g\widehat{𝒪}_g^2`$ (117)
$`=`$ $`n_0\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]\beta _g+\widehat{}_g+{\displaystyle \frac{1}{2n_0}}\text{}\lambda ^2𝒘𝒘\left[\begin{array}{cc}\beta & \mathrm{𝟎}\\ \mathrm{𝟎}& \beta \end{array}\right]\beta _g.`$ (122)
We drop the $`\beta _g`$, as before and then get the following

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\psi \left(𝒓\right)`$ $`=`$ $`\widehat{\mathrm{H}}\psi \left(𝒓\right)`$
$`\widehat{\mathrm{H}}`$ $`=`$ $`n_0\beta +\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`\left(n\left(𝒓\right)n_0\right)\beta \mathrm{i}\text{}\lambda \beta 𝚺𝒖+wide{\displaystyle \frac{1}{2n_0}}\text{}\lambda ^2w^2\beta `$
$`\widehat{𝒪}`$ $`=`$ $`\mathrm{i}\left(M_yp_xM_xp_y\right)`$ (123)
$`=`$ $`\beta \left(𝑴_{}\widehat{𝒑}_{}\right),`$
where, $`w^2=𝒘𝒘`$, the square of the logarithmic gradient of the resistance function. This is how the basic beam optical Hamiltonian (103) gets modified. The next degree of accuracy is achieved by going a step further in the Foldy-Wouthuysen iteration and obtaining the $`\widehat{}_g^{(4)}`$. Then, this would be the higher refined starting beam optical Hamiltonian, further modifying the basic beam optical Hamiltonian (103). This way, we can apply the Foldy-Wouthuysen in cascade to obtain the higher order contributions coming from the logarithmic gradient of the resistance function, to any desired degree of accuracy. We are very unlikely to need any of these contributions, but it is very much possible to keep track of them.

## 4 Applications

In the previous sections we presented the exact matrix representation of the Maxwell equations in a medium with varying permittivity and permeability following the recipe in . From this we derived an exact optical Hamiltonian, which was shown to be in close algebraic analogy with the Dirac equation. This enabled us to apply the machinery of the Foldy-Wouthuysen transformation and we obtained an expansion for the beam-optical Hamiltonian which works to all orders. Formal expressions were obtained for the paraxial Hamiltonian and the leading order aberrating Hamiltonian, without assuming any form for the refractive index. Even at the paraxial level the wavelength-dependent effects manifest by the presence of a matrix term coupled to the logarithmic gradient of the refractive index. This matrix term is very similar to the spin term in the Dirac equation and we call it as the polarizing term in our formalism. The aberrating Hamiltonian contains numerous wavelength-dependent terms in two guises: One of these is the explicit wavelength-dependent terms coming from the commutators inbuilt in the formalism with $`\lambda `$ playing the role played by $`\mathrm{}`$ in quantum mechanics. The other set arises from the the polarizing term.

Now, we apply the formalism to specific examples. One is the medium with constant refractive index. This is perhaps the only problem which can be solved exactly in a closed form expression. This is just to illustrate how the aberration expansion in our formalism can be summed to give the familiar exact result.

The next example is that of the axially symmetric graded index medium. This example enables us to demonstrate the power of the formalism, reproducing the familiar results from the traditional approaches and further giving rise to new results, dependent on the wavelength.

### 4.1 Medium with Constant Refractive Index

Constant refractive index is the simplest possible system. In our formalism, this is perhaps the only case where it is possible to do an exact diagonalization. This is very similar to the exact diagonalization of the free Dirac Hamiltonian. From the experience of the Dirac theory we know that there are hardly any situations where one can do the exact diagonalization. One necessarily has to resort to some approximate diagonalization procedure. The Foldy-Wouthuysen transformation scheme provides the most convenient and accurate diagonalization to any desired degree of accuracy. So we have adopted the Foldy-Wouthuysen scheme in our formalism.

For a medium with constant refractive index, $`n\left(𝒓\right)=n_c`$, we have,

$`\widehat{\mathrm{H}}_c`$ $`=`$ $`n_c\beta +\mathrm{i}\left(M_yp_xM_xp_y\right),`$ (124)
which is exactly diagonalized by the following transform,

$`T^\pm `$ $`=`$ $`\mathrm{exp}\left[\mathrm{i}\left(\pm \mathrm{i}\beta \right)\widehat{𝒪}\theta \right]`$ (125)
$`=`$ $`\mathrm{exp}\left[\mathrm{i}\beta \left(M_yp_xM_xp_y\right)\theta \right]`$
$`=`$ $`\mathrm{cosh}\left(\left|\widehat{𝒑}_{}\right|\theta \right)\mathrm{i}{\displaystyle \frac{\beta \left(M_yp_xM_xp_y\right)}{\left|\widehat{𝒑}_{}\right|}}\mathrm{sinh}\left(\left|\widehat{𝒑}_{}\right|\theta \right).`$
We choose,

$`\mathrm{tanh}\left(2\left|\widehat{𝒑}_{}\right|\theta \right)={\displaystyle \frac{\left|\widehat{𝒑}_{}\right|}{n_c}},`$ (126)
then

$`T^\pm ={\displaystyle \frac{\left(n_c+P_z\right)\mathrm{i}\beta \left(M_yp_xM_xp_y\right)}{\sqrt{2P_z\left(n_c+P_z\right)}}},`$ (127)
where $`P_z=+\sqrt{\left(n_c^2\widehat{𝒑}_{}^2\right)}`$. Then we obtain,

$`\widehat{\mathrm{H}}_c^{\mathrm{diagonal}}`$ $`=`$ $`T^+\widehat{\mathrm{H}}_cT^{}`$ (128)
$`=`$ $`T^+\left\{n_c\beta +\mathrm{i}\left(M_yp_xM_xp_y\right)\right\}T^{}`$
$`=`$ $`\left\{n_c^2\widehat{𝒑}_{}^2\right\}^{\frac{1}{2}}\beta .`$
We next, compare the exact result thus obtained with the approximate ones, obtained through the systematic series procedure we have developed.

$`\widehat{}_c^{(4)}`$ $`=`$ $`n_c\left\{1{\displaystyle \frac{1}{2n_c^2}}\widehat{𝒑}_{}^2{\displaystyle \frac{1}{8n_c^4}}\widehat{𝒑}_{}^4\mathrm{}\right\}\beta `$ (129)
$``$ $`n_c\left\{1{\displaystyle \frac{1}{n_c^2}}\widehat{𝒑}_{}^2\right\}^{\frac{1}{2}}\beta `$
$`=`$ $`\left\{n_c^2\widehat{𝒑}_{}^2\right\}^{\frac{1}{2}}\beta `$
$`=`$ $`\widehat{\mathrm{H}}_c^{\mathrm{diagonal}}.`$
Knowing the Hamiltonian, we can compute the transfer maps. The transfer operator between any pair of points $`\left\{(z^{\prime \prime },z^{})\right|z^{\prime \prime }>z^{}\}`$ on the $`z`$-axis, is formally given by

$`\left|\psi (z^{\prime \prime },z^{})\right|=\widehat{𝒯}(z^{\prime \prime },z^{})|\psi (z^{\prime \prime },z^{}),`$ (130)
with

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\widehat{𝒯}(z^{\prime \prime },z^{})=\widehat{}\widehat{𝒯}(z^{\prime \prime },z^{}),\widehat{𝒯}(z^{\prime \prime },z^{})=\widehat{},`$
$`\widehat{𝒯}(z^{\prime \prime },z^{})=\mathrm{}\left\{\mathrm{exp}\left[{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\widehat{}(z)\right]\right\}`$
$`=\widehat{}{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\widehat{}(z)`$
$`+\left({\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\right)^2{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z{\displaystyle _z^{}^z}𝑑z^{}\widehat{}(z)\widehat{}(z^{})`$
$`+\mathrm{},`$ (131)
where $`\widehat{}`$ is the identity operator and $`\mathrm{}`$ denotes the path-ordered exponential. There is no closed form expression for $`\widehat{𝒯}(z^{\prime \prime },z^{})`$ for an arbitrary choice of the refractive index $`n(𝒓)`$. In such a situation the most convenient form of the expression for the $`z`$-evolution operator $`\widehat{𝒯}(z^{\prime \prime },z^{})`$, or the $`z`$-propagator, is

$$\widehat{𝒯}(z^{\prime \prime },z^{})=\mathrm{exp}\left[\frac{\mathrm{i}}{\text{}\lambda }\widehat{T}(z^{\prime \prime },z^{})\right],$$
(132)
with

$`\widehat{T}(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\widehat{}(z)`$ (133)
$`+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\right){\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z{\displaystyle _z^{}^z}𝑑z^{}[\widehat{}(z),\widehat{}(z^{})]`$
$`+\mathrm{},`$
as given by the Magnus formula , which is described in Appendix-D. We shall be needing these expressions in the next example where the refractive index is not a constant.

Using the procedure outlined above we compute the transfer operator,

$`\widehat{U}_c(z_{\mathrm{out}},z_{\mathrm{in}})=\mathrm{exp}\left[{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\mathrm{\Delta }z_c\right]`$
$`=\mathrm{exp}\left[+{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}n_c\mathrm{\Delta }z\left\{1{\displaystyle \frac{1}{2}}{\displaystyle \frac{\widehat{p}_{}^2}{n_c^2}}{\displaystyle \frac{1}{8}}\left({\displaystyle \frac{\widehat{p}_{}^2}{n_c^2}}\right)^2\mathrm{}\right\}\right],`$ (134)
where, $`\mathrm{\Delta }z=(z_{\mathrm{out}},z_{\mathrm{in}})`$. Using (134), we compute the transfer maps

$`\left(\begin{array}{c}𝒓_{}\\ 𝒑_{}\end{array}\right)_{\mathrm{out}}=\left(\begin{array}{ccc}1& & \frac{1}{\sqrt{n_c^2𝒑_{}^2}}\mathrm{\Delta }z\\ 0& & 1\end{array}\right)\left(\begin{array}{c}𝒓_{}\\ 𝒑_{}\end{array}\right)_{\mathrm{in}}.`$ (141)
The beam-optical Hamiltonian is intrinsically aberrating. Even for the simplest situation of a constant refractive index, we have aberrations to all orders.

### 4.2 Axially Symmetric Graded Index Medium

We just saw the treatment of the medium with a constant refractive index. This is perhaps the only problem which can be solved exactly in a closed form expression. This was just to illustrate how the aberration expansion in our formalism can be obtained. We now consider the next example. The refractive index of an axially symmetric graded-index material can be most generally described by the following polynomial (see, pp. 117 in )

$`n\left(𝒓\right)=n_0+\alpha _2(z)𝒓_{}^2+\alpha _4(z)𝒓_{}^4+\mathrm{},`$ (142)
where, we have assumed the axis of symmetry to coincide with the optic-axis, namely the $`z`$-axis without any loss of generality. We note,

$`\widehat{}`$ $`=`$ $`\{\alpha _2(z)𝒓_{}^2+\alpha _4(z)𝒓_{}^4+\mathrm{},\}\beta \mathrm{i}\text{}\lambda \beta 𝚺𝒖`$
$`\widehat{𝒪}`$ $`=`$ $`\mathrm{i}\left(M_yp_xM_xp_y\right)`$ (143)
$`=`$ $`\beta \left(𝑴_{}\widehat{𝒑}_{}\right)`$
where

$`𝚺𝒖`$ $`=`$ $`{\displaystyle \frac{1}{n_0}}\alpha _2(z)𝚺_{}𝒓_{}{\displaystyle \frac{1}{2n_0}}\left({\displaystyle \frac{d}{dz}}\alpha _2(z)\right)\mathrm{\Sigma }_z𝒓_{}^2`$ (144)
To simplify the formal expression for the beam-optical Hamiltonian $`\widehat{}^{(4)}`$ given in (107)-(108), we make use of the following:

$`\left(𝑴_{}\widehat{𝒑}_{}\right)^2`$ $`=`$ $`\widehat{𝒑}_{}^2,\widehat{𝒪}^2=\widehat{𝒑}_{}^2,{\displaystyle \frac{}{z}}\widehat{𝒪}=0,`$
$`\left(𝑴_{}\widehat{𝒑}_{}\right)𝒓_{}^2\left(𝑴_{}\widehat{𝒑}_{}\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(𝒓_{}^2\widehat{𝒑}_{}^2+\widehat{𝒑}_{}^2𝒓_{}^2\right)+2\text{}\lambda \beta \widehat{L}_z+2\text{}\lambda ^2,`$ (145)
where, $`\widehat{L}_z`$ is the angular momentum. Finally, the beam-optical Hamiltonian to order $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)^2`$ is

$`\widehat{}`$ $`=`$ $`\widehat{H}_{0,p}+\widehat{H}_{0,(4)}`$
$`+\widehat{H}_{0,(2)}^{(\text{}\lambda )}+\widehat{H}_{0,(4)}^{(\text{}\lambda )}`$
$`+\widehat{H}^{(\text{}\lambda ,\sigma )}`$
$`\widehat{H}_{0,p}`$ $`=`$ $`n_0+{\displaystyle \frac{1}{2n_0}}\widehat{𝒑}_{}^2\alpha _2(z)𝒓_{}^2`$
$`\widehat{H}_{0,(4)}`$ $`=`$ $`{\displaystyle \frac{1}{8n_0^3}}\widehat{𝒑}_{}^4`$
$`{\displaystyle \frac{\alpha _2(z)}{4n_0^2}}\left(𝒓_{}^2\widehat{𝒑}_{}^2+\widehat{𝒑}_{}^2𝒓_{}^2\right)`$
$`\alpha _4(z)𝒓_{}^4`$
$`\widehat{H}_{0,(2)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda ^2}{2n_0^2}}\alpha _2(z){\displaystyle \frac{\text{}\lambda }{2n_0^2}}\alpha _2(z)\widehat{L}_z+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2^2(z)𝒓_{}^2`$
$`\widehat{H}_{0,(4)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda }{4n_0^3}}\alpha _2^2(z)\left(𝒓_{}^2\widehat{L}_z+\widehat{L}_z𝒓_{}^2\right)+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)𝒓_{}^4`$
$`\widehat{H}^{(\text{}\lambda ,\sigma )}`$ $`=`$ $`{\displaystyle \frac{\mathrm{i}\text{}\lambda ^3}{2n_0^3}}\left\{{\displaystyle \frac{d}{dz}}\alpha _2(z)\right\}\beta \mathrm{\Sigma }_z`$ (146)
$`+{\displaystyle \frac{\mathrm{i}\text{}\lambda ^2}{4n_0^3}}\alpha _2(z)\left(\mathrm{\Sigma }_xp_y\mathrm{\Sigma }_yp_x\right)`$
$`+{\displaystyle \frac{\mathrm{i}\text{}\lambda ^3}{2n_0^3}}\left\{{\displaystyle \frac{d}{dz}}\alpha _2(z)\right\}\mathrm{\Sigma }_z\widehat{L}_z`$
$`wide+{\displaystyle \frac{\mathrm{i}\text{}\lambda }{4n_0^3}}\alpha _2(z)\beta [𝚺_{}𝒓_{},\widehat{𝒑}_{}^2]_+`$
$`+{\displaystyle \frac{\mathrm{i}\text{}\lambda }{8n_0^3}}\left\{{\displaystyle \frac{d}{dz}}\alpha _2(z)\right\}\beta \mathrm{\Sigma }_z[𝒓_{}^2,\widehat{𝒑}_{}^2]_+`$
$`+\mathrm{}`$
where $`[A,B]_+=(AB+BA)`$ and ‘$`\mathrm{}`$’ are the numerous other terms arising from the polarization term. We have retained only the leading order of such terms above for an illustration. All these matrix terms, related to the polarization will be addressed elsewhere.

The reasons for partitioning the beam-optical Hamiltonian $`\widehat{}`$ in the above manner are as follows. The paraxial Hamiltonian, $`\widehat{H}_{0,p}`$, describes the ideal behaviour. $`\widehat{H}_{0,(4)}`$ is responsible for the third-order aberrations. Both of these Hamiltonians are modified by the wavelength-dependent contributions given in $`\widehat{H}_{0,(2)}^{(\text{}\lambda )}`$ and $`\widehat{H}_{0,(4)}^{(\text{}\lambda )}`$ respectively. Lastly, we have $`\widehat{H}^{(\text{}\lambda ,\sigma )}`$, which is associated with the polarization and shall be examined elsewhere.

#### 4.2.1 Image Rotation

From these sub-Hamiltonians we make several observations:

The term $`\frac{\text{}\lambda }{2n_0^2}\alpha _2(z)\widehat{L}_z`$ which contributes to the paraxial Hamiltonian, gives rise to an image rotation by an angle $`\theta (z)`$:

$`\theta (z^{\prime \prime },z^{})={\displaystyle \frac{\text{}\lambda }{2n_0^2}}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\alpha _2(z).`$ (147)
This image rotation (which need not be small) has no analogue in the square-root approach  and the scalar approach .

#### 4.2.2 Aberrations

The Hamiltonian $`\widehat{H}_{0,(4)}`$ is the one we have in the traditional prescriptions and is responsible for the six aberrations. $`\widehat{H}_{0,(4)}^{(\text{}\lambda )}`$ modifies the above six aberrations by wavelength-dependent contributions and further gives rise to the remaining three aberrations permitted by the axial symmetry. Before proceeding further we enumerate all the nine aberrations permitted by the axial symmetry. The axial symmetry permits exactly nine third-order aberrations which are:

| Symbol | Polynomial | Name |
| --- | --- | --- |
| $`C`$ | $`\widehat{𝒑}_{}^4`$ | Spherical Aberration |
| $`K`$ | $`[\widehat{𝒑}_{}^2,\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)]_+`$ | Coma |
| $`k`$ | $`\widehat{𝒑}_{}^2\widehat{L}_z`$ | Anisotropic Coma |
| $`A`$ | $`\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)^2`$ | Astigmatism |
| $`a`$ | $`\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)\widehat{L}_z`$ | Anisotropic Astigmatism |
| $`F`$ | $`\left(\widehat{𝒑}_{}^2𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^2\right)`$ | Curvature of Field |
| $`D`$ | $`[𝒓_{}^2,\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)]_+`$ | Distortion |
| $`d`$ | $`𝒓_{}^2\widehat{L}_z`$ | Anisotropic Distortion |
| $`E`$ | $`𝒓_{}^4`$ | Nameless? or POCUS |

The name POCUS is used in  on page 137.

The axial symmetry allows only the terms (in the Hamiltonian) which are produced out of, $`\widehat{𝒑}_{}^2`$, $`𝒓_{}^2`$, $`\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)`$ and $`\widehat{L}_z`$. Combinatorially, to fourth-order one would get ten terms including $`\widehat{L}_z^2`$. We have listed nine of them in the table above. The tenth one namely,

$`\widehat{L}_z^2={\displaystyle \frac{1}{2}}\left(\widehat{𝒑}_{}^2𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^2\right){\displaystyle \frac{1}{4}}\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)^2+\text{}\lambda ^2`$ (148)
So, $`\widehat{L}_z^2`$ is not listed separately. Hence, we have only nine third-order aberrations permitted by axial symmetry, as stated earlier.

The paraxial transfer maps are given by

$`\left(\begin{array}{c}𝒓_{}\\ 𝒑_{}\end{array}\right)_{\mathrm{out}}=\left(\begin{array}{cc}P& Q\\ R& S\end{array}\right)\left(\begin{array}{c}𝒓_{}\\ 𝒑_{}\end{array}\right)_{\mathrm{in}},`$ (155)
where $`P`$, $`Q`$, $`R`$ and $`S`$ are the solutions of the paraxial Hamiltonian (146). The symplecticity condition tells us that $`PSQR=1`$. In this particular case from the structure of the paraxial equations, we can further conclude that: $`R=P^{}`$ and $`S=Q^{}`$ where $`()^{}`$ denotes the $`z`$-derivative.

The transfer operator is most accurately and neatly expressed in terms of the paraxial solutions, $`P`$, $`Q`$, $`R`$ and $`S`$, via the interaction picture of the Lie algebraic formulation of light beam optics and charged-particle beam optics .

$`\widehat{𝒯}(z,z_0)`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\widehat{T}(z,z_0)\right],`$ (156)
$`=`$ $`\mathrm{exp}[{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\{C(z^{\prime \prime },z^{})\widehat{𝒑}_{}^4`$
$`+K(z^{\prime \prime },z^{})[\widehat{𝒑}_{}^2,\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)]_+`$
$`+k(z^{\prime \prime },z^{})\widehat{𝒑}_{}^2\widehat{L}_z`$
$`+A(z^{\prime \prime },z^{})\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)^2`$
$`+a(z^{\prime \prime },z^{})\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)\widehat{L}_z`$
$`+F(z^{\prime \prime },z^{})\left(\widehat{𝒑}_{}^2𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^2\right)`$
$`+D(z^{\prime \prime },z^{})[𝒓_{}^2,\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)]_+`$
$`+d(z^{\prime \prime },z^{})𝒓_{}^2\widehat{L}_z`$
$`+E(z^{\prime \prime },z^{})𝒓_{}^4\}].`$
The nine aberration coefficients are given by,

$`C(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}S^4{\displaystyle \frac{\alpha _2(z)}{2n_0^2}}Q^2S^2\alpha _4(z)Q^4`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)Q^4\}`$
$`K(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}RS^3{\displaystyle \frac{\alpha _2(z)}{4n_0^2}}QS(PS+QR)\alpha _4(z)PQ^3`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)PQ^3\}`$
$`k(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda }{2n_0^3}}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\alpha _2^2(z)Q^2`$
$`A(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}R^2S^2{\displaystyle \frac{\alpha _2(z)}{2n_0^2}}PQRS\alpha _4(z)P^2Q^2`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)P^2Q^2\}`$
$`a(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda }{2n_0^3}}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\alpha _2^2(z)PQ`$
$`F(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}R^2S^2{\displaystyle \frac{\alpha _2(z)}{4n_0^2}}(P^2S^2+Q^2R^2)\alpha _4(z)P^2Q^2`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)P^2Q^2\}`$
$`D(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}R^3S{\displaystyle \frac{\alpha _2(z)}{4n_0^2}}PR(PS+QR)\alpha _4(z)P^3Q`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)P^3Q\}`$
$`d(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda }{2n_0^3}}{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\alpha _2^2(z)P^2`$
$`E(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}dz\{{\displaystyle \frac{1}{8n_0^3}}R^4{\displaystyle \frac{\alpha _2(z)}{2n_0^2}}P^2R^2\alpha _4(z)P^4`$ (157)
$`+{\displaystyle \frac{\text{}\lambda ^2}{2n_0^3}}\alpha _2(z)\alpha _4(z)P^4\}.`$
Thus we see that the current approach gives rise to all the nine permissible aberrations. The six aberrations, familiar from the traditional prescriptions get modified by the wavelength-dependent contributions. The extra three ($`k`$, $`a`$ and $`d`$, all anisotropic!) are all pure wavelength-dependent aberrations and totally absent in the traditional square-root approach  and the recently developed scalar approach . A detailed account on the classification of aberrations is available in -.

## 5 Polarization

Let there be Light! (with or/and without polarization) .

## 6 Concluding Remarks

We have developed an exact matrix representation of the Maxwell equations taking into account the spatial and temporal variations of the permittivity and permeability. This representation, using $`8\times 8`$ matrices is the basis for an exact formalism of Maxwell optics presented here. The exact beam optical Hamiltonian, derived from this representation has an algebraic structure in direct correspondence with the Dirac equation of the electron. We exploit this correspondence to adopt the standard machinery, namely the Foldy-Wouthuysen transformation technique of the Dirac theory, to the beam optical formalism. This enabled us to obtain a systematic procedure to obtain the aberration expansion from the beam-optical Hamiltonian to any desired degree of accuracy. We further get the wavelength-dependent contributions at each order, starting with the lowest-order paraxial paraxial Hamiltonian. Formal expressions were obtained for the paraxial and leading order aberrating Hamiltonians, without making any assumption on the form of the varying refractive index.

The beam-optical Hamiltonians also have the wavelength-dependent matrix terms which are associated with the polarization. In this approach we have been able to derive a Hamiltonian which contains both the beam-optics and the polarization.

In Section-IV, we applied the formalism to the specific examples and saw how the beam-optics (paraxial behaviour and the aberrations) gets modified by the wavelength-dependent contributions. First of the two examples is the medium with a constant refractive index. This is perhaps the only problem which can be solved exactly, in a closed form expression. This example is primarily for illustrating certain aspects of the machinery we have used.

The second, and the much more interesting example is that of the axially symmetric graded index medium. For this example, in the traditional approaches one gets only six aberrations. In our formalism we get all the nine aberrations permitted by the axial symmetry. The six aberration coefficients of the traditional approaches get modified by the wavelength-dependent contributions. It is very interesting to note that apart from the wavelength-dependent modifications of the aberrations, this approach also gives rise to the image rotation. This image rotation is proportional to the wavelength and we have derived an explicit relationship for the angle in (147). Such, an image rotation has no analogue/counterpart in any of the traditional prescriptions. It would be worthwhile to experimentally look for the predicted image rotation. The existence of the nine aberrations and image rotation are well-known in axially symmetric magnetic electron lenses, even when treated classically. The quantum treatment of the same system leads to the wavelength-dependent modifications .

The optical Hamiltonian has two components: Beam-Optics and Polarization respectively. We have addressed the former in some detail and the later is in progress. The formalism initiated in this article provides a natural framework for the study of light polarization. This would provide a unified treatment for the beam-optics and the polarization. It also promises a possible generalization of the substitution result in . We shall present this approach elsewhere .

The close analogy between geometrical optics and charged-particle beam optics has been known for too long a time. Until recently it was possible to see this analogy only between the geometrical optics and the classical prescriptions of charge-particle optics. A quantum theory of charged-particle optics was presented in recent years -. With the current development of the non-traditional prescriptions of Helmholtz optics  and the matrix formulation of Maxwell optics presented here, using the rich algebraic machinery of quantum mechanics it is now possible to see a parallel of the analogy at each level. The non-traditional prescription of the Helmholtz optics is in close analogy with the quantum theory of charged-particles based on the Klein-Gordon equation. The matrix formulation of Maxwell optics presented here is in close analogy with the quantum theory of charged-particles based on the Dirac equation . The parallel of these analogies is described in Appendix-E.

An important omission in the present study is the study of the evolution of the fields, $`(𝑬,𝑩)`$, which we shall address in detail elsewhere . Even without the discussion of the fields (as is the case in several other prescriptions) the present study is complete at the Hamiltonian level. We have presented an alternate and exact way of deriving the beam optical Hamiltonian, which reproduces the established results. Furthermore we have derived the extra wavelength-dependent contributions. In the low wavelength limit our formalism reproduces the Lie algebraic formalism of optics. The Foldy-Wouthuysen technique employed by us is ideally suited for the Lie algebraic approach to optics. The present study further strengthens the close analogy between the various prescription of light and charged-particle beam optics .

##

Appendix-A
Riemann-Silberstein Vector

The Riemann-Silberstein complex vector , $`𝑭(𝒓,t)`$ built from the electric field $`𝑫(𝒓,t)`$ and the magnetic filed $`𝑩(𝒓,t)`$ is given by

$`𝑭(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left({\displaystyle \frac{1}{\sqrt{ϵ(𝒓)}}}𝑫(𝒓,t)+\mathrm{i}{\displaystyle \frac{1}{\sqrt{\mu (𝒓)}}}𝑩(𝒓,t)\right)`$ (A.1)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\sqrt{ϵ(𝒓)}𝑬(𝒓,t)+\mathrm{i}{\displaystyle \frac{1}{\sqrt{\mu (𝒓)}}}𝑩(𝒓,t)\right),`$
where $`ϵ(𝒓)`$ is the permittivity of the medium and $`\mu (𝒓)`$ is the permeability of the medium. In vacuum we have $`ϵ_0=8.85\times 10^{12}C^2/N.m^2`$ and $`\mu _0=4\pi \times 10^7N/A^2`$. The Riemann-Silberstein complex vector, $`𝑭(𝒓,t)`$ can also be derived from the potential $`𝒁(𝒓,t)`$, (for example, see ),

$`𝑭(𝒓,t)`$ $`=`$ $`\mathbf{}\times \left\{{\displaystyle \frac{\mathrm{i}}{v}}{\displaystyle \frac{}{t}}𝒁(𝒓,t)+\mathbf{}\times 𝒁(𝒓,t)\right\}.`$ (A.2)
$`𝒁(𝒓,t)`$ is the superpotential and is commonly known as the polarization potential or the Hertz Vector (see ). This further leads to the wave-equation

$`\left\{\mathbf{}^2{\displaystyle \frac{1}{v^2}}{\displaystyle \frac{^2}{t^2}}\right\}𝒁(𝒓,t)`$ $`=`$ $`0.`$ (A.3)
Riemann-Silberstein vector can be used to express many of the quantities associated with the electromagnetic field:

$`\mathrm{Poynting}\mathrm{Vector}:𝑺`$ $`=`$ $`{\displaystyle \frac{1}{\mu }}𝑬\times 𝑩`$
$`=`$ $`\mathrm{i}v\left(𝑭^{}\times 𝑭\right)`$
$`\mathrm{Energy}\mathrm{Density}:u`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(ϵ𝑬𝑬+{\displaystyle \frac{1}{\mu }}𝑩𝑩\right)`$
$`=`$ $`𝑭^{}𝑭`$
$`\mathrm{Momentum}\mathrm{Density}:𝒑_{EB}`$ $`=`$ $`ϵ\left(𝑬\times 𝑩\right)`$
$`=`$ $`{\displaystyle \frac{\mathrm{i}}{v}}\left(𝑭^{}\times 𝑭\right)`$
$`\mathrm{Angular}\mathrm{Momentum}\mathrm{Density}:𝑳_{EB}`$ $`=`$ $`ϵ\left\{𝒓\times \left(𝑬\times 𝑩\right)\right\}`$ (A.4)
$`=`$ $`{\displaystyle \frac{\mathrm{i}}{v}}\left\{𝒓\times \left(𝑭^{}\times 𝑭\right)\right\}.`$
And

$`\mathrm{Total}\mathrm{Energy}:E`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3r\left\{ϵ𝑬𝑬+\frac{1}{\mu }𝑩𝑩\right\}}`$
$`=`$ $`{\displaystyle d^3r\left\{𝑭^{}𝑭\right\}}`$
$`\mathrm{Total}\mathrm{Momentum}:𝑷`$ $`=`$ $`ϵ{\displaystyle d^3r\left\{𝑬\times 𝑩\right\}}`$
$`=`$ $`{\displaystyle \frac{\mathrm{i}}{v}}{\displaystyle d^3r\left\{𝑭^{}\times 𝑭\right\}}`$
$`\mathrm{Total}\mathrm{Angular}\mathrm{Momentum}:𝑴`$ $`=`$ $`ϵ{\displaystyle d^3r\left\{𝒓\times \left(𝑬\times 𝑩\right)\right\}}`$
$`=`$ $`{\displaystyle \frac{\mathrm{i}}{v}}{\displaystyle d^3r\left\{𝒓\times \left(𝑭^{}\times 𝑭\right)\right\}}`$
$`\mathrm{Moment}\mathrm{of}\mathrm{Energy}:𝑵`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3r\left\{𝒓\left(ϵ𝑬𝑬+\frac{1}{\mu }𝑩𝑩\right)\right\}}`$ (A.5)
$`=`$ $`{\displaystyle d^3r\left\{𝒓\left(𝑭^{}𝑭\right)\right\}}`$
In this form these quantities look like the quantum-mechanical expectation values! The use of the Riemann-Silberstein vector as a possible candidate for the photon wavefunction has been advocated for a long time .

##

Appendix-B
An Exact Matrix Representation of the Maxwell Equations in a Medium

Matrix representations of the Maxwell equations are very well-known -. However, all these representations lack an exactness or/and are given in terms of a pair of matrix equations. Some of these representations are in free space. Such a representation is an approximation in a medium with space- and time-dependent permittivity $`ϵ(𝒓,t)`$ and permeability $`\mu (𝒓,t)`$ respectively. Even this approximation is often expressed through a pair of equations using $`3\times 3`$ matrices: one for the curl and one for the divergence which occur in the Maxwell equations. This practice of writing the divergence condition separately is completely avoidable by using $`4\times 4`$ matrices for Maxwell equations in free-space . A single equation using $`4\times 4`$ matrices is necessary and sufficient when $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ are treated as ‘local’ constants .

A treatment taking into account the variations of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ has been presented in . This treatment uses the Riemann-Silberstein vectors, $`𝑭^\pm (𝒓,t)`$ to reexpress the Maxwell equations as four equations: two equations are for the curl and two are for the divergences and there is mixing in $`𝑭^+(𝒓,t)`$ and $`𝑭^{}(𝒓,t)`$. This mixing is very neatly expressed through the two derived functions of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$. These four equations are then expressed as a pair of matrix equations using $`6\times 6`$ matrices: again one for the curl and one for the divergence. Even though this treatment is exact it involves a pair of matrix equations.

Here, we present a treatment which enables us to express the Maxwell equations in a single matrix equation instead of a pair of matrix equations. Our approach is a logical continuation of the treatment in . We use the linear combination of the components of the Riemann-Silberstein vectors, $`𝑭^\pm (𝒓,t)`$ and the final matrix representation is a single equation using $`8\times 8`$ matrices. This representation contains all the four Maxwell equations in presence of sources taking into account the spatial and temporal variations of the permittivity $`ϵ(𝒓,t)`$ and the permeability $`\mu (𝒓,t)`$.

In Section-I we shall summarize the treatment for a homogeneous medium and introduce the required functions and notation. In Section-II we shall present the matrix representation in an inhomogeneous medium, in presence of sources.

### B.1 Homogeneous Medium

We shall start with the Maxwell equations  in an inhomogeneous medium with sources,

$`\mathbf{}𝑫(𝒓,t)=\rho ,`$
$`\mathbf{}\times 𝑯(𝒓,t){\displaystyle \frac{}{t}}𝑫(𝒓,t)=𝑱,`$
$`\mathbf{}\times 𝑬(𝒓,t)+{\displaystyle \frac{}{t}}𝑩(𝒓,t)=0,`$
$`\mathbf{}𝑩(𝒓,t)=0.`$ (B.1)
We assume the media to be linear, that is $`𝑫=ϵ𝑬`$, and $`𝑩=\mu 𝑯`$, where $`ϵ`$ is the permittivity of the medium and $`\mu `$ is the permeability of the medium. In general $`ϵ=ϵ(𝒓,t)`$ and $`\mu =\mu (𝒓,t)`$. In this section we treat them as ‘local’ constants in the various derivations. The magnitude of the velocity of light in the medium is given by $`v(𝒓,t)=\left|𝒗(𝒓,t)\right|=1/\sqrt{ϵ(𝒓,t)\mu (𝒓,t)}`$. In vacuum we have, $`ϵ_0=8.85\times 10^{12}C^2/N.m^2`$ and $`\mu _0=4\pi \times 10^7N/A^2`$.

One possible way to obtain the required matrix representation is to use the Riemann-Silberstein vector  given by

$`𝑭^+(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\sqrt{ϵ(𝒓,t)}𝑬(𝒓,t)+\mathrm{i}{\displaystyle \frac{1}{\sqrt{\mu (𝒓,t)}}}𝑩(𝒓,t)\right)`$
$`𝑭^{}(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\sqrt{ϵ(𝒓,t)}𝑬(𝒓,t)\mathrm{i}{\displaystyle \frac{1}{\sqrt{\mu (𝒓,t)}}}𝑩(𝒓,t)\right).`$ (B.2)
For any homogeneous medium it is equivalent to use either $`𝑭^+(𝒓,t)`$ or $`𝑭^{}(𝒓,t)`$. The two differ by the sign before ‘$`\mathrm{i}`$’ and are not the complex conjugate of one another. We have not assumed any form for $`𝑬(𝒓,t)`$ and $`𝑩(𝒓,t)`$. We will be needing both of them in an inhomogeneous medium, to be considered in detail in Section-II.

If for a certain medium $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ are constants (or can be treated as ‘local’ constants under certain approximations), then the vectors $`𝑭^\pm (𝒓,t)`$ satisfy

$`\mathrm{i}{\displaystyle \frac{}{t}}𝑭^\pm (𝒓,t)`$ $`=`$ $`\pm v\mathbf{}\times 𝑭^\pm (𝒓,t){\displaystyle \frac{1}{\sqrt{2ϵ}}}(\mathrm{i}𝑱)`$
$`\mathbf{}𝑭^\pm (𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2ϵ}}}(\rho ).`$ (B.3)
Thus, by using the Riemann-Silberstein vector it has been possible to reexpress the four Maxwell equations (for a medium with constant $`ϵ`$ and $`\mu `$) as two equations. The first one contains the the two Maxwell equations with curl and the second one contains the two Maxwell with divergences. The first of the two equations in (B.3) can be immediately converted into a $`3\times 3`$ matrix representation. However, this representation does not contain the divergence conditions (the first and the fourth Maxwell equations) contained in the second equation in (B.3). A further compactification is possible only by expressing the Maxwell equations in a $`4\times 4`$ matrix representation. To this end, using the components of the Riemann-Silberstein vector, we define,

$`\mathrm{\Psi }^+(𝒓,t)`$ $`=`$ $`\left[\begin{array}{c}F_x^++\mathrm{i}F_y^+\\ F_z^+\\ F_z^+\\ F_x^++\mathrm{i}F_y^+\end{array}\right],\mathrm{\Psi }^{}(𝒓,t)=\left[\begin{array}{c}F_x^{}\mathrm{i}F_y^{}\\ F_z^{}\\ F_z^{}\\ F_x^{}\mathrm{i}F_y^{}\end{array}\right].`$ (B.12)
The vectors for the sources are

$`W^+`$ $`=`$ $`\left({\displaystyle \frac{1}{\sqrt{2ϵ}}}\right)\left[\begin{array}{c}J_x+\mathrm{i}J_y\\ J_zv\rho \\ J_z+v\rho \\ J_x+\mathrm{i}J_y\end{array}\right],W^{}=\left({\displaystyle \frac{1}{\sqrt{2ϵ}}}\right)\left[\begin{array}{c}J_x\mathrm{i}J_y\\ J_zv\rho \\ J_z+v\rho \\ J_x\mathrm{i}J_y\end{array}\right].`$ (B.21)
Then we obtain

$`{\displaystyle \frac{}{t}}\mathrm{\Psi }^+`$ $`=`$ $`v\left\{𝑴\mathbf{}\right\}\mathrm{\Psi }^+W^+`$
$`{\displaystyle \frac{}{t}}\mathrm{\Psi }^{}`$ $`=`$ $`v\left\{𝑴^{}\mathbf{}\right\}\mathrm{\Psi }^{}W^{},`$ (B.22)
where $`()^{}`$ denotes complex-conjugation and the triplet, $`𝑴=(M_x,M_y,M_z)`$ is expressed in terms of

$`\mathrm{\Omega }`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{𝟎}& 1\mathrm{l}\\ 1\mathrm{l}& \mathrm{𝟎}\end{array}\right],\beta =\left[\begin{array}{cccc}1\mathrm{l}& \mathrm{𝟎}& & \\ \mathrm{𝟎}& 1\mathrm{l}& & \end{array}\right],1\mathrm{l}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].`$ (B.29)
Alternately, we may use the matrix $`J=\mathrm{\Omega }`$. Both differ by a sign. For our purpose it is fine to use either $`\mathrm{\Omega }`$ or $`J`$. However, they have a different meaning: $`J`$ is contravariant and $`\mathrm{\Omega }`$ is covariant; The matrix $`\mathrm{\Omega }`$ corresponds to the Lagrange brackets of classical mechanics and $`J`$ corresponds to the Poisson brackets. An important relation is $`\mathrm{\Omega }=J^1`$. The $`M`$-matrices are:

$`M_x`$ $`=`$ $`\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]=\beta \mathrm{\Omega },`$ (B.34)
$`M_y`$ $`=`$ $`\left[\begin{array}{cccc}0& 0& \mathrm{i}& 0\\ 0& 0& 0& \mathrm{i}\\ \mathrm{i}& 0& 0& 0\\ 0& \mathrm{i}& 0& 0\end{array}\right]=\mathrm{i}\mathrm{\Omega },`$ (B.39)
$`M_z`$ $`=`$ $`\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]=\beta .`$ (B.44)
Each of the four Maxwell equations are easily obtained from the matrix representation in (B.22). This is done by taking the sums and differences of row-I with row-IV and row-II with row-III respectively. The first three give the $`y`$, $`x`$ and $`z`$ components of the curl and the last one gives the divergence conditions present in the evolution equation (B.3).

It is to be noted that the matrices $`𝑴`$ are all non-singular and all are hermitian. Moreover, they satisfy the usual algebra of the Dirac matrices, including,

$`M_x\beta =\beta M_x,`$
$`M_y\beta =\beta M_y,`$
$`M_x^2=M_y^2=M_z^2=I,`$
$`M_xM_y=M_yM_x=\mathrm{i}M_z,`$
$`M_yM_z=M_zM_y=\mathrm{i}M_x,`$
$`M_zM_x=M_xM_z=\mathrm{i}M_y.`$ (B.45)
Before proceeding further we note the following: The pair $`(\mathrm{\Psi }^\pm ,𝑴)`$ are not unique. Different choices of $`\mathrm{\Psi }^\pm `$ would give rise to different $`𝑴`$, such that the triplet $`𝑴`$ continues to to satisfy the algebra of the Dirac matrices in (B.45). We have preferred $`\mathrm{\Psi }^\pm `$ via the the Riemann-Silberstein vector (B.2) in . This vector has certain advantages over the other possible choices. The Riemann-Silberstein vector is well-known in classical electrodynamics and has certain interesting properties and uses .

In deriving the above $`4\times 4`$ matrix representation of the Maxwell equations we have ignored the spatial and temporal derivatives of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ in the first two of the Maxwell equations. We have treated $`ϵ`$ and $`\mu `$ as ‘local’ constants.

### B.2 Inhomogeneous Medium

In the previous section we wrote the evolution equations for the Riemann-Silberstein vector in (B.3), for a medium, treating $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ as ‘local’ constants. From these pairs of equations we wrote the matrix form of the Maxwell equations. In this section we shall write the exact equations taking into account the spatial and temporal variations of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$. It is very much possible to write the required evolution equations using $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$. But we shall follow the procedure in  of using the two derived laboratory functions

$`\mathrm{Velocity}\mathrm{Function}:v(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{ϵ(𝒓,t)\mu (𝒓,t)}}}`$
$`\mathrm{Resistance}\mathrm{Function}:h(𝒓,t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mu (𝒓,t)}{ϵ(𝒓,t)}}}.`$ (B.46)
The function, $`v(𝒓,t)`$ has the dimensions of velocity and the function, $`h(𝒓,t)`$ has the dimensions of resistance (measured in Ohms). We can equivalently use the Conductance Function, $`\kappa (𝒓,t)=1/h(𝒓,t)=ϵ(𝒓,t)/\mu (𝒓,t)`$ (measured in Ohms<sup>-1</sup> or Mhos!) in place of the resistance function, $`h(𝒓,t)`$. These derived functions enable us to understand the dependence of the variations more transparently . Moreover the derived functions are the ones which are measured experimentally. In terms of these functions, $`ϵ=1/\sqrt{vh}`$ and $`\mu =\sqrt{h/v}`$. Using these functions the exact equations satisfied by $`𝑭^\pm (𝒓,t)`$ are

$`\mathrm{i}{\displaystyle \frac{}{t}}𝑭^+(𝒓,t)`$ $`=`$ $`v(𝒓,t)\left(\mathbf{}\times 𝑭^+(𝒓,t)\right)+{\displaystyle \frac{1}{2}}\left(\mathbf{}v(𝒓,t)\times 𝑭^+(𝒓,t)\right)`$
$`+{\displaystyle \frac{v(𝒓,t)}{2h(𝒓)}}\left(\mathbf{}h(𝒓,t)\times 𝑭^{}(𝒓,t)\right){\displaystyle \frac{\mathrm{i}}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}𝑱`$
$`+{\displaystyle \frac{\mathrm{i}}{2}}{\displaystyle \frac{\dot{v}(𝒓,t)}{v(𝒓,t)}}𝑭^+(𝒓,t)+{\displaystyle \frac{\mathrm{i}}{2}}{\displaystyle \frac{\dot{h}(𝒓,t)}{h(𝒓,t)}}𝑭^{}(𝒓,t)`$
$`\mathrm{i}{\displaystyle \frac{}{t}}𝑭^{}(𝒓,t)`$ $`=`$ $`v(𝒓,t)\left(\mathbf{}\times 𝑭^{}(𝒓,t)\right){\displaystyle \frac{1}{2}}\left(\mathbf{}v(𝒓,t)\times 𝑭^{}(𝒓,t)\right)`$
$`{\displaystyle \frac{v(𝒓,t)}{2h(𝒓,t)}}\left(\mathbf{}h(𝒓,t)\times 𝑭^+(𝒓,t)\right){\displaystyle \frac{\mathrm{i}}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}𝑱`$
$`+{\displaystyle \frac{\mathrm{i}}{2}}{\displaystyle \frac{\dot{v}(𝒓,t)}{v(𝒓,t)}}𝑭^{}(𝒓,t)+{\displaystyle \frac{\mathrm{i}}{2}}{\displaystyle \frac{\dot{h}(𝒓,t)}{h(𝒓,t)}}𝑭^+(𝒓,t)`$
$`\mathbf{}𝑭^+(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2v(𝒓,t)}}\left(\mathbf{}v(𝒓,t)𝑭^+(𝒓,t)\right)`$
$`+{\displaystyle \frac{1}{2h(𝒓,t)}}\left(\mathbf{}h(𝒓,t)𝑭^{}(𝒓,t)\right)`$
$`+{\displaystyle \frac{1}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}\rho ,`$
$`\mathbf{}𝑭^{}(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2v(𝒓,t)}}\left(\mathbf{}v(𝒓,t)𝑭^{}(𝒓,t)\right)`$ (B.47)
$`+{\displaystyle \frac{1}{2h(𝒓,t)}}\left(\mathbf{}h(𝒓,t)𝑭^+(𝒓,t)\right)`$
$`+{\displaystyle \frac{1}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}\rho ,`$
where $`\dot{v}=\frac{v}{t}`$ and $`\dot{h}=\frac{h}{t}`$. The evolution equations in (B.47) are exact (for a linear media) and the dependence on the variations of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ has been neatly expressed through the two derived functions. The coupling between $`𝑭^+(𝒓,t)`$ and $`𝑭^{}(𝒓,t)`$. is via the gradient and time-derivative of only one derived function namely, $`h(𝒓,t)`$ or equivalently $`\kappa (𝒓,t)`$. Either of these can be used and both are the directly measured quantities. We further note that the dependence of the coupling is logarithmic

$`{\displaystyle \frac{1}{h(𝒓,t)}}\mathbf{}h(𝒓,t)=\mathbf{}\left\{\mathrm{ln}\left(h(𝒓,t)\right)\right\},{\displaystyle \frac{1}{h(𝒓,t)}}\dot{h}(𝒓,t)={\displaystyle \frac{}{t}}\left\{\mathrm{ln}\left(h(𝒓,t)\right)\right\},`$ (B.48)
where ‘$`\mathrm{ln}`$’ is the natural logarithm.

The coupling can be best summarized by expressing the equations in (B.47) in a (block) matrix form. For this we introduce the following logarithmic function

$`(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{1\mathrm{l}\mathrm{ln}\left(v(𝒓,t)\right)+\sigma _x\mathrm{ln}\left(h(𝒓,t)\right)\right\},`$ (B.49)
where $`\sigma _x`$ is one the triplet of the Pauli matrices

$`𝝈`$ $`=`$ $`\left[\sigma _x=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\sigma _y=\left[\begin{array}{cc}0\hfill & \hfill \mathrm{i}\\ \mathrm{i}\hfill & \hfill 0\end{array}\right],\sigma _z=\left[\begin{array}{cc}1\hfill & \hfill 0\\ 0\hfill & \hfill 1\end{array}\right]\right].`$ (B.56)
Using the above notation the matrix form of the equations in (B.47) is

$`\mathrm{i}\left\{1\mathrm{l}{\displaystyle \frac{}{t}}{\displaystyle \frac{}{t}}\right\}\left[\begin{array}{cc}𝑭^+(𝒓,t)& \\ 𝑭^{}(𝒓,t)& \end{array}\right]`$ $`=`$ $`v(𝒓)\sigma _z\left\{1\mathrm{l}\mathbf{}+\mathbf{}\right\}\times \left[\begin{array}{cc}𝑭^+(𝒓,t)& \\ 𝑭^{}(𝒓,t)& \end{array}\right]`$
$`{\displaystyle \frac{\mathrm{i}}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}𝑱`$
$`\left\{1\mathrm{l}\mathbf{}\mathbf{}\right\}\left[\begin{array}{cc}𝑭^+(𝒓,t)& \\ 𝑭^{}(𝒓,t)& \end{array}\right]`$ $`=`$ $`+{\displaystyle \frac{1}{\sqrt{2}}}\sqrt{v(𝒓,t)h(𝒓,t)}\rho ,`$ (B.64)
where the dot-product and the cross-product are to be understood as

$`\left[\begin{array}{cc}𝑨& 𝑩\\ 𝑪& 𝑫\end{array}\right]\left[\begin{array}{c}𝒖\\ 𝒗\end{array}\right]`$ $`=`$ $`\left[\begin{array}{c}𝑨𝒖+𝑩𝒗\\ 𝑪𝒖+𝑫𝒗\end{array}\right]`$ (B.71)
$`\left[\begin{array}{cc}𝑨& 𝑩\\ 𝑪& 𝑫\end{array}\right]\times \left[\begin{array}{c}𝒖\\ 𝒗\end{array}\right]`$ $`=`$ $`\left[\begin{array}{c}𝑨\times 𝒖+𝑩\times 𝒗\\ 𝑪\times 𝒖+𝑫\times 𝒗\end{array}\right].`$ (B.78)
It is to be noted that the $`6\times 6`$ matrices in the evolution equations in (B.64) are either hermitian or antihermitian. Any dependence on the variations of $`ϵ(𝒓,t)`$ and $`\mu (𝒓,t)`$ is at best ‘weak’. We further note, $`\mathbf{}\left(\mathrm{ln}\left(v(𝒓,t)\right)\right)=\mathbf{}\left(\mathrm{ln}\left(n(𝒓,t)\right)\right)`$ and $`\frac{}{t}\left(\mathrm{ln}\left(v(𝒓,t)\right)\right)=\frac{}{t}\left(\mathrm{ln}\left(n(𝒓,t)\right)\right)`$. In some media, the coupling may vanish ($`\mathbf{}h(𝒓,t)=0`$ and $`\dot{h}(𝒓,t)=0`$) and in the same medium the refractive index, $`n(𝒓,t)=c/v(𝒓,t)`$ may vary ($`\mathbf{}n(𝒓,t)0`$ or/and $`\dot{n}(𝒓,t)0`$). It may be further possible to use the approximations $`\mathbf{}\left(\mathrm{ln}\left(h(𝒓,t)\right)\right)0`$ and $`\frac{}{t}\left(\mathrm{ln}\left(h(𝒓,t)\right)\right)0`$.

We shall be using the following matrices to express the exact representation

$`𝚺=\left[\begin{array}{cc}𝝈& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝝈\end{array}\right],𝜶=\left[\begin{array}{cc}\mathrm{𝟎}& 𝝈\\ 𝝈& \mathrm{𝟎}\end{array}\right],𝑰=\left[\begin{array}{cc}1\mathrm{l}& \mathrm{𝟎}\\ \mathrm{𝟎}& 1\mathrm{l}\end{array}\right],`$ (B.85)
where $`𝚺`$ are the Dirac spin matrices and $`𝜶`$ are the matrices used in the Dirac equation. Then,

$`{\displaystyle \frac{}{t}}\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]{\displaystyle \frac{\dot{v}(𝒓,t)}{2v(𝒓,t)}}\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (B.94)
$`+{\displaystyle \frac{\dot{h}(𝒓,t)}{2h(𝒓,t)}}\left[\begin{array}{cc}\mathrm{𝟎}& \mathrm{i}\beta \alpha _y\\ \mathrm{i}\beta \alpha _y& \mathrm{𝟎}\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (B.99)
$`=v(𝒓,t)\left[\begin{array}{ccc}\left\{𝑴\mathbf{}+𝚺𝒖\right\}& & \mathrm{i}\beta \left(𝚺𝒘\right)\alpha _y\\ \mathrm{i}\beta \left(𝚺^{}𝒘\right)\alpha _y& & \left\{𝑴^{}\mathbf{}+𝚺^{}𝒖\right\}\end{array}\right]\left[\begin{array}{cc}\mathrm{\Psi }^+& \\ \mathrm{\Psi }^{}& \end{array}\right]`$ (B.104)
$`\left[\begin{array}{cc}𝑰& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝑰\end{array}\right]\left[\begin{array}{c}W^+\\ W^{}\end{array}\right],`$ (B.109)
where

$`𝒖(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2v(𝒓,t)}}\mathbf{}v(𝒓,t)={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}v(𝒓,t)\right\}={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}n(𝒓,t)\right\}`$
$`𝒘(𝒓,t)`$ $`=`$ $`{\displaystyle \frac{1}{2h(𝒓,t)}}\mathbf{}h(𝒓,t)={\displaystyle \frac{1}{2}}\mathbf{}\left\{\mathrm{ln}h(𝒓,t)\right\}`$ (B.110)
The above representation contains thirteen $`8\times 8`$ matrices! Ten of these are hermitian. The exceptional ones are the ones that contain the three components of $`𝒘(𝒓,t)`$, the logarithmic gradient of the resistance function. These three matrices, for the resistance function are antihermitian.

We have been able to express the Maxwell equations in a matrix form in a medium with varying permittivity $`ϵ(𝒓,t)`$ and permeability $`\mu (𝒓,t)`$, in presence of sources. We have been able to do so using a single equation instead of a pair of matrix equations. We have used $`8\times 8`$ matrices and have been able to separate the dependence of the coupling between the upper components ($`\mathrm{\Psi }^+`$) and the lower components ($`\mathrm{\Psi }^{}`$) through the two laboratory functions. Moreover the exact matrix representation has an algebraic structure very similar to the Dirac equation. We feel that this representation would be more suitable for some of the studies related to the photon wave function .

##

Appendix-C.
Foldy-Wouthuysen Transformation

In the traditional scheme the purpose of expanding the light optics Hamiltonian $`\widehat{H}=\left(n^2(𝒓)\widehat{𝒑}_{}^2\right)^{1/2}`$ in a series using $`\left(\frac{1}{n_0^2}\widehat{𝒑}_{}^2\right)`$ as the expansion parameter is to understand the propagation of the quasiparaxial beam in terms of a series of approximations (paraxial + nonparaxial). Similar is the situation in the case of the charged-particle optics. Let us recall that in relativistic quantum mechanics too one has a similar problem of understanding the relativistic wave equations as the nonrelativistic approximation plus the relativistic correction terms in the quasirelativistic regime. For the Dirac equation (which is first order in time) this is done most conveniently using the Foldy-Wouthuysen transformation leading to an iterative diagonalization technique.

The main framework of the formalism of optics, used here (and in the charged-particle optics) is based on the transformation technique of the Foldy-Wouthuysen theory which casts the Dirac equation in a form displaying the different interaction terms between the Dirac particle and and an applied electromagnetic field in a nonrelativistic and easily interpretable form (see, , -, for a general discussion of the role of the Foldy-Wouthuysen-type transformations in particle interpretation of relativistic wave equations). The suggestion to employ the Foldy-Wouthuysen Transformation technique in the case of the Helmholtz equation was mentioned in the literature as a remark . It was only in the recent works, that this idea was exploited to analyze the quasiparaxial approximations for for specific beam optical system . The Foldy-Wouthuysen technique is ideally suited for the Lie algebraic approach to optics. With all these plus points, the powerful and ambiguity-free expansion, the Foldy-Wouthuysen Transformation is still little used in optics . In the Foldy-Wouthuysen theory the Dirac equation is decoupled through a canonical transformation into two two-component equations: one reduces to the Pauli equation in the nonrelativistic limit and the other describes the negative-energy states.

Let us describe here briefly the standard Foldy-Wouthuysen theory so that the way it has been adopted for the purposes of the above studies in optics will be clear. Let us consider a charged-particle of rest-mass $`m_0`$, charge $`q`$ in the presence of an electromagnetic field characterized by $`𝑬=\mathbf{}\varphi \frac{}{t}𝑨`$ and $`𝑩=\mathbf{}\times 𝑨`$. Then the Dirac equation is

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }(𝒓,t)`$ $`=`$ $`\widehat{H}_D\mathrm{\Psi }(𝒓,t)`$ (C.1)
$`\widehat{H}_D`$ $`=`$ $`m_0c^2\beta +q\varphi +c𝜶\widehat{𝝅}`$
$`=`$ $`m_0c^2\beta +\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`q\varphi `$
$`\widehat{𝒪}`$ $`=`$ $`c𝜶\widehat{𝝅},`$ (C.2)
where

$`𝜶`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{𝟎}& 𝝈\\ 𝝈& \mathrm{𝟎}\end{array}\right],\beta =\left[\begin{array}{cc}1\mathrm{l}& \mathrm{𝟎}\\ \mathrm{𝟎}& 1\mathrm{l}\end{array}\right],1\mathrm{l}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],`$ (C.6)
$`𝝈`$ $`=`$ $`\left[\sigma _x=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\sigma _y=\left[\begin{array}{cc}0\hfill & \hfill \mathrm{i}\\ \mathrm{i}\hfill & \hfill 0\end{array}\right],\sigma _z=\left[\begin{array}{cc}1\hfill & \hfill 0\\ 0\hfill & \hfill 1\end{array}\right]\right].`$ (C.13)
with $`\widehat{𝝅}=\widehat{𝒑}q𝑨`$, $`\widehat{𝒑}=\mathrm{i}\mathrm{}\mathbf{}`$, and $`\widehat{\pi }^2=\left(\widehat{\pi }_x^2+\widehat{\pi }_y^2+\widehat{\pi }_z^2\right)`$.

In the nonrelativistic situation the upper pair of components of the Dirac Spinor $`\mathrm{\Psi }`$ are large compared to the lower pair of components. The operator $`\widehat{}`$ which does not couple the large and small components of $`\mathrm{\Psi }`$ is called ‘even’ and $`\widehat{𝒪}`$ is called an ‘odd’ operator which couples the large to the small components. Note that

$$\beta \widehat{𝒪}=\widehat{𝒪}\beta ,\beta \widehat{}=\widehat{}\beta .$$
(C.14)
Now, the search is for a unitary transformation, $`\mathrm{\Psi }^{}`$ $`=`$ $`\mathrm{\Psi }`$ $``$ $`\widehat{U}\mathrm{\Psi }`$, such that the equation for $`\mathrm{\Psi }^{}`$ does not contain any odd operator.

In the free particle case (with $`\varphi =0`$ and $`\widehat{𝝅}=\widehat{𝒑}`$) such a Foldy-Wouthuysen transformation is given by

$`\mathrm{\Psi }\mathrm{\Psi }^{}`$ $`=`$ $`\widehat{U}_F\mathrm{\Psi }`$
$`\widehat{U}_F`$ $`=`$ $`e^{\mathrm{i}\widehat{S}}=e^{\beta 𝜶\widehat{𝒑}\theta },\mathrm{tan}\mathrm{\hspace{0.17em}2}|\widehat{𝒑}|\theta ={\displaystyle \frac{|\widehat{𝒑}|}{m_0c}}.`$ (C.15)
This transformation eliminates the odd part completely from the free particle Dirac Hamiltonian reducing it to the diagonal form:

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }^{}`$ $`=`$ $`e^{\mathrm{i}\widehat{S}}\left(m_0c^2\beta +c𝜶\widehat{𝒑}\right)e^{\mathrm{i}\widehat{S}}\mathrm{\Psi }^{}`$ (C.16)
$`=`$ $`\left(\mathrm{cos}|\widehat{𝒑}|\theta +{\displaystyle \frac{\beta 𝜶\widehat{𝒑}}{|\widehat{𝒑}|}}\mathrm{sin}|\widehat{𝒑}|\theta \right)\left(m_0c^2\beta +c𝜶\widehat{𝒑}\right)`$
$`\times \left(\mathrm{cos}|\widehat{𝒑}|\theta {\displaystyle \frac{\beta 𝜶\widehat{𝒑}}{|\widehat{𝒑}|}}\mathrm{sin}|\widehat{𝒑}|\theta \right)\mathrm{\Psi }^{}`$
$`=`$ $`\left(m_0c^2\mathrm{cos}\mathrm{\hspace{0.17em}2}|\widehat{𝒑}|\theta +c|\widehat{𝒑}|\mathrm{sin}\mathrm{\hspace{0.17em}2}|\widehat{𝒑}|\theta \right)\beta \mathrm{\Psi }^{}`$
$`=`$ $`\left(\sqrt{m_0^2c^4+c^2\widehat{p}^2}\right)\beta \mathrm{\Psi }^{}.`$
In the general case, when the electron is in a time-dependent electromagnetic field it is not possible to construct an $`\mathrm{exp}(\mathrm{i}\widehat{S})`$ which removes the odd operators from the transformed Hamiltonian completely. Therefore, one has to be content with a nonrelativistic expansion of the transformed Hamiltonian in a power series in $`1/m_0c^2`$ keeping through any desired order. Note that in the nonrelativistic case, when $`|𝒑|m_0c`$, the transformation operator $`\widehat{U}_F=\mathrm{exp}(\mathrm{i}\widehat{S})`$ with $`\widehat{S}\mathrm{i}\beta \widehat{𝒪}/2m_0c^2`$, where $`\widehat{𝒪}=c𝜶\widehat{𝒑}`$ is the odd part of the free Hamiltonian. So, in the general case we can start with the transformation

$$\mathrm{\Psi }^{(1)}=e^{\mathrm{i}\widehat{S}_1}\mathrm{\Psi },\widehat{S}_1=\frac{\mathrm{i}\beta \widehat{𝒪}}{2m_0c^2}=\frac{\mathrm{i}\beta 𝜶\widehat{𝝅}}{2m_0c}.$$
(C.17)
Then, the equation for $`\mathrm{\Psi }^{(1)}`$ is

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }^{(1)}`$ $`=`$ $`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(e^{\mathrm{i}\widehat{S}_1}\mathrm{\Psi }\right)=\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(e^{\mathrm{i}\widehat{S}_1}\right)\mathrm{\Psi }+e^{\mathrm{i}\widehat{S}_1}\left(\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }\right)`$ (C.18)
$`=`$ $`\left[\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(e^{\mathrm{i}\widehat{S}_1}\right)+e^{\mathrm{i}\widehat{S}_1}\widehat{H}_D\right]\mathrm{\Psi }`$
$`=`$ $`\left[\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(e^{\mathrm{i}\widehat{S}_1}\right)e^{\mathrm{i}\widehat{S}_1}+e^{\mathrm{i}\widehat{S}_1}\widehat{H}_De^{\mathrm{i}\widehat{S}_1}\right]\mathrm{\Psi }^{(1)}`$
$`=`$ $`\left[e^{\mathrm{i}\widehat{S}_1}\widehat{H}_De^{\mathrm{i}\widehat{S}_1}\mathrm{i}\mathrm{}e^{\mathrm{i}\widehat{S}_1}{\displaystyle \frac{}{t}}\left(e^{\mathrm{i}\widehat{S}_1}\right)\right]\mathrm{\Psi }^{(1)}`$
$`=`$ $`\widehat{H}_D^{(1)}\mathrm{\Psi }^{(1)}`$
where we have used the identity $`\frac{}{t}\left(e^{\widehat{A}}\right)e^{\widehat{A}}`$ $`+`$ $`e^{\widehat{A}}\frac{}{t}\left(e^{\widehat{A}}\right)`$ $`=`$ $`\frac{}{t}\widehat{I}`$ $`=0`$.

Now, using the identities

$`e^{\widehat{A}}\widehat{B}e^{\widehat{A}}`$ $`=`$ $`\widehat{B}+[\widehat{A},\widehat{B}]+{\displaystyle \frac{1}{2!}}[\widehat{A},[\widehat{A},\widehat{B}]]+{\displaystyle \frac{1}{3!}}[\widehat{A},[\widehat{A},[\widehat{A},\widehat{B}]]]+\mathrm{}`$ (C.19)
$`e^{\widehat{A}(t)}{\displaystyle \frac{}{t}}\left(e^{\widehat{A}(t)}\right)`$
$`=\left(1+\widehat{A}(t)+{\displaystyle \frac{1}{2!}}\widehat{A}(t)^2+{\displaystyle \frac{1}{3!}}\widehat{A}(t)^3\mathrm{}\right)`$
$`\times {\displaystyle \frac{}{t}}\left(1\widehat{A}(t)+{\displaystyle \frac{1}{2!}}\widehat{A}(t)^2{\displaystyle \frac{1}{3!}}\widehat{A}(t)^3\mathrm{}\right)`$
$`=\left(1+\widehat{A}(t)+{\displaystyle \frac{1}{2!}}\widehat{A}(t)^2+{\displaystyle \frac{1}{3!}}\widehat{A}(t)^3\mathrm{}\right)`$
$`\times ({\displaystyle \frac{\widehat{A}(t)}{t}}+{\displaystyle \frac{1}{2!}}\{{\displaystyle \frac{\widehat{A}(t)}{t}}\widehat{A}(t)+\widehat{A}(t){\displaystyle \frac{\widehat{A}(t)}{t}}\}`$
$`{\displaystyle \frac{1}{3!}}\{{\displaystyle \frac{\widehat{A}(t)}{t}}\widehat{A}(t)^2+\widehat{A}(t){\displaystyle \frac{\widehat{A}(t)}{t}}\widehat{A}(t)`$
$`+\widehat{A}(t)^2{\displaystyle \frac{\widehat{A}(t)}{t}}\}\mathrm{})`$
$`{\displaystyle \frac{\widehat{A}(t)}{t}}{\displaystyle \frac{1}{2!}}[\widehat{A}(t),{\displaystyle \frac{\widehat{A}(t)}{t}}]`$
$`{\displaystyle \frac{1}{3!}}[\widehat{A}(t),[\widehat{A}(t),{\displaystyle \frac{\widehat{A}(t)}{t}}]]`$
$`{\displaystyle \frac{1}{4!}}[\widehat{A}(t),[\widehat{A}(t),[\widehat{A}(t),{\displaystyle \frac{\widehat{A}(t)}{t}}]]],`$
with $`\widehat{A}=\mathrm{i}\widehat{S}_1`$, we find

$`\widehat{H}_D^{(1)}`$ $``$ $`\widehat{H}_D\mathrm{}{\displaystyle \frac{\widehat{S}_1}{t}}+\mathrm{i}[\widehat{S}_1,\widehat{H}_D{\displaystyle \frac{\mathrm{}}{2}}{\displaystyle \frac{\widehat{S}_1}{t}}]`$ (C.20)
$`{\displaystyle \frac{1}{2!}}[\widehat{S}_1,[\widehat{S}_1,\widehat{H}_D{\displaystyle \frac{\mathrm{}}{3}}{\displaystyle \frac{\widehat{S}_1}{t}}]]`$
$`{\displaystyle \frac{\mathrm{i}}{3!}}[\widehat{S}_1,[\widehat{S}_1,[\widehat{S}_1,\widehat{H}_D{\displaystyle \frac{\mathrm{}}{4}}{\displaystyle \frac{\widehat{S}_1}{t}}]]].`$
Substituting in (C.20), $`\widehat{H}_D=m_0c^2\beta +\widehat{}+\widehat{𝒪}`$, simplifying the right hand side using the relations $`\beta \widehat{𝒪}=\widehat{𝒪}\beta `$ and $`\beta \widehat{}=\widehat{}\beta `$ and collecting everything together, we have

$`\widehat{H}_D^{(1)}`$ $``$ $`m_0c^2\beta +\widehat{}_1+\widehat{𝒪}_1`$
$`\widehat{}_1`$ $``$ $`\widehat{}+{\displaystyle \frac{1}{2m_0c^2}}\beta \widehat{𝒪}^2{\displaystyle \frac{1}{8m_0^2c^4}}[\widehat{𝒪},\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}}{t}}\right)]`$
$`{\displaystyle \frac{1}{8m_0^3c^6}}\beta \widehat{𝒪}^4`$
$`\widehat{𝒪}_1`$ $``$ $`{\displaystyle \frac{\beta }{2m_0c^2}}\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}}{t}}\right){\displaystyle \frac{1}{3m_0^2c^4}}\widehat{𝒪}^3,`$ (C.21)
with $`\widehat{}_1`$ and $`\widehat{𝒪}_1`$ obeying the relations $`\beta \widehat{𝒪}_1=\widehat{𝒪}_1\beta `$ and $`\beta \widehat{}_1=\widehat{}_1\beta `$ exactly like $`\widehat{}`$ and $`\widehat{𝒪}`$. It is seen that while the term $`\widehat{𝒪}`$ in $`\widehat{H}_D`$ is of order zero with respect to the expansion parameter $`1/m_0c^2`$ (i.e., $`\widehat{𝒪}`$ $`=`$ $`O\left(\left(1/m_0c^2\right)^0\right)`$ the odd part of $`\widehat{H}_D^{(1)}`$, namely $`\widehat{𝒪}_1`$, contains only terms of order $`1/m_0c^2`$ and higher powers of $`1/m_0c^2`$ (i.e., $`\widehat{𝒪}_1=O\left(\left(1/m_0c^2\right)\right)`$).

To reduce the strength of the odd terms further in the transformed Hamiltonian a second Foldy-Wouthuysen transformation is applied with the same prescription:

$`\mathrm{\Psi }^{(2)}`$ $`=`$ $`e^{\mathrm{i}\widehat{S}_2}\mathrm{\Psi }^{(1)},`$
$`\widehat{S}_2`$ $`=`$ $`{\displaystyle \frac{\mathrm{i}\beta \widehat{𝒪}_1}{2m_0c^2}}`$ (C.22)
$`=`$ $`{\displaystyle \frac{\mathrm{i}\beta }{2m_0c^2}}\left[{\displaystyle \frac{\beta }{2m_0c^2}}\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}}{t}}\right){\displaystyle \frac{1}{3m_0^2c^4}}\widehat{𝒪}^3\right].`$
After this transformation,

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }^{(2)}`$ $`=`$ $`\widehat{H}_D^{(2)}\mathrm{\Psi }^{(2)},\widehat{H}_D^{(2)}=m_0c^2\beta +\widehat{}_2+\widehat{𝒪}_2`$
$`\widehat{}_2wide`$ $``$ $`\widehat{}_1,\widehat{𝒪}_2{\displaystyle \frac{\beta }{2m_0c^2}}\left([\widehat{𝒪}_1,\widehat{}_1]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}_1}{t}}\right),`$ (C.23)
where, now, $`\widehat{𝒪}_2=O\left(\left(1/m_0c^2\right)^2\right)`$. After the third transformation

$$\mathrm{\Psi }^{(3)}=e^{\mathrm{i}\widehat{S}_3}\mathrm{\Psi }^{(2)},\widehat{S}_3=\frac{\mathrm{i}\beta \widehat{𝒪}_2}{2m_0c^2}$$

we have

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }^{(3)}`$ $`=`$ $`\widehat{H}_D^{(3)}\mathrm{\Psi }^{(3)},\widehat{H}_D^{(3)}=m_0c^2\beta +\widehat{}_3+\widehat{𝒪}_3`$
$`\widehat{}_3`$ $``$ $`\widehat{}_2\widehat{}_1,\widehat{𝒪}_3{\displaystyle \frac{\beta }{2m_0c^2}}\left([\widehat{𝒪}_2,\widehat{}_2]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}_2}{t}}\right),`$ (C.24)
where $`\widehat{𝒪}_3=O\left(\left(1/m_0c^2\right)^3\right)`$. So, neglecting $`\widehat{𝒪}_3`$,

$`\widehat{H}_D^{(3)}`$ $``$ $`m_0c^2\beta +\widehat{}+{\displaystyle \frac{1}{2m_0c^2}}\beta \widehat{𝒪}^2`$ (C.25)
$`{\displaystyle \frac{1}{8m_0^2c^4}}[\widehat{𝒪},\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}}{t}}\right)]`$
$`{\displaystyle \frac{1}{8m_0^3c^6}}\beta \left\{\widehat{𝒪}^4+\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}}{t}}\right)^2\right\}`$
It may be noted that starting with the second transformation successive $`(\widehat{},\widehat{𝒪})`$ pairs can be obtained recursively using the rule

$`\widehat{}_j`$ $`=`$ $`\widehat{}_1\left(\widehat{}\widehat{}_{j1},\widehat{𝒪}\widehat{𝒪}_{j1}\right)`$
$`\widehat{𝒪}_j`$ $`=`$ $`\widehat{𝒪}_1\left(\widehat{}\widehat{}_{j1},\widehat{𝒪}\widehat{𝒪}_{j1}\right),j>1,`$ (C.26)
and retaining only the relevant terms of desired order at each step.

With $`\widehat{}=q\varphi `$ and $`\widehat{𝒪}=c𝜶\widehat{𝝅}`$, the final reduced Hamiltonian (C.25) is, to the order calculated,

$`\widehat{H}_D^{(3)}`$ $`=`$ $`\beta \left(m_0c^2+{\displaystyle \frac{\widehat{\pi }^2}{2m_0}}{\displaystyle \frac{\widehat{p}^4}{8m_0^3c^6}}\right)+q\varphi {\displaystyle \frac{q\mathrm{}}{2m_0c}}\beta 𝚺𝑩`$ (C.27)
$`{\displaystyle \frac{\mathrm{i}q\mathrm{}^2}{8m_0^2c^2}}𝚺\mathrm{curl}𝑬{\displaystyle \frac{q\mathrm{}}{4m_0^2c^2}}𝚺𝑬\times \widehat{𝒑}`$
$`{\displaystyle \frac{q\mathrm{}^2}{8m_0^2c^2}}\mathrm{div}𝑬,`$
with the individual terms having direct physical interpretations. The terms in the first parenthesis result from the expansion of $`\sqrt{m_0^2c^4+c^2\widehat{\pi }^2}`$ showing the effect of the relativistic mass increase. The second and third terms are the electrostatic and magnetic dipole energies. The next two terms, taken together (for hermiticity), contain the spin-orbit interaction. The last term, the so-called Darwin term, is attributed to the zitterbewegung (trembling motion) of the Dirac particle: because of the rapid coordinate fluctuations over distances of the order of the Compton wavelength ($`2\pi \mathrm{}/m_0c`$) the particle sees a somewhat smeared out electric potential.

It is clear that the Foldy-Wouthuysen transformation technique expands the Dirac Hamiltonian as a power series in the parameter $`1/m_0c^2`$ enabling the use of a systematic approximation procedure for studying the deviations from the nonrelativistic situation. We note the analogy between the nonrelativistic particle dynamics and paraxial optics:

The Analogy

| Standard Dirac Equation | Beam Optical Form |
| --- | --- |
| $`m_0c^2\beta +\widehat{}_D+\widehat{𝒪}_D`$ | $`n_0\beta +\widehat{}+\widehat{𝒪}`$ |
| $`m_0c^2`$ | $`n_0`$ |
| Positive Energy | Forward Propagation |
| Nonrelativistic, $`\left|𝝅\right|m_0c`$ | Paraxial Beam, $`\left|\widehat{𝒑}_{}\right|n_0`$ |
| Non relativistic Motion | Paraxial Behavior |
| \+ Relativistic Corrections | \+ Aberration Corrections |

Noting the above analogy, the idea of Foldy-Wouthuysen form of the Dirac theory has been adopted to study the paraxial optics and deviations from it by first casting the Maxwell equations in a spinor form resembling exactly the Dirac equation (C.1, C.2) in all respects: i.e., a multicomponent $`\mathrm{\Psi }`$ having the upper half of its components large compared to the lower components and the Hamiltonian having an even part $`(\widehat{})`$, an odd part $`(\widehat{𝒪})`$, a suitable expansion parameter, ($`|\widehat{𝒑}_{}|/n_01`$) characterizing the dominant forward propagation and a leading term with a $`\beta `$ coefficient commuting with $`\widehat{}`$ and anticommuting with $`\widehat{𝒪}`$. The additional feature of our formalism is to return finally to the original representation after making an extra approximation, dropping $`\beta `$ from the final reduced optical Hamiltonian, taking into account the fact that we are primarily interested only in the forward-propagating beam.

##

Appendix-D
The Magnus Formula

The Magnus formula is the continuous analogue of the famous Baker-Campbell-Hausdorff (BCH) formula

$$\mathrm{e}^{\widehat{A}}\mathrm{e}^{\widehat{B}}=\mathrm{e}^{\widehat{A}+\widehat{B}+\frac{1}{2}[\widehat{A},\widehat{B}]+\frac{1}{12}\left\{[[\widehat{A},\widehat{A}],\widehat{B}]+[[\widehat{A},\widehat{B}],\widehat{B}]\right\}+\mathrm{}}.$$
(D.1)
Let it be required to solve the differential equation

$$\frac{}{t}u(t)=\widehat{A}(t)u(t)$$
(D.2)
to get $`u(T)`$ at $`T>t_0`$, given the value of $`u(t_0)`$; the operator $`\widehat{A}`$ can represent any linear operation. For an infinitesimal $`\mathrm{\Delta }t`$, we can write

$$u(t_0+\mathrm{\Delta }t)=e^{\mathrm{\Delta }t\widehat{A}(t_0)}u(t_0).$$
(D.3)
Iterating this solution we have

$`u(t_0+2\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{e}^{\mathrm{\Delta }t\widehat{A}(t_0+\mathrm{\Delta }t)}\mathrm{e}^{\mathrm{\Delta }t\widehat{A}(t_0)}u(t_0)`$
$`u(t_0+3\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{e}^{\mathrm{\Delta }t\widehat{A}(t_0+2\mathrm{\Delta }t)}\mathrm{e}^{\mathrm{\Delta }t\widehat{A}(t_0+\mathrm{\Delta }t)}\mathrm{e}^{\mathrm{\Delta }t\widehat{A}(t_0)}u(t_0)`$ (D.4)
$`\mathrm{}\mathrm{and}\mathrm{so}\mathrm{on}.`$
If $`T=t_0+N\mathrm{\Delta }t`$ we would have

$$u(T)=\left\{\underset{n=0}{\overset{N1}{}}e^{\mathrm{\Delta }t\widehat{A}(t_0+n\mathrm{\Delta }t)}\right\}u(t_0).$$
(D.5)
Thus, $`u(T)`$ is given by computing the product in (D.5) using successively the BCH-formula (D.1) and considering the limit $`\mathrm{\Delta }t0,N\mathrm{}`$ such that $`N\mathrm{\Delta }t=Tt_0`$. The resulting expression is the Magnus formula (Magnus, ) :

$`u(T)`$ $`=`$ $`\widehat{𝒯}(T,t_0)u(t_0)`$
$`𝒯(T,t_0)`$ $`=`$ $`\mathrm{exp}\{{\displaystyle _{t_0}^T}dt_1\widehat{A}(t_1)`$ (D.6)
$`+{\displaystyle \frac{1}{2}}{\displaystyle _{t_0}^T}𝑑t_2{\displaystyle _{t_0}^{t_2}}𝑑t_1[\widehat{A}(t_2),\widehat{A}(t_1)]`$
$`+{\displaystyle \frac{1}{6}}{\displaystyle _{t_0}^T}dt_3{\displaystyle _{t_0}^{t_3}}dt_2{\displaystyle _{t_0}^{t_2}}dt_1([[\widehat{A}(t_3),\widehat{A}(t_2)],\widehat{A}(t_1)]`$
$`+[[\widehat{A}(t_1),\widehat{A}(t_2)],\widehat{A}(t_3)])+\mathrm{}\}.`$
To see how the equation (D.6) is obtained let us substitute the assumed form of the solution, $`u(t)=\widehat{𝒯}(t,t_0)u\left(t_0\right)`$, in (D.2). Then, it is seen that $`\widehat{𝒯}(t,t_0)`$ obeys the equation

$$\frac{}{t}\widehat{𝒯}(t,t_0)=\widehat{A}(t)𝒯(t,t_0),\widehat{𝒯}(t_0,t_0)=\widehat{I}.$$
(D.7)
Introducing an iteration parameter $`\lambda `$ in (D.7), let

$`{\displaystyle \frac{}{t}}\widehat{𝒯}(t,t_0;\lambda )`$ $`=`$ $`\lambda \widehat{A}(t)\widehat{𝒯}(t,t_0;\lambda ),`$ (D.8)
$`\widehat{𝒯}(t_0,t_0;\lambda )`$ $`=`$ $`\widehat{I},\widehat{𝒯}(t,t_0;1)=\widehat{𝒯}(t,t_0).`$ (D.9)
Assume a solution of (A8) to be of the form

$$\widehat{𝒯}(t,t_0;\lambda )=\mathrm{e}^{\mathrm{\Omega }(t,t_0;\lambda )}$$
(D.10)
with

$$\mathrm{\Omega }(t,t_0;\lambda )=\underset{n=1}{\overset{\mathrm{}}{}}\lambda ^n\mathrm{\Delta }_n(t,t_0),\mathrm{\Delta }_n(t_0,t_0)=0\mathrm{for}\mathrm{all}n.$$
(D.11)
Now, using the identity (se Wilcox, )

$$\frac{}{t}\mathrm{e}^{\mathrm{\Omega }(t,t_0;\lambda )}=\left\{_0^1𝑑se^{s\mathrm{\Omega }(t,t_0;\lambda )}\frac{}{t}\mathrm{\Omega }(t,t_0;\lambda )\mathrm{e}^{s\mathrm{\Omega }(t,t_0;\lambda )}\right\}\mathrm{e}^{\mathrm{\Omega }(t,\lambda )},$$
(D.12)
one has

$$_0^1𝑑se^{s\mathrm{\Omega }(t,t_0;\lambda )}\frac{}{t}\mathrm{\Omega }(t,t_0;\lambda )e^{s\mathrm{\Omega }(t,t_0;\lambda )}=\lambda \widehat{A}(t).$$
(D.13)
Substituting in (D.13) the series expression for $`\mathrm{\Omega }(t,t_0;\lambda )`$ (D.11), expanding the left hand side using the first identity in (D.8), integrating and equating the coefficients of $`\lambda ^j`$ on both sides, we get, recursively, the equations for $`\mathrm{\Delta }_1(t,t_0)`$, $`\mathrm{\Delta }_2(t,t_0),\mathrm{},`$ etc. For $`j=1`$

$$\frac{}{t}\mathrm{\Delta }_1(t,t_0)=\widehat{A}(t),\mathrm{\Delta }_1(t_0,t_0)=0$$
(D.14)
and hence

$$\mathrm{\Delta }_1(t,t_0)=_{t_0}^t𝑑t_1\widehat{A}(t_1).$$
(D.15)
For $`j=2`$

$$\frac{}{t}\mathrm{\Delta }_2(t,t_0)+\frac{1}{2}[\mathrm{\Delta }_1(t,t_0),\frac{}{t}\mathrm{\Delta }_1(t,t_0)]=0,\mathrm{\Delta }_2(t_0,t_0)=0$$
(D.16)
and hence

$$\mathrm{\Delta }_2(t,t_0)=\frac{1}{2}_{t_0}^t𝑑t_2_{t_0}^{t_2}𝑑t_1[\widehat{A}(t_2),\widehat{A}(t_1)].$$
(D.17)
Similarly,

$`\mathrm{\Delta }_3(t,t_0)`$ $`=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle _{t_0}^t}dt_1{\displaystyle _{t_0}^{t_1}}dt_2{\displaystyle _{t_0}^{t_2}}dt_3\{[[\widehat{A}(t_1),\widehat{A}(t_2)],\widehat{A}(t_3)]`$ (D.18)
$`+[[\widehat{A}(t_3),\widehat{A}(t_2)],\widehat{A}(t_1)]\}.`$
Then, the Magnus formula in (D.6) follows from (D.9)-(D.11). Equation (133) we have, in the context of $`z`$-evolution follows from the above discussion with the identification $`tz`$, $`t_0z^{(1)}`$, $`Tz^{(2)}`$ and $`\widehat{A}(t)\frac{\mathrm{i}}{\mathrm{}}\widehat{}_o(z)`$.

For more details on the exponential solutions of linear differential equations, related operator techniques and applications to physical problems the reader is referred to Wilcox , Bellman and Vasudevan , Dattoli et al. , and references therein.

##

Appendix-E
Analogies between light optics and charged-particle optics: Recent Developments

Historically, variational principles have played a fundamental role in the evolution of mathematical models in classical physics, and many equations can be derived by using them. Here the relevant examples are Fermat’s principle in optics and Maupertuis’ principle in mechanics. The beginning of the analogy between geometrical optics and mechanics is usually attributed to Descartes (1637), but actually it can traced back to Ibn Al-Haitham Alhazen (0965-1037) . The analogy between the trajectory of material particles in potential fields and the path of light rays in media with continuously variable refractive index was formalized by Hamilton in 1833. This Hamiltonian analogy lead to the development of electron optics in 1920s, when Busch derived the focusing action and a lens-like action of the axially symmetric magnetic field using the methodology of geometrical optics. Around the same time Louis de Broglie associated his now famous wavelength to moving particles. Schrödinger extended the analogy by passing from geometrical optics to wave optics through his wave equation incorporating the de Broglie wavelength. This analogy played a fundamental role in the early development of quantum mechanics. The analogy, on the other hand, lead to the development of practical electron optics and one of the early inventions was the electron microscope by Ernst Ruska. A detailed account of Hamilton’s analogy is available in -.

Until very recently, it was possible to see this analogy only between the geometrical-optic and classical prescriptions of electron optics. The reasons being that, the quantum theories of charged-particle beam optics have been under development only for about a decade - with the very expected feature of wavelength-dependent effects, which have no analogue in the traditional descriptions of light beam optics. With the current development of the non-traditional prescriptions of Helmholtz optics  and the matrix formulation of Maxwell optics, accompanied with wavelength-dependent effects, it is seen that the analogy between the two systems persists. The non-traditional prescription of Helmholtz optics is in close analogy with the quantum theory of charged-particle beam optics based on the Klein-Gordon equation. The matrix formulation of Maxwell optics is in close analogy with the quantum theory of charged-particle beam optics based on the Dirac equation. This analogy is summarized in the table of Hamiltonians. In this short note it is difficult to present the derivation of the various Hamiltonians which are available in the references. We shall briefly consider an outline of the quantum prescriptions and the non-traditional prescriptions respectively. A complete coverage to the new field of Quantum Aspects of Beam Physics (QABP), can be found in the proceedings of the series of meetings under the same name .

### E.1 Quantum Formalism of Charged-Particle Beam Optics

The classical treatment of charged-particle beam optics has been extremely successful in the designing and working of numerous optical devices, from electron microscopes to very large particle accelerators. It is natural, however to look for a prescription based on the quantum theory, since any physical system is quantum mechanical at the fundamental level! Such a prescription is sure to explain the grand success of the classical theories. It is sure to help get a deeper understanding and lead to better designing of charged-particle beam devices.

The starting point of the quantum prescription of charged particle beam optics is to build a theory based on the basic equations of quantum mechanics (Schrödinger, Klein-Gordon, Dirac) appropriate to the situation under study. In order to analyze the evolution of the beam parameters of the various individual beam optical elements (quadrupoles, bending magnets, $`\mathrm{}`$) along the optic axis of the system, the first step is to start with the basic time-dependent equations of quantum mechanics and then obtain an equation of the form

$$\mathrm{i}\mathrm{}\frac{}{s}\psi (x,y;s)=\widehat{}(x,y;s)\psi (x,y;s),$$
(E.1)
where $`(x,y;s)`$ constitute a curvilinear coordinate system, adapted to the geometry of the system. Eq. (E.1) is the basic equation in the quantum formalism, called as the beam-optical equation; $``$ and $`\psi `$ as the beam-optical Hamiltonian and the beam wavefunction respectively. The second step requires obtaining a relationship between any relevant observable $`\{O(s)\}`$ at the transverse-plane at $`s`$ and the observable $`\{O(s_{\mathrm{in}})\}`$ at the transverse plane at $`s_{\mathrm{in}}`$, where $`s_{\mathrm{in}}`$ is some input reference point. This is achieved by the integration of the beam-optical equation in (E.1)

$`\psi (x,y;s)`$ $`=`$ $`\widehat{U}(s,s_{\mathrm{in}})\psi (x,y;s_{\mathrm{in}}),`$ (E.2)
which gives the required transfer maps

$`O\left(s_{\mathrm{in}}\right)O\left(s\right)`$ $`=`$ $`\psi (x,y;s)\left|O\right|\psi (x,y;s),`$ (E.3)
$`=`$ $`\psi (x,y;s_{\mathrm{in}})\left|\widehat{U}^{}O\widehat{U}\right|\psi (x,y;s_{\mathrm{in}}).`$
The two-step algorithm stated above gives an over-simplified picture of the quantum formalism. There are several crucial points to be noted. The first step in the algorithm of obtaining the beam-optical equation is not to be treated as a mere transformation which eliminates $`t`$ in preference to a variable $`s`$ along the optic axis. A clever set of transforms are required which not only eliminate the variable $`t`$ in preference to $`s`$ but also give us the $`s`$-dependent equation which has a close physical and mathematical correspondence with the original $`t`$-dependent equation of standard time-dependent quantum mechanics. The imposition of this stringent requirement on the construction of the beam-optical equation ensures the execution of the second-step of the algorithm. The beam-optical equation is such that all the required rich machinery of quantum mechanics becomes applicable to the computation of the transfer maps that characterize the optical system. This describes the essential scheme of obtaining the quantum formalism. The rest is mostly mathematical detail which is inbuilt in the powerful algebraic machinery of the algorithm, accompanied with some reasonable assumptions and approximations dictated by the physical considerations. The nature of these approximations can be best summarized in the optical terminology as a systematic procedure of expanding the beam optical Hamiltonian in a power series of $`|\widehat{𝝅}_{}/p_0|`$, where $`p_0`$ is the design (or average) momentum of beam particles moving predominantly along the direction of the optic axis and $`\widehat{𝝅}_{}`$ is the small transverse kinetic momentum. The leading order approximation along with $`|\widehat{𝝅}_{}/p_0|1`$, constitutes the paraxial or ideal behaviour and higher order terms in the expansion give rise to the nonlinear or aberrating behaviour. It is seen that the paraxial and aberrating behaviour get modified by the quantum contributions which are in powers of the de Broglie wavelength ($`\text{}\lambda _0=\mathrm{}/p_0`$). The classical limit of the quantum formalism reproduces the well known Lie algebraic formalism of charged-particle beam optics .

### E.2 Light Optics: Various Prescriptions

The traditional scalar wave theory of optics (including aberrations to all orders) is based on the beam-optical Hamiltonian derived by using Fermat’s principle. This approach is purely geometrical and works adequately in the scalar regime. The other approach is based on the square-root of the Helmholtz operator, which is derived from the Maxwell equations . This approach works to all orders and the resulting expansion is no different from the one obtained using the geometrical approach of Fermat’s principle. As for the polarization: a systematic procedure for the passage from scalar to vector wave optics to handle paraxial beam propagation problems, completely taking into account the way in which the Maxwell equations couple the spatial variation and polarization of light waves, has been formulated by analyzing the basic Poincaré invariance of the system, and this procedure has been successfully used to clarify several issues in Maxwell optics -.

In the above approaches, the beam-optics and the polarization are studied separately, using very different machineries. The derivation of the Helmholtz equation from the Maxwell equations is an approximation as one neglects the spatial and temporal derivatives of the permittivity and permeability of the medium. Any prescription based on the Helmholtz equation is bound to be an approximation, irrespective of how good it may be in certain situations. It is very natural to look for a prescription based fully on the Maxwell equations, which is sure to provide a deeper understanding of beam-optics and light polarization in a unified manner.

The two-step algorithm used in the construction of the quantum theories of charged-particle beam optics is very much applicable in light optics! But there are some very significant conceptual differences to be borne in mind. When going beyond Fermat’s principle the whole of optics is completely governed by the Maxwell equations, and there are no other equations, unlike in quantum mechanics, where there are separate equations for, spin-$`1/2`$, spin-$`1`$, $`\mathrm{}`$.

Maxwell’s equations are linear (in time and space derivatives) but coupled in the fields. The decoupling leads to the Helmholtz equation which is quadratic in derivatives. In the specific context of beam optics, purely from a calculational point of view, the starting equations are the Helmholtz equation governing scalar optics and for a more accurate prescription one uses the full set of Maxwell equations, leading to vector optics. In the context of the two-step algorithm, the Helmholtz equation and the Maxwell equations in a matrix representation can be treated as the ‘basic’ equations, analogue of the basic equations of quantum mechanics. This works perfectly fine from a calculational point of view in the scheme of the algorithm we have.

Exploiting the similarity between the Helmholtz wave equation and the Klein-Gordon equation, the former is linearized using the Feshbach-Villars procedure used for the linearization of the Klein-Gordon equation. Then the Foldy-Wouthuysen iterative diagonalization technique is applied to obtain a Hamiltonian description for a system with varying refractive index. This technique is an alternative to the conventional method of series expansion of the radical. Besides reproducing all the traditional quasiparaxial terms, this method leads to additional terms, which are dependent on the wavelength, in the optical Hamiltonian. This is the non-traditional prescription of scalar optics.

The Maxwell equations can be cast into an exact matrix form taking into account the spatial and temporal variations of the permittivity and permeability. The derived representation using $`8\times 8`$ matrices has a close algebraic analogy with the Dirac equation, enabling the use of the rich machinery of the Dirac electron theory. The beam optical Hamiltonian derived from this representation reproduces the Hamiltonians obtained in the traditional prescription along with wavelength-dependent matrix terms, which we have named as the polarization terms. These polarization terms are very similar to the spin terms in the Dirac electron theory and the spin-precession terms in the beam-optical version of the Thomas-BMT equation . The matrix formulation provides a unified treatment of beam optics and light polarization. Some well known results of light polarization are obtained as the paraxial limit of the matrix formulation -. The traditional beam optics is completely obtained from our approach in the limit of small wavelength, $`\text{}\lambda 0`$, which we call as the traditional limit of our formalisms. This is analogous to the classical limit obtained by taking $`\mathrm{}0`$, in the quantum prescriptions.

From the Hamiltonians in the Table we make the following observations: The classical/traditional Hamiltonians of particle/light optics are modified by wavelength-dependent contributions in the quantum/non-traditional prescriptions respectively. The algebraic forms of these modifications in each row is very similar. This should not come as a big surprise. The starting equations have one-to-one algebraic correspondence: Helmholtz $``$ Klein-Gordon; Matrix form of Maxwell $``$ Dirac equation. Lastly, the de Broglie wavelength, $`\text{}\lambda _0`$, and $`\lambda `$ have an analogous status, and the classical/traditional limit is obtained by taking $`\text{}\lambda _00`$ and $`\text{}\lambda 0`$ respectively. The parallel of the analogies between the two systems is sure to provide us with more insights.

##

Appendix-F.
An Invitation to the Experimentalists

It would be worthwhile to experimentally look for the predicted image rotation and the wavelength-dependent modifications of the aberration coefficients.

##

Table A.
Hamiltonians in Different Prescriptions

The following are the Hamiltonians, in the different prescriptions of light beam optics and charged-particle beam optics for magnetic systems. $`\widehat{H}_{0,p}`$ are the paraxial Hamiltonians, with lowest order wavelength-dependent contributions.

| Light Beam Optics | Charged-Particle Beam Optics |
| --- | --- |
| Fermat’s Principle $`=\left\{n^2\left(𝒓\right)𝒑_{}^2\right\}^{1/2}`$ | Maupertuis’ Principle $`=\left\{p_0^2𝝅_{}^2\right\}^{1/2}qA_z`$ |
| Non-Traditional Helmholtz $`\widehat{H}_{0,p}=n\left(𝒓\right)+\frac{1}{2n_0}\widehat{𝒑}_{}^2\frac{\mathrm{i}\text{}\lambda }{16n_0^3}[\widehat{𝒑}_{}^2,\frac{}{z}n\left(𝒓\right)]`$ | Klein-Gordon Formalism $`\widehat{H}_{0,p}=p_0qA_z+\frac{1}{2p_0}\widehat{𝝅}_{}^2+\frac{\mathrm{i}\mathrm{}}{16p_0^4}[\widehat{𝝅}_{}^2,\frac{}{z}\widehat{𝝅}_{}^2]`$ |
| Maxwell, Matrix $`\widehat{H}_{0,p}=n\left(𝒓\right)+\frac{1}{2n_0}\widehat{𝒑}_{}^2\mathrm{i}\text{}\lambda \beta 𝚺𝒖\text{}+\frac{1}{2n_0}\text{}\lambda ^2w^2\beta `$ | Dirac Formalism $`\widehat{H}_{0,p}=p_0qA_z+\frac{1}{2p_0}\widehat{𝝅}_{}^2\frac{\mathrm{}}{2p_0}\left\{\mu \gamma 𝚺_{}𝑩_{}+\left(q+\mu \right)\mathrm{\Sigma }_zB_z\right\}+\mathrm{i}\frac{\mathrm{}}{m_0c}ϵB_z`$ |

Notation

| $`\mathrm{Refractive}\mathrm{Index},n\left(𝒓\right)=c\sqrt{ϵ\left(𝒓\right)\mu \left(𝒓\right)}\mathrm{Resistance},h\left(𝒓\right)=\sqrt{\mu \left(𝒓\right)/ϵ\left(𝒓\right)}𝒖\left(𝒓\right)=\frac{1}{2n\left(𝒓\right)}\mathbf{}n\left(𝒓\right)𝒘\left(𝒓\right)=\frac{1}{2h\left(𝒓\right)}\mathbf{}h\left(𝒓\right)`$ $`𝚺`$ and $`\beta `$ are the Dirac matrices. | $`\widehat{𝝅}_{}=\widehat{𝒑}_{}q𝑨_{}\mu _a\mathrm{anomalous}\mathrm{magnetic}\mathrm{moment}.ϵ_a\mathrm{anomalous}\mathrm{electric}\mathrm{moment}.\mu =2m_0\mu _a/\mathrm{},ϵ=2m_0ϵ_a/\mathrm{}\gamma =E/m_0c^2`$ |
| --- | --- |