text:
# Contents

## 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; Then one makes 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 -.

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. 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 Feschbach-Villars technique for linearizing the Klein-Gordon equation is summarized in Appendix-A. The Foldy-Wouthuysen transformation technique is outlined in Appendix-B.

The other application, B. Axially Symmetric Graded Index Medium is used to demonstrate the power of the formalism. The traditional approaches give six aberrations. Our formalism modifies these six aberration coefficients by wavelength-dependent contributions.

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 in the limit $`\lambda _00`$ where $`\lambda _0=\mathrm{}/p_0`$ is the de Broglie wavelength and $`p_0`$ is the design momentum of the system under study.

## 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 the refractive index, $`n(𝒓)=n(x,y,z)`$. 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|p_z1`$ and $`\left|n(𝒓)n_0\right|n_0`$. The refractive index is the order of unity. Let us further assume that the refractive index varies smoothly around the constant background value $`n_0`$ without any abrupt jumps or discontinuities. 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.

Here, we follow a procedure similar to the one used for linearizing the Klein-Gordon equation via the Feshbach-Villars linearizing procedure . The resulting Feshbach-Villars-like form has an algebraic structure very similar to the Dirac equation. This enables us to make an expansion using the Foldy-Wouthuysen transformation technique well-known in the Dirac electron theory . The resulting expansion reproduces the above expansion in (3) as it should. Furthermore it gives rise to a set of wavelength-dependent contributions. The formalism presented here is an elaboration of the recent work which provides an alternative to the traditional square-root technique of obtaining the optical Hamiltonian .

Let us start with the wave-equation in the rectilinear coordinate system.

$`\left\{\mathbf{}^2{\displaystyle \frac{n^2(𝒓)}{v^2}}{\displaystyle \frac{^2}{t^2}}\right\}\mathrm{\Psi }=0.`$ (4)
Let

$`\mathrm{\Psi }=\psi (𝒓)e^{\mathrm{i}\omega t},\omega >0,`$ (5)
then

$`\left\{\mathbf{}^2+{\displaystyle \frac{n^2(𝒓)}{v^2}}\omega ^2\right\}\psi (𝒓)=0.`$ (6)
At this stage we introduce the wavization,

$`\mathrm{i}\text{}\lambda \mathbf{}_{}\widehat{𝒑}_{},\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}p_z,`$ (7)
where $`\lambda `$ is the reduced wavelength given by $`\text{}\lambda =\lambda /2\pi `$, $`c=\text{}\lambda \omega `$ and $`n(𝒓)=c/v(𝒓)`$. It is to be noted that $`pqqp=\mathrm{i}\text{}\lambda `$. This is similar to the commutation relation 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 formalism is completely obtained from our formalism in the limit $`\text{}\lambda 0`$. Then, we get,

$`\left\{\left(\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\right)^2+\left(\widehat{𝒑}_{}^2n^2(𝒓)\right)\right\}\psi (𝒓)=0.`$ (8)
Next, we linearize Eq. (8) following a procedure similar to, the one which gives the Feshbach-Villars  form of the Klein-Gordon equation. To this end, let

$`\left(\begin{array}{c}\psi _1(𝒓)\\ \psi _2(𝒓)\end{array}\right)=\left(\begin{array}{c}\psi (𝒓)\\ \mathrm{i}\frac{\text{}\lambda }{n_0}\frac{}{z}\psi (𝒓)\end{array}\right).`$ (13)
Then, Eq. (8) is equivalent to

$`\mathrm{i}{\displaystyle \frac{\text{}\lambda }{n_0}}{\displaystyle \frac{}{z}}\left(\begin{array}{c}\psi _1(𝒓)\\ \psi _2(𝒓)\end{array}\right)=\left(\begin{array}{ccc}0& & 1\\ \frac{1}{n_0^2}\left(n^2(𝒓)\widehat{𝒑}_{}^2\right)& & 0\end{array}\right)\left(\begin{array}{c}\psi _1(𝒓)\\ \psi _2(𝒓)\end{array}\right)`$ (20)
Next, we make the transformation,

$`\left(\begin{array}{c}\psi _1\left(𝒓\right)\\ \psi _2\left(𝒓\right)\end{array}\right)\mathrm{\Psi }^{(1)}`$ $`=`$ $`\left(\begin{array}{c}\psi _+\left(𝒓\right)\\ \psi _{}\left(𝒓\right)\end{array}\right)=M\left(\begin{array}{c}\psi _1(𝒓)\\ \psi _2(𝒓)\end{array}\right)`$ (27)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}\psi _1(𝒓)+\psi _2(𝒓)\\ \psi _1(𝒓)\psi _2(𝒓)\end{array}\right)`$ (30)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}\psi (𝒓)\mathrm{i}\frac{\text{}\lambda }{n_0}\frac{}{z}\psi (𝒓)\\ \psi (𝒓)+\mathrm{i}\frac{\text{}\lambda }{n_0}\frac{}{z}\psi (𝒓)\end{array}\right)`$ (33)
where

