Behavior of sectional curvature under metric deformations Metric deformation:
Let $(M,g_0)$ be a Riemannian manifold and consider a (sufficiently smooth) deformation of $g_0$, $$g_t=g_0+th+O(t^2), \quad 0< t<\varepsilon $$ where $h$ is some symmetric (0,2)-tensor. A natural (and important) question is how the sectional curvatures of $g_0$ change under this deformation, e.g., what is the infinitesimal change in terms of $h$. More precisely, given two $g_0$-orthonormal vectors $X$ and $Y$ in $T_pM$, define the (unnormalized) sectional curvature $$k(t)=g_t(R_t(X,Y)Y,X),$$ where we are using the appropriate sign convention on $R$.

Q: What is the explicit formula for $k'(0)=\frac{d}{dt}k(t)\big|_{t=0}$?


Possible (but different?) answers:
I have found a few papers with an answer, but (understandably) none provide the complete argument. Unfortunately, it seems like some of them are really different, and it would be very helpful if someone could point out if they coincide for some (possibly silly) reason I am not seeing.

*

*Berger'66 (Trois remarques sur les variétés riemanniennes à courbure positive)/Bourguignon, Deschamps, Sentenac'72 (Conjecture de H. Hopf sur les produits de variétés):
$$k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)$$


*Strake'87 (Curvature increasing metric variations):
$$k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+h(R_0(X,Y)Y,X)$$
$$-k(0)(h(X,X)+h(Y,Y))$$
3. Topping'06 (Lectures on Ricci Flow): $$k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+\tfrac12h(R_0(X,Y)Y,X)$$
$$-\tfrac{1}{2}h(R_0(X,Y)X,Y)$$
EDIT: I should point out that, although formula 3) as I wrote above is NOT a correct expression for $k'(0)$, I (embarrassingly) misinterpreted it from Prop. 2.3.5 in Topping's lecture notes -- which actually contains the correct formula (matching Vitali's answer below). He gives a general expression for $\frac{d}{dt}g_t(R^t(X,Y)Z,W)\big|_{t=0}$, and by using the Ricci identity it becomes clear that his formula is indeed the same as Vitali's. I sincerely apologize for the confusion.
Note that all answers above coincide if $(M,g_0)$ has non-negative sectional curvature and $X$ and $Y$ span a plane of zero $g_0$-curvature. As Strake remarks, the self-adjoint endomorphism $A_X Z=R(Z,X)X$ is positive-semidefinite and $g(A_X Y,Y)=0$. (Nevertheless, to the best of my understanding, THESE HYPOTHESES ARE NOT ASSUMED in the references in 1). Also, it seems to me that answers 2 and 3 DO NOT COINCIDE.
 A: I guess I say this a lot, but this is something you really should work out yourself. You might get lost in the calculations the first few or many times you do it, but after a while you should get the hang of it.
It's probably best to do it first in local co-ordinates. It's relatively straightforward to compute the variation of the Christoffel symbol and the Riemann curvature tensor written with respect to the co-ordinates. It is therefore relatively easy to compute the first variation of $R(X,Y)Y\cdot X$, where $X$ and $Y$ are assumed to be fixed, independent of the variation.
Finally, it is important to note that even if you assume $X$ and $Y$ are orthonormal for the original metric, they do not remain so under the variation. So you have to normalize $R_t(X,Y)Y\cdot X$ to get the sectional curvature. So you have to compute the variation of the normalization, too.
If you do this enough times, you will know which of the formulas above are right. That's a lot better than taking anyone else's word for it. And you will develop a much greater facility for doing such calculations.
A: Formula 2) is the correct one in general except it's the derivative of the sectional curvature i.e of $\frac{k_t(X,Y)}{|X\wedge Y|^2_t}$  (and not just of $k_t(X,Y)$) for an orthonormal frame $X,Y$ with respect to the original metric. This accounts for the last term in formula 2.
For $k'(0)$ itself the correct formula is
$$
k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+h(R(X,Y)Y,X)
$$
As Deane said, this is a straightforward calculation. But since the OP seems to be struggling with it I'll supply some details. Let $\nabla^t$ be the Levi-Civita connection for $g_t$. Then $\nabla^t_XY=\nabla_XY+tS(X,Y)+O(t^2)$ where $S$ is a (2,1)-tensor and $\nabla=\nabla^0$.
Let's first compute $k'(0)$ in terms of $S$. As usual let's work in normal coordinates around $p$ with $X,Y$ coordinate fields.
We have $$\langle R^t(X,Y)Y,X\rangle_t=\langle R^t(X,Y)Y,X\rangle_0+t\cdot h(R(X,Y)Y,X)+O(t^2)=$$
$$=\langle \nabla^t_X\nabla^t_YY,X\rangle_0-\langle \nabla^t_Y\nabla^t_XY,X\rangle_0+t\cdot h(R(X,Y)Y,X)+O(t^2)$$
 Next we expand the first term
$$\langle \nabla^t_X\nabla^t_YY,X\rangle_0=\langle \nabla^t_X(\nabla_YY+tS(Y,Y)),X\rangle_0+O(t^2)=$$
$$\langle \nabla_X\nabla_YY,X\rangle_0+t\langle S(X,\nabla_YY)+\nabla_XS(Y,Y), X\rangle_0+O(t^2)=$$
$$  \langle \nabla_X\nabla_YY,X\rangle_0+t\langle \nabla_XS(Y,Y), X\rangle_0+O(t^2)$$
where in the last equality we used that $\nabla_YY(p)=0$.
After a similar computation for $\langle \nabla^t_Y\nabla^t_XY,X\rangle_0$ we get that
$$k'(0)=\langle\nabla_XS(Y,Y)-\nabla_YS(X,Y), X\rangle_0+h(R(X,Y)Y,X)$$
$$=\nabla_X S(Y,Y,X)-\nabla_YS(X,Y,X)+h(R(X,Y)Y,X)$$
where we lowered the index and turned $S$ into a $(3,0)$-tensor $S(X,Y,Z)=\langle S(X,Y),Z\rangle_0$
Lastly, recall that
$\langle \nabla_XY,Z\rangle=\frac{1}{2}[X\langle Y,Z\rangle+Y\langle X,Z\rangle -Z\langle X, Y\rangle]$ for coordinate fileds. This easily gives
$S(Y,Y,X)=\frac{1}{2}[Yh(X,Y)+Yh(X,Y)-Xh(Y,Y)]=\nabla_Yh(X,Y)-\frac{1}{2}\nabla_Xh(Y,Y)$ and 
$S(X,Y,X)=\frac{1}{2}\nabla_Yh(X,X)$ written invariantly as tensors.
Altogether this gives
$$
k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+h(R(X,Y)Y,X)
$$ as promised.
Formulas 1) and 3) are not true in general but might be true in the specific circumstances where they are applied in the papers in question. 
