I am currently in the process of reading Riemannian Geometry by Gallot, Hulin and Lafontaine in order to learn more about differential and, well, Riemannian geometry. I'm more of a functional analyst so I'm trying to get accustomed to all the algebraic and geometric concepts and arguments, but not without issues...
Theorem $2.51$ is the theorem about the existence and uniqueness of the Levi-Civita connection, denoted by $D$, on a (pseudo-)Riemannian manifold $(M,g)$ (only the non-degenerescence of $g$ is used, hence $g$ can have any signature). In this book it is proved using a coordinate-free approach via the Koszul formula: $$2g(D_XY, Z) = L_X(g(Y,Z)) + L_Y(g(Z,X)) - L_Z(g(X,Y)) + g([X,Y], Z) - g([Y, Z], X) - g([X,Z], Y)$$
I have no problem with the uniqueness part, but I have trouble proving that this expression will serve to define a connection, i.e. a $\mathbb{R}$-bilinear map $D: \Gamma(TM)^2 \to \Gamma(TM)$ which is $\mathcal{C}^\infty(M)$-linear in the first variable and satisfies the Leibniz formula in respect to the second variable.
The advice from the authors given in that part of the proof suggests to check that the RHS of the formula, let's call it $P_{X,Y}(Z)$, is $\mathcal{C}^\infty(M)$-linear in regards to $Z$. However, what I have using the properties of the bracket and the metric and what I learned about Lie derivatives is this: $$P_{X,Y}(fZ) = f P_{X,Y}(Z) + fL_Z(g(X,Y)) - L_{fZ}(g(X,Y))$$ This should be good if the Lie derivative $L_X$ is $\mathcal{C}^\infty(M)$-linear in $X$ but at this step I wasn't too sure, because there is no formula given for $L_{fZ}$ relative to $L_Z$ and $f$ as far as I could see, be it in this book or on the net, including MSE.
Ideally I would have liked to prove it coordinate-free but here's my simple-and-to-the-point attempt using local charts:
Let $(x_1, \cdots, x_n)$ be local coordinates on a chart $U$ on $M$ and $h \in \mathcal{C}^\infty(M)$. In these coordinates, we have: $$L_{fX}(h) = \sum_{i = 1}^n (fX)_i \frac{\partial h}{\partial x_i} = \sum_{i = 1}^n f X_i \frac{\partial h}{\partial x_i} = f L_X(h)$$
There was a formula on Wikipedia and on MSE for the action of the Lie derivative on differential forms: $$L_{fX}(\omega) = fL_X(\omega) + df \wedge \iota_X\omega$$ but I'm not sure whether this works on functions considering that the interior product $\iota_X$ is the contraction $c_{1,1}(X \otimes \cdot)$ which eats contravariant tensors of type $(0, q)$ for $q \geq 1$, and thus not functions. What does worry me a bit is that this does not seem to be $\mathcal{C}^\infty(M)$-linear in regard to $X$?
Hence my question is:
Is the action of the Lie derivative on functions $\mathcal{C}^\infty(M)$-linear in $X$? If not, what is wrong with what I've written?