Proof of Validity of Thick Lens Model by Hamiltonian Formalism or Otherwise The paraxial imaging geometry of any axially symmetric optical system can be described by three real parameters: the axial positions of the two principal planes and one focal length (i.e. $f_f=f_b$ in the diagrams below - or in the case where the input and output refractive indices are different, the focal lengths are different but their ratio is the ratio between the refractive indices). A ray entering one principal plane "teleports" to the other to emerge at the same transverse position in the other.
This is an amazingly simple (i.e. only three parameters) model for any combination of axially symmetric optical devices. Is there a simple proof that it always works, maybe grounded on general, Hamiltonian optics or Fermat's principle? Or, are there cases which it is not adequate for?


 A: There are very recent devices discovered which this model may not be adequate for (more later), but the model is nonetheless very general. It works for any axially symmetric optical system made of the following subsystems fulfilling:


*

*Abrupt interfaces between different dielectric mediums (lens and mirror surfaces) are allowed;

*Aside from "jump discontinuities" of refractive index or its derivatives at interfaces between subsystems in lenses and mirrors in 1., inhomogeneous mediums must have refractive indices that are twice differentiable functions of position and Fermat's principle must govern ray propagation. In short, a Hamiltonian description with $\mathcal{H} = \frac{1}{2}\,g^{jk}\,p_j\,p_k$ must prevail in these regions.


Its simplicity arises from two symmetries that constrain away most of the general situation's degrees of freedom:


*

*The postulated rotational symmetry about an optical axis; together with

*The symplectic symmetry of a general passive optical system's transfer matrix, which follows from the validity of the Hamiltonian formulation aside from the discontinuities allowed above.


Each of these symmetries constrains away six of the sixteen possible degrees of freedom from the $4\times4$ ray transfer matrix, leaving four free parameters. One of these four parameters is a rotation angle (about the optical axis) that does not affect whether the system's imaging planes can be modelled by focal lengths and principal plane positions. Moreover, the fourth parameter can be nonzero only in anisotropic systems. The remaining three are equivalent to the diagrams shown by the OP.

Simple Demonstration
If we make the intuitively reasonable (but actually not fully correct)
 postulate that axial symmetry means that meridional rays are confined to meridional planes and thus that the system's behavior is fully modelled by the propagation of a general, meridional ray in any chosen meridional plane, then this means that our system can be modelled by a two dimensional optical phase space with $2\times2$ Ray Transfer Matrix. The validity of the Hamiltonian approach then guarantees that the ray transfer matrix has a unity determinant:
$$\left(\begin{array}{cc}a&b\\c&d\end{array}\right)\in\mathrm{Sp}(2,\,\mathbb{R})\cong \mathrm{SL}(2,\,\mathbb{R})\tag{8}$$
and this, with a few diagrams and a bit of algebra, can be readily shown to be logically equivalent to the diagrams. The diagrams imply such a transfer matrix, and a three independent parameter unimodular transfer matrix implies the diagrams. For instance, the axial positions of the principal planes $z_1$ and $z_2$ relative to the input and output and the two focal lengths are:  
$$\begin{array}{ll}\displaystyle{z_1 = f_b-d_b = \frac{d-1}{c}\,n_i}&\displaystyle{f_b = -\frac{n_i}{c}}\\\displaystyle{z_2 = d_f-f_f = \frac{1-a}{c}\,n_o}&\displaystyle{f_f = -\frac{n_o}{c}}\end{array}\tag{9}$$
where $n_i$ and $n_o$ are the refractive indices at the input and output.
Full Discussion: The Gory Details
Actually the assumption of confinement of meridional rays to meridional planes, although physically, intuitively reasonable and valid for the vast majority of passive optical systems, is not always correct: some anisotropic materials can violate it. So we have to account for the rotational symmetry more fully and rigorously. 
As I have outlined in recent answers here, here and here, the Hamiltonian formulation of the consequences of Fermat's least time principle shows that passive optical systems impart a symplectomorphism on the phase space state of the light passing through them. When we calculate image positions through a paraxial model, we linearize this symplectomorphism about a chief, or reference ray in phase space and calculate what happens to rays that are near to the chief ray in phase space. The states of these nearby rays are mapped by the Jacobi matrix of the symplectomorphism, and the Jacobi matrix $\Gamma\in\mathrm{Sp}(4,\,\mathbb{R})$ is always a symplectic matrix:
$$\Gamma^T\,\Omega\, \Gamma = \Omega;\quad\Omega = \left(\begin{array}{c|c}0&-\mathrm{id}\\\hline\mathrm{id}&0\end{array}\right)\tag{1}$$ 
where $\mathrm{id}$ is the $2\times2$ identity. These Jacobi matrices are the Ray Transfer Matrices and they act on the optical states of rays:
$$X = \left(\begin{array}{c}x\\y\\n\,\gamma_x\\n\,\gamma_y\end{array}\right)\tag{2}$$
where the $z$ axis is along the system's optical axis, $(x,\,y)$ is the position in the transverse plane of the ray (at the system input, and, for the image state, the position in the output plane), $(\gamma_x,\,\gamma_y)$ are the $x$ and $y$ components of the direction cosines of the ray's direction and $n$ is the refractive index at the reference ray's position in the transverse plane. The quantities $p_x=n\,\gamma_x$ and $p_y=n\,\gamma_y$ are the optical momentums conjugate (in the sense of Hamiltonian mechanics) to the positions $x$ and $y$; interestingly, they are indeed equivalent (modulo scaling by the constant $\hbar\,\omega/c$) to the $x$ and $y$ components of the photonic momentum $\hbar\,\vec{k}$, where $\vec{k}$ is the wavevector, but this  fact is an aside. (2) describes our points in optical phase space. 
An important aside: The fact that the Jacobi matrices of the transformations are symplectic means, in particular, that their determinant is unity so that the volume of the system of rays describing a light state is unchanged by the general, non-linearized, non-paraxial system symplectomorphism: this is Liouville's theorem applied to optics, where it is known as the law of conservation of étendue (Randall Munroe's Fire From Moonlight is also an obligatory read for any discussion of étendue).
From (1), the Lie algebra $\mathfrak{sp}(4,\,\mathbb{R})$ of $\mathrm{Sp}(4,\,\mathbb{R})$ is the following set of $4\times4$ matrices:
$$\mathfrak{sp}(4,\,\mathbb{R}) = \left\{G|\;G\text{ is }4\times 4;\,G^T\,\Omega + \Omega\,G=0\right\}\tag{3}$$
which can be shown to be equivalent to the following structure for $G$:
$$G = \left(\begin{array}{c|c}A&B\\\hline C&-A^T\end{array}\right);\quad B=B^T;\;C=C^T\tag{4}$$
where the partitions are all $2\times2$. This form shows that the general ray transfer matrix has "only" ten independent parameters, rather than sixteen. 
Now we impose our postulate of rotational symmetry. Any $\Gamma\in\mathrm{Sp}(4,\,\mathbb{R})$ is a finite product of the form $\Gamma = \exp(X_1)\,\exp(X_2),\,\cdots$ with $X_i\in\mathfrak{sp}(4,\,\mathbb{R})$. Rotational symmetry means that if we rotate our system through any angle about the optical axis, its transfer matrix will be unchanged. Therefore:
$$\begin{array}{cll}&\Gamma = \mathrm{R}\,\Gamma\,\mathrm{R}^{-1}&\mathrm{R}=\left(\begin{array}{c|c}\begin{array}{cc}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{array}&0\\\hline0&\begin{array}{cc}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{array}\end{array}\right)\\
\Leftrightarrow &\Gamma = \exp\left(\mathrm{R}\,X_2\,\mathrm{R}^{-1}\right)\exp\left(\mathrm{R}\,X_1\,\mathrm{R}^{-1}\right)\cdots&\\
\Leftrightarrow & [r,\,X_j]=0 & r=\log \mathrm{R} =\theta\,\left(\begin{array}{c|c}i&0\\\hline0&i\end{array}\right);\;\theta\in\mathbb{R} \\&&i=\left(\begin{array}{cc}0&-1\\1&0\end{array}\right)\end{array}\tag{5}$$
where $\mathrm{R}\in \mathrm{SO}(2)$ corresponds to a system rotation about the $z$ axis. 
If we combine the last line $[r,\,X_j]=0$ of (5) with (2), we see that the Lie algebra of our set of axially symmetric optical transfer matrices is the set of all matrices of the form:
$$\left(\begin{array}{c|c}\theta\,i + \tilde{a}\,\mathrm{id}&\tilde{b}\,\mathrm{id}\\\tilde{c}\,\mathrm{id}&\theta\,i-\tilde{a}\,\mathrm{id} \end{array}\right)\tag{6}$$
where $\mathrm{id}$ is, again, the $2\times2$ identity, $\theta,\,\tilde{a},\,\tilde{b},\,\tilde{c}\in\mathbb{R}$ and $i$ is the $2\times 2$ infinitessimal rotation as in (5). This Lie algebra is isomorphic to $\mathfrak{so}(2)\oplus\mathfrak{sp}(2,\,\mathbb{R})$; note that it has a central element $r$ (defined in (5)), so it has a one dimensional continuous center. Furthermore, the complement to this center is itself a Lie algebra, hence the algebra splits as a direct sum. Therefore, the most general possible transfer matrix is of the form:
$$\Gamma = R(\theta)\,\left(\begin{array}{c|c}a\,\mathrm{id}&b\,\mathrm{id}\\\hline c\,\mathrm{id}&d\,\mathrm{id}\end{array}\right)=\left(\begin{array}{c|c}a\,\mathrm{id}&b\,\mathrm{id}\\\hline c\,\mathrm{id}&d\,\mathrm{id}\end{array}\right)\,R(\theta);\quad \det(\Gamma)=a\,d-b\,c=1\tag{7}$$
What (7) says is that we can model any axially symmetric system as a rotation about the axis followed by a transformation that is equivalent to simply modelling the rays on one, meridional plane (e.g. the $x-z$ plane) by the usual, 2-dimensional $ABCD$ transfer matrix method. Intuitively, this is not surprising, but what perhaps is a little surprising is that our reasoning above shows that skew rays (those not confined to a meridional plane) don't change the analysis. Note that, for a general paraxial optical system with a general symplectic transfer matrix the skew rays are different, and rays in the meridional and sagittal planes have different focal lengths - this is the phenomenon of astigmatism. 
The rotation $R(\theta)$ does not affect our ability to describe the system by principal planes and focal lengths. It seems unphysical and indeed it doesn't show up in everyday lens systems, where rays confined to meridional planes stay on those planes. But it is a theoretical possibility: axially symmetric lenses made of an anisotropic material can in theory have their anisotropy arranged so that the arrangement is both axially symmetric and bears the rays along the system in helical paths, so that a rotation is combined with the usual thick lens model. I'm not aware of any real devices having been built to do this, however.
So now we are in a position to explain the diagrams. Given the above, a system's imaging plane is defined by the four parameters in the 2-dimensional ray transfer matrix describing the transfer of rays confined to a meridional plane, and we're back to the "Simple Demonstration" cited above.
Lastly, one should say something about recent exceptions to these ideas that have been studied. There are some recently discovered diffractive systems that do not impart a symplectomorphism on the state of the light passing through them. See, for example,  
Xiaomin Lai, Shaoqun Zeng, Xiaohua Lv, Jing Yuan and Ling Fu, "Violation of the Lagrange invariant in an optical imaging system“, Optics Letters, 38, #11,  pp1896-1898 2013
The optical invariant (aka Helmholtz invariant or Lagrange invariant) is the symplectic form on optical phase space. Therefore, the self interfering holography systems described in the paper above are not fully modelled by (4) and, in theory, one could build systems that had object position dependent principal plane positions and/ or focal lengths using such systems. These interesting systems do still conserve étendue, however; the Jacobi matrix of the transformation, although not symplectic, is still unimodular. 
