  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.