$`M=M^1={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}1\hfill & \hfill 1\\ 1\hfill & \hfill 1\end{array}\right),detM=1.`$ (36)
It is to be noted that the transformation matrix $`M`$ is independent of $`z`$. For a monochromatic quasiparaxial beam (in forward direction), with leading $`z`$-dependence $`\psi (𝒓)\mathrm{exp}\{\mathrm{i}n(𝒓)z/\text{}\lambda \}`$. Then

$`\psi _+`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left\{1+{\displaystyle \frac{n(𝒓)}{n_0}}\right\}\psi (𝒓)`$
$`\psi _{}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left\{1{\displaystyle \frac{n(𝒓)}{n_0}}\right\}\psi (𝒓)`$ (37)
Since, $`\left|n(𝒓)n_0\right|n_0`$, we have $`\psi _+\psi _{}`$.

Consequently, Eq. (8) can be written as

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\left(\begin{array}{c}\psi _+(𝒓)\\ \psi _{}(𝒓)\end{array}\right)`$ $`=`$ $`\widehat{\mathrm{H}}\left(\begin{array}{c}\psi _+(𝒓)\\ \psi _{}(𝒓)\end{array}\right),`$ (42)
$`\widehat{\mathrm{H}}`$ $`=`$ $`n_0\sigma _z+\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`{\displaystyle \frac{1}{2n_0}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2\left(𝒓\right)\right)\right\}\sigma _z`$
$`\widehat{𝒪}`$ $`=`$ $`{\displaystyle \frac{1}{2n_0}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}\left(\mathrm{i}\sigma _y\right),`$ (43)
where $`\sigma _y`$ and $`\sigma _z`$ are, respectively, the $`y`$ and $`z`$ components of the triplet of 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).$$
(44)
It is to be noted that the even-part and odd-part in Hamiltonian (43) differ only by a Pauli matrix. This simplifies the commutations a lot as we shall see, shortly. The details of the Feshbach-Villars linearizing procedure for the Klein-Gordon equation are available in Appendix-A.

The square of the Hamiltonian is

$`\widehat{\mathrm{H}}^2`$ $`=`$ $`\left\{\left(n^2\left(𝒓\right)\widehat{𝒑}_{}^2\right)\right\},`$ (45)
as expected. Thus we have have taken the square-root is a different way. This has certain distinct advantages over the traditional procedure of directly taking the square-root.

The purpose of casting Eq. (8) in the form of Eq. (43) will be obvious now, when we compare the latter with the form of the Dirac equation

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(\begin{array}{c}\mathrm{\Psi }_u\\ \mathrm{\Psi }_l\end{array}\right)`$ $`=`$ $`\widehat{H}_D\left(\begin{array}{c}\mathrm{\Psi }_u\\ \mathrm{\Psi }_l\end{array}\right)`$ (50)
$`\widehat{H}_D`$ $`=`$ $`m_0c^2\beta +\widehat{}_D+\widehat{𝒪}_D`$
$`\widehat{}_D`$ $`=`$ $`q\varphi `$
$`\widehat{𝒪}_D`$ $`=`$ $`c𝜶\widehat{𝝅},`$ (51)
where $`u`$ and $`l`$ stand for the upper and lower components respectively and

$`𝜶`$ $`=`$ $`\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].`$ (58)
To proceed further, we note the striking similarities between Eq. (43) and Eq. (51). In the nonrelativistic positive energy case, the upper components $`\mathrm{\Psi }_u`$ are large compared to the lower components $`\mathrm{\Psi }_l`$. The odd ($`\widehat{𝒪}`$) part of ($`\widehat{H}_Dm_0c^2\beta `$), anticommuting with $`\beta `$ couples the large $`\mathrm{\Psi }_u`$ to $`\mathrm{\Psi }_l`$ while the even ($`\widehat{}`$) part commuting with $`\beta `$, does not couple them. Using this fact, the well known Foldy-Wouthuysen formalism of the Dirac electron theory (see, e.g.) employs a series of transformations on Eq. (51) to reach a representation in which the Hamiltonian is a sum of the nonrelativistic part and a series of relativistic correction terms; $`\left|c\widehat{𝝅}\right|/m_0c^2`$ serves as the expansion parameter and the nonrelativistic part corresponds to an approximation of order up to $`\left|c\widehat{𝝅}\right|/m_0c^2`$. The terms of higher order in $`\left|c\widehat{𝝅}\right|/m_0c^2`$ constitute the relativistic corrections. Examining Eq. (43) we conclude $`\psi _+\psi _{}`$, and the odd operator $`\widehat{𝒪}`$, anticommuting with $`\sigma _z`$, couples the large $`\psi _+`$ with the small $`\psi _{}`$, while the even operator $`\widehat{}`$ does not make such a coupling. This spontaneously suggests that a Foldy-Wouthuysen-like technique can be used to transform Eq. (43) into a representation in which the corresponding beam optical Hamiltonian is a series with expansion parameter $`\left|\widehat{𝒑}_{}\right|/n_0`$. The correspondence between the beam optical Hamiltonian (43) and the Dirac electron theory is summarized in the following table:

