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.