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.