The Analogy

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

Application of the Foldy-Wouthuysen-like technique to Eq. (43) involves a series of transformations on it and after the required number of transformations, depending on the degree of accuracy, Eq. (43) is transformed into a form in which the residual odd part can be neglected and hence the upper and lower components ($`\psi _+`$ and $`\psi _{}`$) are effectively decoupled. In this representation the larger component ($`\psi _+`$) corresponds to the beam moving in the $`+z`$-direction and the smaller component ($`\psi _{}`$) corresponds to the backward moving component of the beam.

Using the correspondence between Eq. (43) and Eq. (51) the Foldy-Wouthuysen expansion given formally in terms of $`\widehat{}`$ and $`\widehat{𝒪}`$ leads to the Hamiltonian

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}|\psi `$ $`=`$ $`\widehat{}^{(2)}|\psi ,`$
$`\widehat{}^{(2)}`$ $`=`$ $`n_0\sigma _z+\widehat{}{\displaystyle \frac{1}{2n_0}}\sigma _z\widehat{𝒪}^2,`$ (59)
To simplify the formal Hamiltonian we use, $`\widehat{𝒪}^2=\frac{1}{4n_0^2}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}^2`$ and recall that $`\widehat{}=\frac{1}{2n_0}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}\sigma _z`$. Dropping the $`\sigma _z`$ the formal Hamiltonian in in (59) is expressed in terms of the phase-space variables as:

$`\widehat{}^{(2)}`$ $`=`$ $`n_0+{\displaystyle \frac{1}{2n_0}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}`$ (60)
$`+{\displaystyle \frac{1}{8n_0^3}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}^2.`$
The details of the Foldy-Wouthuysen iterative procedure are described in detain in Appendix-B. The lowest order Hamiltonian obtained in this procedure agrees with the traditional approaches.

To go beyond the expansions in (60) one goes a step further in the Foldy-Wouthuysen iterative procedure. To next-to-leading order the Hamiltonian is formally given by

$`\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}|\psi `$ $`=`$ $`\widehat{}^{(4)}|\psi ,`$
$`\widehat{}^{(4)}`$ $`=`$ $`n_0\sigma _z+\widehat{}{\displaystyle \frac{1}{2n_0}}\sigma _z\widehat{𝒪}^2`$ (61)
$`{\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}}\sigma _z\left\{\widehat{𝒪}^4+\left([\widehat{𝒪},\widehat{}]+\mathrm{i}\text{}\lambda {\displaystyle \frac{}{z}}\widehat{𝒪}\right)^2\right\}.`$
As before we drop the $`\sigma _z`$ and the resulting Hamiltonian in the phase-space variable is

$`\widehat{}^{(4)}`$ $`=`$ $`n_0+{\displaystyle \frac{1}{2n_0}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}`$ (62)
$`+{\displaystyle \frac{1}{8n_0^3}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}^2`$
$`{\displaystyle \frac{\mathrm{i}\text{}\lambda }{32n_0^4}}[\widehat{𝒑}_{}^2,{\displaystyle \frac{}{z}}\left(n^2(𝒓)\right)]`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}\left({\displaystyle \frac{}{z}}\left(n^2(𝒓)\right)\right)^2`$
$`+{\displaystyle \frac{1}{16n_0^5}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}^3`$
$`+{\displaystyle \frac{5}{128n_0^7}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n^2(𝒓)\right)\right\}^4.`$
The Hamiltonian thus derived has all the terms which one gets in the traditional square-root approaches. In addition we also get the wavelength-dependent contributions.

The details of the various transforms and the beam optical formalism being discussed here turns out to be a simplified analog of the more general formalism recently developed for the quantum theory of charged-particle beam optics -, both in the scalar and the spinor cases, respectively. A very detailed description of these transforms and techniques is available in .

Now, we can compare the above Hamiltonians with the conventional Hamiltonian given by the square-root approach . The square-root approach does not give all the terms, such as the one involving the commutator of $`𝒑_{}^2`$ with $`\frac{}{z}\left(n^2(𝒓)\right)`$. Our procedure of linearization and expansion in powers of $`\left|\widehat{𝒑}_{}\right|/n_0`$ gives all the terms which one gets by the square-root expansion of (3) and some additional terms, which are the wavelength-dependent terms. Such, wavelength-dependent terms can in no way be obtained by any of the conventional prescriptions, starting with the Helmholtz equation (6).

