Reexpress $\frac{\partial^2 z}{\partial X_i \partial X_j}$, where $z(X, Y(X))$ subject to $\frac{\partial z}{\partial Y} = 0$ as a function of $X,Y$? Consider $C_{ij}$, the second derivative of $z(X, Y(X))$ with respect to $X$:
$$C_{ij} = \frac{\partial^2z}{\partial X_i \partial X_j}$$
where $Y$ satisfies the set of equations
$$\left. \frac{\partial z}{\partial Y_l}\right |_X = 0 \text{ for all } l$$
I'm trying to derive the following expression, as stated in https://dx.doi.org/10.1103/physrevb.1.3599 (equations 24 through 27):
$$C_{ij} = \frac{\partial^2z}{\partial X_i \partial X_j} = \left.\frac{\partial^2z}{\partial X_i \partial X_j}\right|_Y - \frac{\partial^2z}{\partial X_i \partial Y_k} R_{kl} \frac{\partial^2z}{\partial Y_l \partial X_j} $$
where
$$R_{kl} \frac{\partial^2z}{\partial Y_l \partial Y_m} = \delta_{km}$$
(I've changed the variables from $U$ to $z$, $\varepsilon$ to X and $Q$ to $Y$ to make it clear what is a scalar and what is a vector).
The motivation for this problem is that z is a potential that can vary with $X$ and $Y$. When $X$ is changed, $Y$ adjusts to minimize the value of $z$. We know how $z$ changes as a function of $X$ and $Y$.
I assume that I need to use the chain rule and Lagrange multipliers but have little experience with the latter. I would very much like to understand how this is done.

I'm guessing $\def\R{\mathbb{R}}z:\R^n\times\R^m\to\R$.
Define the functions $F_i:\R^n\times\R^m\to\R$ as
$$
\begin{align}
F_i(x,y) = \partial_{Y_i}z(x,y).
\end{align}
$$
Now this "$y$ adjusts to minimize $z(x,y)$" means that you have a function $\phi:\R^m\to\R^n$ such that
$$
F_i(x,\phi(x)) = 0.
$$
This $\phi$ is to be thought as the "smart $y$ that adjusts itself".
Differentiating this equation, we have, by the chain rule,
$$
\partial_{X_j}F_i(x,\phi(x))
+
\partial_{Y_k}F_i(x,\phi(x))
\partial_{X_j}\phi^k(x)
=
0
$$
where $\phi^k$ are the components of $\phi$ and I used the summation convention to sum over $k$.
Solving for $\partial_{X_j}\phi^k(x)$ we get
\begin{align}
\partial_{Y_k}F_i(x,\phi(x))
\partial_{X_j}\phi^k(x)
&=
-\partial_{X_j}F_i(x,\phi(x)) \\
\delta_{kl}
\partial_{X_j}\phi^k(x)
&=
-
R_{il}
\partial_{X_j}F_i(x,\phi(x)) \\
\partial_{X_j}\phi^l(x)
&=
-
R_{il}
\partial_{X_j}F_i(x,\phi(x)) \\
\partial_{X_j}\phi^k(x)
&=
-
R_{ik}
\partial_{X_j}F_i(x,\phi(x)) \tag 1
\end{align}
where $R$ is the inverse matrix of $\partial_YF=\partial_Y\partial_Yz$, as in your question.
Now put $f(x)=z(x,\phi(x))$ ("the $z$ minimized on its second parameter").
Again by the chain rule, we have
\begin{align}
\partial_{X_i}f(x) 
&=
\partial_{X_i}z(x,\phi(x))
+
\partial_{Y_k}z(x,\phi(x))
\partial_{X_i}\phi^k(x) \\
&=
\partial_{X_i}z(x,\phi(x))
\end{align}
where the second term vanished by the condition $F_i(x,\phi(x))=0$.
Taking another derivative, again by the chain rule, you have
\begin{align}
\partial_{X_j}\partial_{X_i}f(x) 
&=
\partial_{X_j}\partial_{X_i}z(x,\phi(x))
+
\partial_{Y_k}\partial_{X_i}z(x,\phi(x))
\partial_{X_j}\phi^k(x).
\end{align}
Now substitute $\partial_{X_j}\phi^k(x)$ with equation $(1)$ and you are done.
