Clarification on the principles of (matrix) geometrical optics as presented in the first pages of "Symplectic Techniques in Physics"

Question

I found the book Symplectic Techniques in Physics by V. Guillemin and S. Sternberg in my university library and started skimming it. The first chapter aims to guide the reader through the historical path from Hamilton's work in geometrical optics to the modern symplectic formulation of mechanics. Their presentation of the rudiments of optics is a bit sketchy (I would say, obscure)... Help!

They start by posing the Big Problem of the Gaussian approximation of geometrical optics. We have an optical system, which we consider symmetric around some axis (the optical axis, which will be our $ z $-axis), and our aim is (semi-quoting) "to trace the trajectory of a light ray as it passes through the various refracting surfaces of the optical system (or is reflected by the reflecting surfaces)". We will stick with trajectories lying in the same plane.

Assuming light travels in straight lines, we need a way of specifying straight lines. The way that the book presents as the most convenient is the following. To specify a line, we pick a coordinate $ z_0 $ on the optical axis (equivalently, we pick a reference plane perpendicular to the optical axis at some $ z_0 $), and we provide the height $ q $ of the line above the axis at $ z_0 $ and the angle $ \theta $ that the line makes with the optical axis (positive if a counterclockwise rotation clashes the positive $ z $-direction of the axis into the ray).

Our problem now becomes to determine the pair $ (q_2, \theta_2) $ that corresponds to the line emerging from the optical system, where $ q_2 $ is given with respect to some reference plane at a point $ z_2 $ placed after the optical system, to the pair $ (q_1, \theta_1) $ corresponding to the ray incoming into the system, where $ z_1 $ is chosen before the system.

Question 1. I tried to determine $ (q_2,\theta_2) $ from $ (q_1,\theta_1) $ explicitly in the concrete case of Fig. 26-6 of the first volume of the Feynman Lectures. Here a light beam is made to pass through a prism with plane parallel faces, set at an angle to the beam. The prism simply displaces the beam parallel to itself of a quantity $ \Delta h $ that can be computed explicitly. The point is that I don't get how phrase the problem in the $ (q,\theta) $ form. What would be a smart/canonical way of choosing the reference planes $ z_1 $, $ z_2 $ (and eventually some $ z_3 $)?

Let's go on and make the following substitution. If we are given a ray $ (q,\theta) $ with respect to some reference plane, we put $$ p = n\theta ,\tag{1}$$ where $ n $ is the refracting index at the reference plane. We thus describe a light ray as a pair $ (q,p) $, and our problem becomes the determination of $ (q_2,p_2) $ as a function of $ (q_1,p_1) $.

The authors claim that "from our approximations" (i.e., paraxial approximation) it follows that $$ \begin{pmatrix} p_2\\ q_2 \end{pmatrix} = M_{21} \begin{pmatrix} p_1\\ q_1 \end{pmatrix} \label{matrix}\tag{$\ast$} $$ where $ M_{21} $ is some matrix with determinant 1, i.e., $ (q_2,p_2) $ is a linear function of $ (q_1,p_1) $.

Question 2. Where does formula $ \eqref{matrix} $ comes from? I don't see how I could have derived it, probably because it's not clear to me how we are setting up the problem.

Answer 1

1. We start with the Hamiltonian action for optical length in 3D $$ \begin{align} S_H[{\bf r},{\bf p},e]~=~&\int_{\lambda_i}^{\lambda_f}\! d\lambda ~ L_H, \cr L_H~:=~&\underbrace{{\bf p}\cdot\dot{\bf r}}_{\text{presympl. pot.}} - \underbrace{\frac{e}{2}({\bf p}^2-n({\bf r})^2)}_{\text{Hamiltonian}},\end{align}\tag{A} $$ of 7 variables: 3-position ${\bf r}=(x,y,z)$, 3-momentum ${\bf p}=(p_x,p_y,p_z)$ and einbein $e>0$; cf. my Phys.SE answer here.

2. Let us write the 3-momentum ${\bf p}$ in spherical coordinates $(p_r,\theta,\varphi)$. If we eliminate/integrate out $e$ and $p_r$ in the Hamiltonian Lagrangian (A), we get $$ L_H \quad\stackrel{e,~p_r}{\longrightarrow}\quad n({{\bf r}})\frac{{\bf p}}{|{\bf p}|}\cdot\dot{\bf r} ~=~n({{\bf r}})\left( \sin\theta\cos\varphi \dot{x}+\sin\theta\sin\varphi \dot{y}+\cos\theta\dot{z}\right), \tag{B}$$ which depends on 5 variables $(x,y,z;\theta,\phi)$.

3. The actions (A) & (B) have a world-line reparametrization gauge symmetry. Let us for simplicity go to the axial gauge $z=\lambda$, i.e. along the optical $z$-axis. Then the Hamiltonian Lagrangian (B) becomes $$\left. L_H\right|_{z=\lambda} ~\stackrel{(B)}{=}~\underbrace{n({{\bf r}})\sin\theta\left(\cos\varphi \dot{x}+\sin\varphi \dot{y} \right)}_{\text{symplectic potential}} + \underbrace{ n({{\bf r}}) \cos\theta }_{\text{minus Hamiltonian}},\tag{C} $$ which depends on 4 phase space variables $(x,y;\theta,\phi)$ in a 4D phase space.

4. From the symplectic potential (C), we can identify the canonical momenta $$\begin{align} p_x~=~&n({{\bf r}})\sin\theta\cos\varphi, \cr p_y~=~&n({{\bf r}})\sin\theta\sin\varphi, \end{align}\tag{D} $$ cf. e.g. my Phys.SE answer here.

5. Now let us return to OP's questions. OP's formula (1) $$p_x~\stackrel{(D)}{\approx}~n({{\bf r}})\theta\tag{E}$$ is the paraxial approximation $\sin\theta\approx\theta\ll 1$ dimensionally reduced to $\varphi=0$.

6. OP's formula (*) corresponds to a linear canonical transformation. A symplectomorphism has determinant is 1, cf. my Math.SE answer here. The corresponding matrix is a ray transfer matrix. It may typical be decomposed in elementary optical operations.