Proof of gauge invariance of the massless Fierz-Pauli action (follow-up)

This question is a follow-up to Proof of gauge invariance of the massless Fierz-Pauli action.

One representation of the Fierz-Pauli action (up to a prefactor) is, $$S[h] =\int dx\left\{\underbrace{\frac{1}{2}(\partial_\lambda h^{\mu\nu})(\partial^\lambda h_{\mu\nu})}_{=:A}-\underbrace{\frac{1}{2}(\partial_\lambda h)(\partial^\lambda h)}_{=:B}-\underbrace{(\partial_\lambda h^{\lambda\nu})(\partial^\mu h_{\mu\nu})}_{=:C}+\underbrace{(\partial^\nu h)(\partial^\mu h_{\mu\nu})}_{=:D}\right\}.\tag{1}$$

We now want to show that $$S[h]$$ is invariant under the gauge transformation, $$h_{\mu\nu}\rightarrow h_{\mu\nu}+\delta h_{\mu\nu},\tag{2}$$ wherein $$\delta h_{\mu\nu}=\partial_\mu\xi_\nu+\partial_\nu\xi_\mu$$. We demand that $$\xi_\mu(x_\nu)$$ falls of rapidly at the respective boundaries of the action.

i) Why is it sufficient to only consider invariance of gauge transformations up to the first-order? Even if we consider the weak gravity regime $$h_{\mu\nu}\ll1$$, I don't see how this should lead to $$\delta h_{\mu\nu}\ll 1$$.

We now start to show first-order invariance by applying the gauge transformation, Eq. (2), to the terms $$A, B, C, D$$.

\begin{align} A &\to\frac{1}{2}(\partial_\lambda h^{\mu\nu}+\partial_\lambda \delta h^{\mu\nu})(\partial^\lambda h_{\mu\nu}+\partial^\lambda\delta h_{\mu\nu})\\ &=\underbrace{\frac{1}{2}(\partial_\lambda h^{\mu\nu})(\partial^\lambda h_{\mu\nu})}_{=A}+\underbrace{(\partial_\lambda h^{\mu\nu})(\partial^\lambda\delta h_{\mu\nu})}_{=\delta A}+\mathcal{O}(\delta h_{\mu\nu}^2)\\ B &\to\frac{1}{2}(\partial_\lambda h+\partial_\lambda \delta h)(\partial^\lambda h+\partial^\lambda\delta h)\\ &=\underbrace{\frac{1}{2}(\partial_\lambda h)(\partial^\lambda h)}_{=B}+\underbrace{(\partial_\lambda h)(\partial^\lambda\delta h)}_{=:\delta B}+\mathcal{O}(\delta h_{\mu\nu}^2)\\ C &\to(\partial_\lambda h^{\lambda\nu}+\partial_\lambda\delta h^{\lambda\nu})(\partial^\mu h_{\mu\nu}+\partial^\mu\delta h_{\mu\nu})\\ &=\underbrace{(\partial_\lambda h^{\lambda\nu})(\partial^\mu h_{\mu\nu})}_{=C}+\underbrace{2(\partial_\lambda h^{\lambda\nu})(\partial^\mu\delta h_{\mu\nu})}_{=:\delta C}+\mathcal{O}(\delta h_{\mu\nu}^2)\\ D &\to (\partial^\nu h+\partial^\nu\delta h)(\partial^\mu h_{\mu\nu}+\partial^\mu\delta h_{\mu\nu})\\ &=\underbrace{(\partial^\nu h)(\partial^\mu h_{\mu\nu})}_{=D}+2\underbrace{(\partial^\nu h)(\partial^\mu \delta h_{\mu\nu})}_{=:\delta D}+\mathcal{O}(\delta h_{\mu\nu}^2) \end{align} ii) Are these results correct so far? How do I show $$(\partial^\nu h)(\partial^\mu \delta h_{\mu\nu})=(\partial^\nu\delta h)(\partial^\mu h_{\mu\nu})$$?