## 3 Applications

In the previous sections we presented an alternative to the square-root expansion and 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 commutator term, which does not vanish for a varying refractive index.

Now, we apply the formalism to specific examples. First 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.

### 3.1 Medium with Constant Refractive Index

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

$`\widehat{\mathrm{H}}_c`$ $`=`$ $`n_0\sigma _z+𝒟\sigma _z+𝒟\left(\mathrm{i}\sigma _y\right)`$
$`𝒟`$ $`=`$ $`{\displaystyle \frac{1}{2n_0}}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n_c^2\right)\right\}.`$ (63)
The Hamiltonian in (63) can be exactly diagonalized by the following transform,

$`T^\pm `$ $`=`$ $`\mathrm{exp}\left[\mathrm{i}\left(\pm \mathrm{i}\sigma _z\right)\widehat{𝒪}\theta \right]`$ (64)
$`=`$ $`\mathrm{exp}\left[\sigma _x𝒟\theta \right]`$
$`=`$ $`\mathrm{cosh}\left(𝒟\theta \right)\sigma _x\mathrm{sinh}\left(𝒟\theta \right).`$
We choose,

$`\mathrm{tanh}\left(2𝒟\theta \right)={\displaystyle \frac{𝒟}{n_0𝒟}}={\displaystyle \frac{n_0^2\left(n_c^2\widehat{𝒑}_{}^2\right)}{n_0^2+\left(n_c^2+\widehat{𝒑}_{}^2\right)}}<1,`$ (65)
then we obtain,

$`\widehat{\mathrm{H}}_c^{\mathrm{diagonal}}`$ $`=`$ $`T^+\widehat{\mathrm{H}}_cT^{}`$ (66)
$`=`$ $`T^+\left\{n_0\sigma _z+𝒟\sigma _z+𝒟\left(\mathrm{i}\sigma _y\right)\right\}T^{}`$
$`=`$ $`\sigma _z\left\{n_0^22n_0𝒟\right\}^{\frac{1}{2}}`$
$`=`$ $`\sigma _z\left\{n_c^2\widehat{𝒑}_{}^2\right\}^{\frac{1}{2}}`$
We next, compare the exact result thus obtained with the approximate one, obtained through the systematic series procedure we have developed. We define $`P=\frac{1}{n_0^2}\left\{\widehat{𝒑}_{}^2+\left(n_0^2n_c^2\right)\right\}`$. Then,

$`\widehat{}_c^{(4)}`$ $`=`$ $`n_0\left\{1{\displaystyle \frac{1}{2}}P{\displaystyle \frac{1}{8}}P^2{\displaystyle \frac{1}{16}}P^3{\displaystyle \frac{5}{128}}P^4\right\}\sigma _z`$ (67)
$``$ $`n_0\left\{1P^2\right\}^{\frac{1}{2}}`$
$`=`$ $`\left\{n_c^2\widehat{𝒑}_{}^2\right\}^{\frac{1}{2}}`$
$`=`$ $`\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^{}),`$ (68)
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{},`$ (69)
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],$$
(70)
with

$`\widehat{T}(z^{\prime \prime },z^{})`$ $`=`$ $`{\displaystyle _z^{}^{z^{\prime \prime }}}𝑑z\widehat{}(z)`$ (71)
$`+{\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-C. 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],`$
$`\mathrm{\Delta }z=z_{\mathrm{out}}z_{\mathrm{in}},`$ (72)
Using (72), 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}}.`$ (79)
The beam-optical Hamiltonian is intrinsically aberrating. Even for simplest situation of a constant refractive index, we have aberrations to all orders!

### 3.2 Axially Symmetric Graded Index Medium

The refractive index of an axially symmetric graded-index material can be most generally described by the following polynomial (see, pp. 117 in )

$`n(𝒓)=n_0+\alpha _2(z)𝒓_{}^2+\alpha _4(z)𝒓_{}^4+\alpha _6(z)𝒓_{}^6+\alpha _8(z)𝒓_{}^8+\mathrm{},`$ (80)
where, we have assumed the axis of symmetry to coincide with the optic-axis, namely the $`z`$-axis without any loss of generality. To write the beam-optical Hamiltonians we introduce the following notation

$`\widehat{T}`$ $`=`$ $`\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)`$
$`w_1(z)`$ $`=`$ $`{\displaystyle \frac{d}{dz}}\left\{2n_0\alpha _2(z)\right\}`$
$`w_2(z)`$ $`=`$ $`{\displaystyle \frac{d}{dz}}\left\{\alpha _2^2(z)+2n_0\alpha _4(z)\right\}`$
$`w_3(z)`$ $`=`$ $`{\displaystyle \frac{d}{dz}}\left\{2n_0\alpha _6(z)+2\alpha _2(z)\alpha _4(z)\right\}`$
$`w_4(z)`$ $`=`$ $`{\displaystyle \frac{d}{dz}}\left\{\alpha _4^2(z)+2\alpha _2(z)\alpha _6(z)+2n_0\alpha _8(z)\right\}`$ (81)
We also use, $`[A,B]_+=\left(AB+BA\right)`$. The beam-optical Hamiltonian is

