Solution (source code)

= Solution

Let the symmetric metric perturbation $h^{\alpha\beta}$ be varied with compact support. Integrating the mixed derivative product by parts makes its integral equal to that of $A_\sigma A^\sigma$. Thus the supplied <massless Fierz-Pauli action> has, up to a boundary term, density
$$
\mathcal L=-\frac14\partial_\mu h^{\rho\sigma}\partial^\mu h_{\rho\sigma}
+\frac12A_\sigma A^\sigma
-\frac12A^\nu\partial_\nu h
+\frac14\partial_\mu h\partial^\mu h.
$$
For example, the equivalence of the mixed term follows from commuting flat-space partial derivatives after moving one derivative off $h^{\rho\sigma}$. Here $\delta h=\eta_{\alpha\beta}\delta h^{\alpha\beta}$ and $\delta A^\nu=\partial_\mu\delta h^{\mu\nu}$.

Integrating each variation by parts, the four terms contribute respectively
$$
\begin{aligned}
\delta S=\int d^4x\,\delta h^{\alpha\beta}\bigg[
&\frac12\Box h_{\alpha\beta}
-\frac12(\partial_\alpha A_\beta+\partial_\beta A_\alpha)\\
&+\frac12\partial_\alpha\partial_\beta h
+\frac12\eta_{\alpha\beta}\partial_\mu A^\mu
-\frac12\eta_{\alpha\beta}\Box h\bigg].
\end{aligned}
$$
The symmetrization in the second term is required because the varied field is symmetric. Comparing this coefficient with the <Einstein tensor> perturbation above gives
$$
\boxed{\delta S=-\int d^4x\,\delta h^{\alpha\beta}\delta G_{\alpha\beta}.}
$$
Thus arbitrary compactly supported variations give \b[$\delta G_{\alpha\beta}=0$], exactly the vacuum <Linearized Einstein equations>. The minus sign and overall normalization of the action do not change those equations; boundary conditions justify the discarded total derivatives.