Using the previous results, we find, $$S[h+\delta h]-S[h] =\int dx\left\{\delta A-\delta B-\delta C+\delta D\right\}+\mathcal{O}(\delta h^2).\tag{3}$$ Only $$\delta B$$ and $$\delta D$$ contain $$h$$, therefore, both should cancel (up to a constant) and we can consider them separate, \begin{align} \int dx\left\{\delta D-\delta B\right\} &=\int dx\left\{2(\partial^\nu h)(\partial^\mu\delta h_{\mu\nu})-(\partial_\lambda h)(\partial^\lambda\delta h) \right\}\\ &=\int dx(\partial^\lambda h)\left\{2(\partial^\mu\delta h_{\mu\lambda})-(\partial_\lambda\delta h) \right\}\\ &=\int dx(\partial^\lambda h)\left\{2(\partial^\mu(\partial_\mu\xi_\lambda+\partial_\lambda\xi_\mu)-\partial_\lambda(2\partial^\mu\xi_\mu) \right\}\\ &=2\int dx(\partial^\lambda h)(\partial^2\xi_\lambda).\tag{4} \end{align} Next, we examine the other two terms, \begin{align} \int dx\left\{\delta A-\delta C\right\} &=\int dx\left\{(\partial_\lambda h^{\mu\nu})(\partial^\lambda\delta h_{\mu\nu})-2(\partial_\lambda h^{\lambda\nu})(\partial^\mu \delta h_{\mu\nu})\right\}\\ &=\int dx\left\{-h^{\mu\nu}(\partial^2\delta h_{\mu\nu})+2h^{\lambda\nu}(\partial_\lambda\partial^\mu \delta h_{\mu\nu})\right\}\\ &=\int dxh^{\mu\nu}\left\{-\partial^2\delta h_{\mu\nu}+2\partial_\mu\partial^\lambda \delta h_{\lambda\nu}\right\}\\ &=\int dxh^{\mu\nu}\left\{-\partial^2(\partial_\mu\xi_\nu+\partial_\nu\xi_\mu)+2\partial_\mu\partial^\lambda (\partial_\lambda\xi_\nu+\partial_\nu\xi_\lambda)\right\}\\ &=\int dxh^{\mu\nu}\left\{\partial_\mu\partial^2\xi_\nu-\partial^2\partial_\nu\xi_\mu+2\partial_\mu\partial_\nu(\partial^\lambda\xi_\lambda)\right\},\tag{5} \end{align} wherein we used partial integration for the second equal and index relabeling for the third equal.

Comparing Eq. (4) and Eq. (5), we see that terms don't add up to a constant or divergence. iii) Where have I made mistakes?

i) Our gauge transformation is a linear transformation and therefore can be considered to form a Lie group. From Lie groups we know, that it is sufficient to show invariance only up to the first-order as we can always dissect transforms "large" in magnitude (think $$\delta h\gg1$$) into infinitesimal steps. If someone can put this in a more rigorous language, please do so!
ii)+iii) Actually, $$(\partial^\nu h)(\partial^\mu \delta h_{\mu\nu})\neq(\partial^\nu\delta h)(\partial^\mu h_{\mu\nu})$$, thus, we must correct the transformation of term $$D$$ to, $$D \to (\partial^\nu h+\partial^\nu\delta h)(\partial^\mu h_{\mu\nu}+\partial^\mu\delta h_{\mu\nu})\\ =\underbrace{(\partial^\nu h)(\partial^\mu h_{\mu\nu})}_{=D}+\underbrace{(\partial^\nu \delta h)(\partial^\mu h_{\mu\nu})+ (\partial^\nu h)(\partial^\mu \delta h_{\mu\nu})}_{=:\delta D}+\mathcal{O}(\delta h_{\mu\nu}^2).$$ Now, Eq. (4) reads, \begin{align} \int dx\left\{\delta D-\delta B\right\} &= \int dx(\partial^\nu h)\left\{\partial^\mu\delta h_{\mu\nu}-\partial_\nu \delta h\right\}+\int dx (\partial^\nu \delta h)(\partial^\mu h_{\mu\nu})\\ &=\underbrace{-\int dx h\partial^2\left\{\partial^\nu\xi_\nu-\partial^\mu\xi_\mu\right\}}_{=0}-\int dx h_{\mu\nu}(\partial^\mu\partial^\nu\delta h). \end{align} Adding up Eq. (5) and the corrected version of Eq. (4), we find that the transformed action up to first-order indeed vanishes, \begin{align} \int dx\delta S &=\int dx h^{\mu\nu}\left\{\partial_\mu\partial^2\xi_\nu-\partial^2\partial_\nu\xi_\mu+\underbrace{2\partial_\mu\partial_\nu(\partial^\lambda\xi_\lambda)-2\partial_\mu\partial_\nu(\partial^\lambda\xi_\lambda)}_{=0}\right\}\\ &=\int dx h^{\mu\nu}\partial^2\partial_\mu\partial^2\xi_\nu-\int dx h^{\nu\mu}\partial^2\partial_\mu\xi_\nu\\ &=\int dx h^{\mu\nu}\partial^2\partial_\mu\partial^2\xi_\nu-\int dx h^{\mu\nu}\partial^2\partial_\mu\xi_\nu =0, \end{align} where we have used in the last steps that $$h^{\mu\nu}=h^{\nu\mu}$$ and that we can relabel summed indices.