$`\widehat{}`$ $`=`$ $`\widehat{H}_{0,p}+\widehat{H}_{0,(4)}+\widehat{H}_{0,(6)}+\widehat{H}_{0,(8)}`$
$`+\widehat{H}_{0,(2)}^{(\text{}\lambda )}+\widehat{H}_{0,(4)}^{(\text{}\lambda )}+\widehat{H}_{0,(6)}^{(\text{}\lambda )}+\widehat{H}_{0,(8)}^{(\text{}\lambda )}`$
$`\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(\widehat{𝒑}_{}^2𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^2\right)\alpha _4(z)𝒓_{}^4`$
$`\widehat{H}_{0,(6)}`$ $`=`$ $`{\displaystyle \frac{1}{16n_0^5}}\widehat{𝒑}_{}^6{\displaystyle \frac{\alpha _2(z)}{8n_0^4}}\left\{\left(\widehat{𝒑}_{}^4𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^4\right)+\widehat{𝒑}_{}^2𝒓_{}^2\widehat{𝒑}_{}^2\right\}`$
$`+{\displaystyle \frac{1}{8n_0^3}}\left\{\left(\alpha _2^2(z)2n_0\alpha _4(z)\right)\left(\widehat{𝒑}_{}^2𝒓_{}^4+𝒓_{}^4\widehat{𝒑}_{}^2\right)+2\alpha _2^2(z)𝒓_{}^2\widehat{𝒑}_{}^2𝒓_{}^2\right\}`$
$`\alpha _6(z)𝒓_{}^6`$
$`\widehat{H}_{0,(8)}`$ $`=`$ $`{\displaystyle \frac{5}{128n_0^7}}\widehat{𝒑}_{}^8{\displaystyle \frac{5\alpha _2(z)}{64n_0^6}}[\widehat{𝒑}_{}^4,\left[\widehat{𝒑}_{}^2𝒓_{}^2\right]_+]_+`$
$`+{\displaystyle \frac{1}{32n_0^5}}\{(3\alpha _2^2(z)4n_0\alpha _4(z))[\widehat{𝒑}_{}^4,𝒓_{}^4]_++5\alpha _2^2(z)[\widehat{𝒑}_{}^2,𝒓_{}^2]_+^2`$
$`(2\alpha _2^2(z)+4n_0\alpha _4(z))\widehat{𝒑}_{}^2𝒓_{}^4\widehat{𝒑}_{}^2\}`$
$`+{\displaystyle \frac{1}{16n_0^4}}\{4(\alpha _2^3(z)+n_0\alpha _2(z)\alpha _4(z)+n_0^2\alpha _6(z))[\widehat{𝒑}_{}^2,𝒓_{}^6]_+`$
$`5\alpha _2^3(z)[𝒓_{}^4,[\widehat{𝒑}_{}^2,𝒓_{}^2]_+]_+`$
$`+(2\alpha _2^3(z)+4n_0\alpha _2(z)\alpha _4(z))[𝒓_{}^2,𝒓_{}^2\widehat{𝒑}_{}^2𝒓_{}^2]_+\}`$
$`\alpha _8(z)𝒓_{}^8`$
$`\widehat{H}_{0,(2)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda ^2}{16n_0^4}}\left\{{\displaystyle \frac{d}{dz}}\left(n_0\alpha _2(z)\right)\right\}\widehat{T}`$
$`\widehat{H}_{0,(4)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda ^2}{32n_0^4}}w_2(z)\left(𝒓_{}^2\widehat{T}+\widehat{T}𝒓_{}^2\right)+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(z)𝒓_{}^4`$
$`\widehat{H}_{0,(6)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{3\text{}\lambda ^2}{32n_0^4}}w_3(z)\left(𝒓_{}^4\widehat{T}+\widehat{T}𝒓_{}^4\right)+{\displaystyle \frac{\text{}\lambda ^2}{16n_0^5}}w_1(z)w_2(z)𝒓_{}^6`$
$`\widehat{H}_{0,(8)}^{(\text{}\lambda )}`$ $`=`$ $`{\displaystyle \frac{\text{}\lambda ^2}{8n_0^4}}w_4(z)\left(𝒓_{}^6\widehat{T}+\widehat{T}𝒓_{}^6\right)+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}\left\{w_2^2(z)+2w_1(z)w_3(z)\right\}𝒓_{}^8`$ (82)
The reason for partitioning $`\widehat{}`$ in the above manner will be clear as we proceed.

