In the paraxial regime of Newtonian optics, propagation of an ensemble of rays is represented by a symplectic ABCD transfer matrix defined on a reduced phase space. Here, we present its analogue for general relativity.
Starting from simultaneously applied null geodesic actions for two curves, we obtain a geodesic deviation action up to quadratic order. We achieve this by following a preexisting method constructed via Synge's world function.
We find the corresponding Hamiltonian function and the reduced phase space coordinates that are composed of the components of the Jacobi fields projected on an observational screen. Our thin ray bundle transfer matrix is then obtained through the matrix representation of the Lie operator associated with this quadratic Hamiltonian.
Moreover, Etherington's distance reciprocity between any two points is shown to be equivalent to the symplecticity conditions of our ray bundle transfer matrix. We further interpret the bundle propagation as a free canonical transformation with a generating function that is equal to the geodesic deviation action.
We present it in the form of matrix inner products. A phase space distribution function and the associated Liouville equation is also provided.
Finally, we briefly sketch the potential applications of our construction. Those include reduced phase space and null bundle averaging; factorization of light propagation in any spacetime uniquely into its thin lens, pure magnifier and fractional Fourier transformer components; wavization of the ray bundle; reduced polarization optics and autonomization of the bundle propagation on the phase space to find its invariants and obtain the stability analysis.