Is there a more efficient way to compute the second functional derivative of Ricci tensor, \begin{equation} \frac{\delta^2 R_{\mu \nu}(x)}{\delta g^{\alpha \beta}(y) \delta g^{\gamma \epsilon}(z)} \end{equation} rather than resorting to a brute-force calculation of Christoffel symbol derivatives?
I know that the first variation of Ricci tensor is the following: \begin{equation} \delta R_{\mu \nu}=\frac{1}{2} g^{\eta \xi}\left[\nabla_\eta \nabla_\nu \delta g_{\mu \xi}+\nabla_\eta \nabla_\mu \delta g_{\nu \xi}-\nabla_\mu \nabla_\nu \delta g_{\eta \xi}-\nabla_\eta \nabla_{\xi} \delta g_{\mu \nu}\right] \end{equation} Using the relation $\frac{\partial g_{\mu \nu}}{\partial g^{\sigma \rho}}=-\frac{1}{2}\left(g_{\mu \rho} g_{\sigma v}+g_{\mu \sigma} g_{\rho \nu}\right)$ one obtains: \begin{equation} \begin{aligned} \frac{\delta R_{\mu \nu}(x)}{\delta g^{\alpha \beta}(y)} &= - \frac{1}{4} \left[ g_{\mu \alpha} \nabla_\beta \nabla_\nu + g_{\mu \beta} \nabla_\alpha \nabla_\nu + g_{\nu \alpha} \nabla_\beta \nabla_\mu + g_{\nu \beta} \nabla_\alpha \nabla_\mu - 2 g_{\alpha \beta} \nabla_\mu \nabla_\nu \\ - (g_{\mu \alpha} g_{\nu \beta} + g_{\mu \beta} g_{\nu \alpha}) \nabla^2 \right] \delta(x-y) \end{aligned} \end{equation} However, taking the second derivative seems cumbersome to me. Are there alternative approaches for computing higher-order functional derivatives of curvature tensors?