The paraxial transfer maps are formally 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}},`$ (89)
where $`P`$, $`Q`$, $`R`$ and $`S`$ are the solutions of the paraxial Hamiltonian in (82).

The transfer operator is most accurately expressed in terms of the the paraxial solutions, $`P`$, $`Q`$, $`R`$ and $`S`$, via the interaction picture .

$`\widehat{𝒯}(z,z_0)`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{\mathrm{i}}{\text{}\lambda }}\widehat{T}(z,z_0)\right],`$ (90)
$`=`$ $`\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)]_+`$
$`+A(z^{\prime \prime },z^{})\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)^2`$
$`+F(z^{\prime \prime },z^{})\left(\widehat{𝒑}_{}^2𝒓_{}^2+𝒓_{}^2\widehat{𝒑}_{}^2\right)`$
$`+D(z^{\prime \prime },z^{})[𝒓_{}^2,\left(\widehat{𝒑}_{}𝒓_{}+𝒓_{}\widehat{𝒑}_{}\right)]_+`$
$`+E(z^{\prime \prime },z^{})𝒓_{}^4\}],`$
The six 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}{8n_0^4}}w_2(z)Q^3S+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(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}{32n_0^4}}w_2(z)\left(Q^2(PS+QR)+2PQ^2S\right)`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(z)PQ^3\}`$
$`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}{16n_0^4}}w_2(z)\left(PQ(PS+QR)\right)`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(z)P^2Q^2\}`$
$`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}{16n_0^4}}w_2(z)\left(PQ(PS+QR)\right)`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(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}{32n_0^4}}w_2(z)\left(P^2(PS+QR)+2P^2QR\right)`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(z)P^3Q\}`$
$`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`$ (91)
$`{\displaystyle \frac{\text{}\lambda ^2}{8n_0^4}}w_2(z)\left(P^3R\right)`$
$`+{\displaystyle \frac{\text{}\lambda ^2}{32n_0^5}}w_1^2(z)P^4\}.`$
Thus we see that the transfer operator and the aberration coefficients are modified by $`\lambda `$-dependent contributions.

The sixth and eighth order Hamiltonians are modified by the presence of wavelength-dependent terms. These will in turn modify the the fifth and seventh order aberrations respectively -.

## 4 Concluding Remarks

We exploited the similarities between the Helmholtz equation and the Klein-Gordon equation to obtain an alternate prescription for the aberration expansion. In this prescription we followed a procedure due to Feschbach-Villars for linearizing the Klein-Gordon equation. After casting the Helmholtz equation to this linear form, it was further possible to use the Foldy-Wouthuysen transformation technique of the Dirac electron theory. This enabled us to obtain the beam-optical Hamiltonian to any desired degree of accuracy. We further get the wavelength-dependent contributions to 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 refractive index.

As an example we considered 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 was primarily for illustrating certain aspects of the machinery we have used.

The second, and the more interesting example is that of the axially symmetric graded index medium. For this system we derived the beam-optical Hamiltonians to eighth order. At each order we find the wavelength-dependent contributions. The fourth order Hamiltonian was used to obtain the six, third order aberrations coefficients which get modified by the wavelength-dependent contributions. Explicit relations for these coefficients were presented. In the limit $`\text{}\lambda 0`$, the alternate prescription here, reproduces the very well known Lie Algebraic Formalism of Light Optics. It would be worthwhile to look for the extra wavelength-dependent contributions experimentally.

The close analogy between geometrical optics and charged-particle has been known for too long a time. Until recently it was possible to see this analogy only between the geometrical optics and 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, 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 . We shall examine the parallel of these analogies in Appendix-D and summarize the Hamiltonians in the various prescriptions in Table-A.

##

Appendix A.
The Feshbach-Villars Form of the Klein-Gordon Equation

The method we have followed to cast the time-independent Klein-Gordon equation into a beam optical form linear in $`\frac{}{z}`$, suitable for a systematic study, through successive approximations, using the Foldy-Wouthuysen-like transformation technique borrowed from the Dirac theory, is similar to the way the time-dependent Klein-Gordon equation is transformed (Feshbach and Villars, ) to the Schrö-dinger form, containing only first-order time derivative, in order to study its nonrelativistic limit using the Foldy-Wouthuysen technique (see, e.g., Bjorken and Drell, ).

Defining

$$\mathrm{\Phi }=\frac{}{t}\mathrm{\Psi },$$
(A.1)
the free particle Klein-Gordon equation is written as

$$\frac{}{t}\mathrm{\Phi }=\left(c^2^2\frac{m_{0}^{}{}_{}{}^{2}c^4}{\mathrm{}^2}\right)\mathrm{\Psi }.$$
(A.2)
Introducing the linear combinations

$$\mathrm{\Psi }_+=\frac{1}{2}\left(\mathrm{\Psi }+\frac{\mathrm{i}\mathrm{}}{m_0c^2}\mathrm{\Phi }\right),\mathrm{\Psi }_{}=\frac{1}{2}\left(\mathrm{\Psi }\frac{\mathrm{i}\mathrm{}}{m_0c^2}\mathrm{\Phi }\right)$$
(A.3)
the Klein-Gordon equation is seen to be equivalent to a pair of of coupled differential equations:

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }_+`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2^2}{2m_0}}\left(\mathrm{\Psi }_++\mathrm{\Psi }_{}\right)+m_0c^2\mathrm{\Psi }_+`$
$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }_{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2^2}{2m_0}}\left(\mathrm{\Psi }_++\mathrm{\Psi }_{}\right)m_0c^2\mathrm{\Psi }_{}.`$ (A.4)
Equation (A.4) can be written in a two-component language as

