Solution (source code)

= Solution

Reality of $E_{ij}^T$ and the given polarization relation imply $\psi_{(p)}(-\mathbf k)=\psi_{(p)}(\mathbf k)^*$. Insert the Fourier expansion twice into each quadratic term. The spatial integral gives
$$
\int d^3x\,e^{i(\mathbf k+\mathbf k')\cdot\mathbf x}
=(2\pi)^3\delta^3(\mathbf k+\mathbf k').
$$
This cancels the two factors of $(2\pi)^{-3/2}$. Set $\mathbf k'=-\mathbf k$ and use $M_{ij}^{(q)}(-\mathbf k)=M_{ij}^{(q)}(\mathbf k)^*$. The unit-normalization condition eliminates cross-polarization terms, giving
$$
\int d^3x\,E_{ij}^{T\prime}E^{Tij\prime}
=\sum_p\int d^3k\,|\psi_{(p)}'(\mathbf k)|^2.
$$
For the gradient term, $(ik_l)(ik_l')$ becomes $k^2$ after the delta-function substitution. Thus
$$
\boxed{S^{(2)}=\frac{M_{\rm Pl}^2}{2}\sum_{p=\pm2}
\int d\eta\,d^3k\,a^2\left[
|\psi_{(p)}'|^2-k^2|\psi_{(p)}|^2
\right]}.
$$
The full momentum-space integral and its reality relation are those of real fields; there is no additional independent copy of the negative-wavevector modes.