The excerpt below is from Box 14.1 of Misner Thorne and Wheeler's Gravitation. This equation is presented in the quoted text:
$$\lim_{\epsilon\to0}\frac{6}{\epsilon^{2}}\left(1-\frac{\text{circumference}}{2\pi\epsilon}\right)=\kappa_{1}\kappa_{2}$$
The authors do not, however, derive the expression, and their verbal description of the procedure for deriving it is vague. How is it established? I'm adding the first paragraph of Box 14.1, since it might provide a clue as to how they derived the expression in question.
This is an illustration of how I see the situation. The surface is given by the binary quadratic form of the second degree Taylor series expansion about a critical point (MacLauren series?). The principle curves are axis aligned. The closed curve on the surface is the graph of the image of the unit circle under the mapping $q.$ That is $\left\{\cos \theta,\sin \theta,q[\cos \theta,\sin \theta]\right\}$. This is not the "circle" for which we seek the circumference. That set of points is determined by constant "radial" pathlength within the surface.
Since this is currently a "back-burner" question, I can't invest the time to produce a more accurate illustration in the near future. But it is also one of those problems that keeps me up at night. I feel like I "should" be able to solve it.
This is an illustration showing intrinsic "circles" on a quadratic surface. The blue curves are simply the projection of the red "circles" onto the $X\times{Y}$ plane.
Given the equation \begin{equation} z=\frac{1}{2} a x^2+b x y+\frac{1}{2} c y^2 \end{equation} It is possible to define the parametrization \begin{equation} \textbf{X}(x,y) =(x,y,\frac{1}{2} a x^2+b x y+\frac{1}{2} c y^2) \end{equation} Without loss of generality we can assume the point P on the surface is the origin of an orthonormal basis $(e_1,e_2,e_3)$, with $e_1$ and $e_2$ belonging to the tangent plane and $e_3$ the normal vector at P. We can then differentiate $\textbf{X}$ with respect to $x$ and $y$ and get the following expressions: \begin{equation} \frac{\partial \textbf{X}}{\partial x} =(1,0,a x+b y) \\ \frac{\partial \textbf{X}}{\partial y} =(0,1,b x+c y) \\ \end{equation} When computing in P we have \begin{equation} \frac{\partial \textbf{X}}{\partial x}(0,0) =(1,0,0)=e_1 \\ \frac{\partial \textbf{X}}{\partial y}(0,0) =(0,1,0)=e_2 \\ \end{equation} The first fundamental form in P is: \begin{equation} g_{11}(0,0)=1, g_{12}(0,0)=0, g_{22}(0,0)=1 \end{equation} Based on the above relation the quadratic expression $z$ can be written as \begin{equation} \textbf{X}(x,y) =(x,y,\frac{1}{2} L_{11}(0,0) x^2+L_{12}(0,0) x y+\frac{1}{2} L_{22}(0,0) y^2) \end{equation} Since $L_{ij}=\textbf{N} \cdot \textbf{X}_{ij}$ and the normal vector at P is $e_3=(0,0,1)$ we have
\begin{equation} L_{11}(0,0)= \frac{\partial^2 \textbf{X}}{\partial x^2}(0,0), L_{12}(0,0)=\frac{\partial^2 \textbf{X}}{\partial x \partial y}(0,0), L_{22}(0,0)=\frac{\partial^2 \textbf{X}}{\partial y^2}(0,0) \end{equation} And the second fundammental form in P is : \begin{equation} L_{11}(0,0)=a, L_{12}(0,0)=b, L_{22}(0,0)=c \end{equation} We can write $z$ as: \begin{equation} z =\frac{1}{2} \frac{\partial^2 \textbf{X}}{\partial x^2}(0,0) x^2+\frac{\partial^2 \textbf{X}}{\partial x \partial y}(0,0) x y+\frac{1}{2} \frac{\partial^2 \textbf{X}}{\partial y^2}(0,0) y^2 \end{equation} Since $z(0,0) = 0$, and the first partial derivative of $z$ at P are zero, the above expression represent the Taylor expansion of $2^{nd}$ degree of the function $z$. Considering the values of $L_{ij}$ previouly computed, $z$ can be written as a function of $\kappa_{1}$ and $\kappa_{2}$, and thus a second order approximation of the surface near P is: \begin{equation} (\xi, \eta, \frac{1}{2} \kappa_{1} \xi ^2+\frac{1}{2} \kappa_{2} \eta ^2) \end{equation} This is because for non umbelical points, it is possible to prove that there exists two orthogonal directions where the normal curvature has a maximum and minimum and $L_{11}=\kappa_{1}$, $L_{12}=0$, and $L_{22}=\kappa_{2}$. The product of $\kappa_{1}$ and $\kappa_{2}$ is equal to the ratio of the determinants of the second and first fundamental form. In our case the determinant of the first fundamental form in P is 1 so that \begin{equation} K=\kappa_{1} \kappa_{2} = \frac{det(L)}{det(g)}= det(L) = a c -b ^2 \end{equation} where $K$ is the Gaussian curvature. As you wrote in the text you added \begin{equation} ds^2=d\xi^2+/d\eta^2+(\kappa_{1} \xi d\xi+ \kappa_{2} \eta d \eta)^2 \end{equation} is the metric corresponding to the embedded viewpoint. Now we want to compute the Gaussian curvature $K$ i.e. the product of $\kappa_1$ and $\kappa_2$. To do that we consider a circle of radius $\epsilon$ around C. This circle is a function defined on $[0, \epsilon] \times [0, 2 \pi]$. If we call this function $\sqrt{g}(\epsilon,\theta)$ we must have: \begin{equation} \sqrt{g}(\epsilon,0)=\sqrt{g}(\epsilon,2 \pi) \end{equation} and \begin{equation} \sqrt{g}(0,\theta)=0 \end{equation} Now we can use the Theorema Egregium by Gauss who states that Gaussian curvature $K$ can be determined entirely by measuring angles, distances and their rates on a surface, i.e. can be determined from the first fundamental form, without requiring any knowledge on how the surface is embedded in the 3D Euclidean space.The theorema Egregium produces a formula for $K$ that in the case of orthogonal parametrization ($g_{12}=g_{21}=0$) of a neighborhood of a surface reduces to: \begin{equation} K=-\frac{1}{\sqrt{g_{11} g_{22}}}\left(\frac{\partial}{\partial x_1}\left(\frac{1}{\sqrt{g_{11}}} \frac{\partial \sqrt{g_{22}}}{\partial x_1}\right)+\frac{\partial}{\partial x_2}\left( \frac{1}{\sqrt{g_{22}}} \frac{\partial \sqrt{g_{11}}}{\partial x_2}\right)\right) \end{equation} If we consider the $x_1$ coordinate lines being parametrized by the arc length, the previous equation simplifies to: \begin{equation} K=-\frac{1}{\sqrt{g_{22}}} \frac{\partial^2 \sqrt{g_{22}}}{\partial {x_1}^2} \end{equation} We can now write: \begin{equation} \frac{\partial^2 \sqrt{g_{22}(\epsilon,\theta)}}{\partial {\epsilon}^2}=-\sqrt{g_{22}(\epsilon,\theta)} K(\epsilon,\theta) \end{equation} And \begin{equation} \frac{\partial^2 \sqrt{g_{22}(0,\theta)}}{\partial {\epsilon}^2}=-\lim_{\epsilon \to 0} \sqrt{g_{22}(\epsilon,\theta)}K(\epsilon,\theta)=0 \end{equation} As a matter of fact \begin{equation} \frac{\partial \sqrt{g_{22}(0,\theta)}}{\partial {\epsilon}}=\lim_{\epsilon \to 0}\sqrt{g_{22}(\epsilon,\theta)}=1 K(\epsilon,\theta) \end{equation} Continuing the process we can compute the third partial derivative as: \begin{equation} \frac{\partial^3 \sqrt{g_{22}(\epsilon,\theta)}}{\partial {\epsilon}^3}=-\frac{\partial \sqrt{g_{22}(\epsilon,\theta)}}{\partial {\epsilon}} K(\epsilon,\theta)-\sqrt{g_{22}(\epsilon,\theta)} \frac{\partial K(\epsilon,\theta)}{\partial {\epsilon}} \end{equation} In the limit where $\epsilon \to 0$ we can write \begin{equation} \frac{\partial^3 \sqrt{g_{22}(0,\theta)}}{\partial {\epsilon}^3}=-K(0, \theta) \end{equation} Till now we compute the derivative of the function $\sqrt{g_{22}(\epsilon,\theta)}$. We can use tham to approximate $\sqrt{g_{22}(\epsilon,\theta)}$ by using the Taylor series. By doing that we have: \begin{equation} \sqrt{g_{22}(\epsilon,\theta)}=\epsilon-K(0,\theta) \frac{\epsilon^3}{6}+..... \end{equation} So the Gaussian curvature is: \begin{equation} K(0,\theta)= \frac{6}{\epsilon^2} \left(1-\frac{\sqrt{g_{22}(\epsilon,\theta)}}{\epsilon} \right) \end{equation} The circumference is obtained by integrating $\sqrt{g_{22}(\epsilon,\theta)}$ over $\theta$ between $0$ and $2 \pi$. I suspect there is a better way to prove the equation based on the Gauss-Bonnet theorem. I don't have the time to dig it out, however since we are talking about geodesic circles, the geodesic curvature is by definition zero. The Gauss-Bonnet theorem then allows us to relate the Gaussian curvature to the area of the circle, which in turn is related to its circumference.