$$\mathrm{i}\mathrm{}\frac{}{t}\left(\begin{array}{c}\mathrm{\Psi }_+\\ \mathrm{\Psi }_{}\end{array}\right)=\widehat{H}_0^{FV}\left(\begin{array}{c}\mathrm{\Psi }_+\\ \mathrm{\Psi }_{}\end{array}\right),$$
(A.5)
with the Feshbach-Villars Hamiltonian for the free particle, $`\widehat{H}_0^{FV}`$, given by

$`\widehat{H}_0^{FV}`$ $`=`$ $`\left(\begin{array}{cc}m_0c^2+\frac{\widehat{p}^2}{2m_0}& \frac{\widehat{p}^2}{2m_0}\\ \frac{\widehat{p}^2}{2m_0}& m_0c^2\frac{\widehat{p}^2}{2m_0}\end{array}\right)`$ (A.8)
$`=`$ $`m_0c^2\sigma _z+{\displaystyle \frac{\widehat{p}^2}{2m_0}}\sigma _z+\mathrm{i}{\displaystyle \frac{\widehat{p}^2}{2m_0}}\sigma _y.`$ (A.9)
For a free nonrelativistic particle with kinetic energy $`m_0c^2`$, it is seen that $`\mathrm{\Psi }_+`$ is large compared to $`\mathrm{\Psi }_{}`$.

In presence of an electromagnetic field, the interaction is introduced through the minimal coupling

$$\widehat{𝒑}\widehat{𝝅}=\widehat{𝒑}q𝑨,\mathrm{i}\mathrm{}\frac{}{t}\mathrm{i}\mathrm{}\frac{}{t}q\varphi .$$
(A.10)
The corresponding Feshbach-Villars form of the Klein-Gordon equation becomes

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\left(\begin{array}{c}\mathrm{\Psi }_+\\ \mathrm{\Psi }_{}\end{array}\right)`$ $`=`$ $`\widehat{H}^{FV}\left(\begin{array}{c}\mathrm{\Psi }_+\\ \mathrm{\Psi }_{}\end{array}\right)`$ (A.15)
$`\left(\begin{array}{c}\mathrm{\Psi }_+\\ \\ \mathrm{\Psi }_{}\end{array}\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\begin{array}{c}\mathrm{\Psi }+\frac{1}{m_0c^2}\left(\mathrm{i}\mathrm{}\frac{}{t}q\varphi \right)\mathrm{\Psi }\\ \\ \mathrm{\Psi }\frac{1}{m_0c^2}\left(\mathrm{i}\mathrm{}\frac{}{t}q\varphi \right)\mathrm{\Psi }\end{array}\right)`$ (A.22)
$`\widehat{H}^{FV}`$ $`=`$ $`m_0c^2\sigma _z+\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`q\varphi +{\displaystyle \frac{\widehat{\pi }^2}{2m_0}}\sigma _z,\widehat{𝒪}=i{\displaystyle \frac{\widehat{\pi }^2}{2m_0}}\sigma _y.`$ (A.23)
As in the free-particle case, in the nonrelativistic situation $`\mathrm{\Psi }_+`$ is large compared to $`\mathrm{\Psi }_{}`$. The even term $`\widehat{}`$ does not couple $`\mathrm{\Psi }_+`$ and $`\mathrm{\Psi }_{}`$ whereas $`\widehat{𝒪}`$ is odd which couples $`\mathrm{\Psi }_+`$ and $`\mathrm{\Psi }_{}`$. Starting from (A.23), the nonrelativistic limit of the Klein-Gordon equation, with various correction terms, can be understood using the Foldy-Wouthuysen technique (see, e.g., Bjorken and Drell, ).

