The massless Fierz-Pauli action describing a spin-2 field $h_{\mu\nu}$ is (up to a prefactor) given by, $$ S[h]=\int dx h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu} h_{\mu\nu},\tag{1} $$ wherein we define the differential operator, $$ \zeta_{\alpha\beta}^{\mu\nu}=\square\left(P_\alpha^\mu P_\beta^\nu-P_{\alpha\beta}P^{\mu\nu}\right),\tag{2} $$ with projection tensor $P_{\mu\nu}=\eta_{\mu\nu}-\partial^{-2}\partial_\mu\partial_\nu$ and d'Alembertian $\square$.
Many references, e.g. Hinterbichler. 2011, claim that such action, Eq. (2), is invariant under the gauge transformation, $$ h_{\mu\nu}\to h_{\mu\nu}+\delta h_{\mu\nu}=h_{\mu\nu}+\partial_\mu\xi_\nu+\partial_\nu\xi_\mu, $$ wherein we demand $\xi_\mu(x_\mu)$ to be continous differentiable and to fall of sufficient fast at infinity such that boundary terms vanish.
How do I prove the claimed gauge invariance?
We claim a theory invariant under a specific transformation if the equations of motion (EOMs) remain unchanged. From classical mechanics, we know that the EOMs remain unchanged if the action is changed by a total time derivative or a constant term as these drop out of the Euler-Lagrange equations which lead to the EOMs. I believe the time derivative is not relevant if we consider spacetime as we cannot easily separate time from space, thus in our case we are left to show, $$ S[h+\delta h]-S[h]=\text{const}.\tag{3} $$
When inserting Eq. (1) into Eq. (3) I struggle with the final steps. Furthermore, I would be grateful for tricks on how to simplify my calculations.
Calculations
We insert Eq. (1) into Eq. (3) and find that the term without $\delta h$ cancels out, $$ \begin{align} S[h+\delta h]-S[h] &=\int dx (h^{\alpha\beta}+\delta h^{\alpha\beta})\zeta_{\alpha\beta}^{\mu\nu}(h_{\mu\nu}+\delta h_{\mu\nu})-\int dx h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu} h_{\mu\nu}\\ &=\int dx \left\{h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu}\delta h_{\mu\nu}+\delta h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu} h_{\mu\nu}+\delta h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu}\delta h_{\mu\nu}\right\}.\tag{A.1} \end{align} $$ We note that the first two terms need to cancel each other as the these are the only terms that contain $h_{\mu\nu}$. Consequently, the third term has to be a constant.
We perform partial integration on the second term in Eq. (A.1), $$ \int dx\delta h^{\alpha\beta} \zeta_{\alpha\beta}^{\mu\nu} h_{\mu\nu} =-\int dx h_{\mu\nu}\left(\zeta_{\alpha\beta}^{\mu\nu}\delta h^{\alpha\beta}\right),\tag{A.2} $$ where we used that $\xi_\mu$ falls of rapidly towards the boundaries. That said, I am not sure if it is justified to use partial integration with $\zeta$ as the differential.
Using the Minkowski metric, we can raise and lower indices, $$ h_{\mu\nu}\zeta_{\alpha\beta}^{\mu\nu}\delta h^{\alpha\beta} =h^{\sigma\rho}\left(\eta_{\mu\sigma}\eta_{\nu\rho}\zeta_{\alpha\beta}^{\mu\nu}\eta^{\alpha\lambda}\eta^{\beta\gamma}\right)\delta h^{\alpha\beta} =h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu}\delta h_{\mu\nu}.\tag{A.3} $$ In the last step we relabeled the indices such that they match the first term in Eq. (A.1).
We are left with the third term, $$ \int dx\delta h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu}\delta h_{\mu\nu} =\int dx (\partial^\alpha\xi^\beta)\zeta_{\alpha\beta}^{\mu\nu}(\partial_\mu\xi_\nu+\partial_\nu\xi_\mu)+\int dx (\partial^\beta\xi^\alpha)\zeta_{\alpha\beta}^{\mu\nu}(\partial_\mu\xi_\nu+\partial_\nu\xi_\mu). \tag{A.4} $$ Because of the tensor symmetry $\zeta_{\alpha\beta}^{\nu\mu}=\zeta_{\alpha\beta}^{\mu\nu}=\zeta_{\beta\alpha}^{\mu\nu}$, we can sum the terms in Eq. (A.4) to, $$ \int dx\delta h^{\alpha\beta}\zeta_{\alpha\beta}^{\mu\nu}\delta h_{\mu\nu} =4\int dx (\partial^\alpha\xi^\beta)\zeta_{\alpha\beta}^{\mu\nu}(\partial_\mu\xi_\nu).\tag{A.5} $$ At this point, I don't see any obvious operations on how to show that (A.5) is constant.
