Product rule of variations I am deriving the Einstein equation using the Einstein-Hilbert action:  
It is obvious that the variation in the Riemann Tensor is calculated from a variational product rule. What is not obvious to me is why variations obey this rule, and I'll like an explanation.
 A: The easiest way to go about this is to think of variation in terms of the background field method. Here, you write the metric as some background plus perturbation:
\begin{equation}
g_{\mu \nu} = \bar{g}_{\mu \nu} + h_{\mu \nu}
\end{equation}
and then you expand every quantity that depends on $g_{\mu \nu}$ as a power series in $h$:
$$
\Gamma = \bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots \\
R = \bar{R} + R^{(1)} + R^{(2)} \cdots
$$
where the superscript $^{(n)}$ labels the $\mathcal{O} (h^n)$ perturbation. You can then easily see that
$$
\Gamma \Gamma = (\bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots)(\bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots)
$$
Extracting the $\mathcal{O}(h)$ variation, like in your problem, we automatically get both the contributions:
$$
(\Gamma \Gamma)^{(1)} = \bar{\Gamma} \Gamma^{(1)} +  \Gamma^{(1)} \bar{\Gamma}
$$
A: It's a fundamental property of infinitesimal variations $\delta$ that they are linear derivations, i.e. they obey Leibniz product rule. Short of a rigorous definition of infinitesimals, the proof is basically to consider variations as 1-parameter families of functions, and then differentiate wrt. the parameter, cf. above comment by Bence Racsko.