It is clear from the above that we have just adopted the above technique for studying the $`z`$-evolution of the Klein-Gordon wavefunction of a charged-particle beam in an optical system comprising a static electromagnetic field. The additional feature of our formalism is the extra approximation of dropping $`\sigma _z`$ in an intermediate stage to take into account the fact that we are interested only in the forward-propagating beam along the $`z`$-direction.

##

Appendix-B.
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). 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)`$ (B.1)
$`\widehat{H}_D`$ $`=`$ $`m_0c^2\beta +q\varphi +c𝜶\widehat{𝝅}`$
$`=`$ $`m_0c^2\beta +\widehat{}+\widehat{𝒪}`$
$`\widehat{}`$ $`=`$ $`q\varphi `$
$`\widehat{𝒪}`$ $`=`$ $`c𝜶\widehat{𝝅},`$ (B.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],`$ (B.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].`$ (B.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 .$$
(B.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}}.`$ (B.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 }^{}`$ (B.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}.$$
(B.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)`$ (B.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{}`$ (B.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}}]`$ (B.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 (B.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,`$ (B.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}}`$ (B.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{}_2`$ $``$ $`\widehat{}_1,\widehat{𝒪}_2{\displaystyle \frac{\beta }{2m_0c^2}}\left([\widehat{𝒪}_1,\widehat{}_1]+\mathrm{i}\mathrm{}{\displaystyle \frac{\widehat{𝒪}_1}{t}}\right),`$ (B.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},$$
(B.24)
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),`$ (B.25)
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`$ (B.26)
$`{\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,`$ (B.27)
and retaining only the relevant terms of desired order at each step.

With $`\widehat{}=q\varphi `$ and $`\widehat{𝒪}=c𝜶\widehat{𝝅}`$, the final reduced Hamiltonian (B.26) 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 𝚺𝑩`$ (B.28)
$`{\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\sigma _z+\widehat{}+\widehat{𝒪}`$ |
| $`m_0c^2`$ | $`n_0`$ |
| Positive Energy | Forward Propagation |
| Nonrelativistic, $`\left|\widehat{𝝅}\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 (B.1, B.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-C
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{}}.$$
(C.1)
Let it be required to solve the differential equation

$$\frac{}{t}u(t)=\widehat{A}(t)u(t)$$
(C.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).$$
(C.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)`$ (C.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).$$
(C.5)
Thus, $`u(T)`$ is given by computing the product in (C.5) using successively the BCH-formula (C.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)`$ (C.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 (C.6) is obtained let us substitute the assumed form of the solution, $`u(t)=\widehat{𝒯}(t,t_0)u\left(t_0\right)`$, in (C.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}.$$
(C.7)
Introducing an iteration parameter $`\lambda `$ in (C.7), let

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

$$\widehat{𝒯}(t,t_0;\lambda )=\mathrm{e}^{\mathrm{\Omega }(t,t_0;\lambda )}$$
(C.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.$$
(C.11)
Now, using the identity (see, 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 )},$$
(C.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).$$
(C.13)
Substituting in (A13) the series expression for $`\mathrm{\Omega }(t,t_0;\lambda )`$ (C.11), expanding the left hand side using the first identity in (C8), 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$$
(C.14)
and hence

$$\mathrm{\Delta }_1(t,t_0)=_{t_0}^t𝑑t_1\widehat{A}(t_1).$$
(C.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$$
(C.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)].$$
(C.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)]`$ (C.18)
$`+[[\widehat{A}(t_3),\widehat{A}(t_2)],\widehat{A}(t_1)]\}.`$
Then, the Magnus formula in (C.6) follows from (C.9)-(C.11). Equation 71 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-D
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. The 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 from the quantum theory of charged-particle beam optics, which are all 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 .

### D.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 and may also help get a deeper understanding and to lead to better designing of charged-particle beam devices.

The starting point to obtain a quantum prescription of charged particle beam optics is to build a theory based on the basic equations (Schrödinger, Klein-Gordon, Dirac) of quantum mechanics 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),$$
(D.1)
where $`(x,y;s)`$ constitute a curvilinear coordinate system, adapted to the geometry of the system. Eq. (D.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 (D.1)

$`\psi (x,y;s)`$ $`=`$ $`\widehat{U}(s,s_{\mathrm{in}})\psi (x,y;s_{\mathrm{in}}),`$ (D.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),`$ (D.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 analogy 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 ($`\overline{\lambda }_0=\mathrm{}/p_0`$). The classical limit of the quantum formalism reproduces the well known Lie algebraic formalism  of charged-particle beam optics.

### D.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 .

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 born 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 complete 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 are 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 algebraically 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 a 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, $`\overline{\lambda }_0`$, and $`\lambda `$ have an analogous status, and the classical/traditional limit is obtained by taking $`\overline{\lambda }_00`$ and $`\text{}\lambda 0`$ respectively. The parallel of the analogies between the two systems is sure to provide us with more insights.

##

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`$ |
| --- | --- |