The Laplacian of an expression involving the Ricci tensor While doing some computations on a compact Riemannian manifold I have reached the following expression:
$$ \Delta_y \big( Ric_y (\exp_y ^{-1} x, \exp_y ^{-1} x) \big) (x)$$
where $\Delta_y$ is the Laplacian with respect to $y$, $Ric_y$ is the Ricci tensor in $y$ and $\exp_y : U_y \subset TyM \to M$ is the exponential map at $y$. Note that, if $r = d(x,y)$, then each argument of $Ric_y$ is at least $O(r)$ for $r \to 0$, if not even better. Thus, the whole expression inside the Laplacian is at least $O(r^2)$.
Can you tell whether the expression above is $0$ or not? Explicitly doing the calculations does not seem a feasible choice. (In particular, if the whole expression inside the Laplacian were $O(r ^{2 + \varepsilon})$ for some $\varepsilon > 0$, the answer would be affirmative.)
 A: In general, I don't think so. If my interpretation is correct (see my comment above), then if you take an Einstein manifold with $Ric = \lambda g$ you immediately see that the best you can expect is $O(r^{2})$. In fact, if you take the standard sphere $\mathbb{S}^n$, the Taylor expansion is precisely $r^2 + \ldots$ which suggests that the Laplacian is exactly $n\neq 0$. 
It is late here so I don't feel like doing the computation, but thinking about it a bit more it seems to me that (meaning, I am making an educated guess) what you wrote down should be in fact equal to the scalar curvature at $x$. 

The vector $\exp_{y}^{-1}x$ is the vector $v$ at $T_yM$ such that the geodesic with initial position $y$ and initial velocity $v$ reaches $x$ after time 1. Since the geodesic flow is time-symmetric, you see that letting $w = \exp_{x}^{-1} y\in T_x M$, the vector $v$, after parallel transport to $x$ along the geodesic thus generated, is equal to $-w$. Now, in geodesic normal coordinates the radial lines are by definition unit-speed geodesics. So in the local coordinate chart given by the exponential map at $x$, we have that the vector $\exp_{y}^{-1} x$ is the vector $-y$ based at $y$ ($x$ is the origin).  
To be even more precise, this means that fixing geodesic normal coordinates at $x$ with $x$ the origin, and identifying $y$ with its image in $T_xM$, we have that the scalar function $Ric_y(\exp_{y}^{-1} x, \exp_{y}^{-1}x)$ can be alternatively expressed as
$$ Ric(y)_{ij} y^i y^j \tag{1}$$
where $Ric(y)_{ij}$ is the components of the Ricci tensor, expressed in this coordinate system, evaluated at the point $y$. The local expression of the Laplacian is
$$ \frac{1}{\sqrt{|g|}} \partial_k g^{kl} \sqrt{|g|} \partial_l = g^{kl}\partial^2_{kl} + O(\partial) $$ 
where the $O(\partial)$ denotes terms that are first order in derivative. Hitting this to the function in (1) you have exactly that the $O(\partial)$ terms vanish when $y = 0$ since (1) vanishes quadratically. So you have
$$ \triangle_y Ric(y)_{ij} y^i y^j |_{y = 0} = g(y)^{kl} \partial^2_{kl} Ric(y)_{ij} y^i y^j |_{y = 0} $$
Using again the quadratic vanishing you have that both of the derivatives must fall on the $y^i$ and $y^j$ terms, otherwise the expression vanishes if $g$ is smooth. There are two possibilities in terms of order of derivation, so you get
$$ \triangle_y Ric(y)_{ij} y^i y^j |_{y = 0} = \delta^{kl} Ric(0)_{ij} \left( \delta_k^i \delta_l^j + \delta_k^j \delta_l^i\right) = 2 Ric(0)_{ij} \delta^{ij} = 2 S $$
where $S$ is the scalar curvature at $x$. 
