Deriving the Fierz-Pauli propagator Fierz-Pauli action can be written as
$$ 
\tag{1} S = \int d^Dx \; \frac{1}{2} h_{\mu \nu} \mathcal{O}^{\mu \nu, \alpha \beta}h_{\alpha \beta}
$$
where the operator
$$
\tag{2} \mathcal{O}^{\mu \nu} {}_{\alpha \beta} = \left( \eta^{( \mu}{}_{\alpha}\eta^{\nu )}{}_{\beta} - \eta^{\mu \nu} \eta_{\alpha \beta}\right) \left(\Box - m^2 \right) - 2\partial^{(\mu}\partial_{(\alpha}\eta^{\nu )}{}_{\beta )} + \partial^{\mu}\partial^{\nu}\eta_{\alpha \beta} + \partial_{\alpha}\partial_{\beta}\eta^{\mu \nu}
$$
satisfies the symmetries 
$$
\tag{3} \mathcal{O}^{\mu \nu , \alpha \beta} = \mathcal{O}^{\nu \mu , \alpha \beta} = \mathcal{O}^{\mu \nu , \beta \alpha} = \mathcal{O}^{\alpha \beta , \mu \nu}.
$$
Now, to find the propagator one rewrites this operator in momentum space as
$$
\tag{4} \mathcal{O}^{\mu \nu}{}_{\alpha \beta}(\partial \to ip) = -\left( \eta^{( \mu}{}_{\alpha}\eta^{\nu )}{}_{\beta} - \eta^{\mu \nu} \eta_{\alpha \beta}\right) \left(p^2 + m^2 \right) \\ + 2p^{(\mu}p_{(\alpha}\eta^{\nu )}{}_{\beta )} - p^{\mu}p^{\nu}\eta_{\alpha \beta} - p_{\alpha}p_{\beta}\eta^{\mu \nu}.
$$
I know that propagator $ \mathcal{D}_{\alpha \beta, \sigma \lambda}$ satisfies the same symmetries above so it will have the same structure as $\mathcal{O}^{\mu \nu, \alpha \beta}$. I also know that their multiplication should give me the identity operator, namely
$$
\tag{5} \mathcal{O}^{\mu \nu, \alpha \beta}\mathcal{D}_{\alpha \beta, \sigma \lambda} = \frac{i}{2}\left(\delta^{\mu}_{\sigma} \delta^{\nu}_{\lambda} + \delta^{\nu}_{\sigma}\delta^{\mu}_{\lambda}\right).
$$
To find the propagator, I should assume a form for it, put that form into the above equation and find the coefficients. My question is, what should be the form of the propagator in that case?
I suppose, after determining the form of the propagator, calculations will be lengthy. So, before trying something wrong and wasting my time, I want to learn the correct form (the correct assumption) for this propagator. I'll be glad if anybody can help.
 A: The propagator has to be built out of $\eta^{\mu\nu}$ and $p^\mu$, in hind sight it will prove more convenient to build it equivalently out of $D^{\mu\nu}\equiv \eta^{\mu\nu} - p^\mu p^\nu / m^2$ and $p^\mu$ (these two are easily understood to be equivalent expansions).
Now there are some symmetry considerations; exchange of $\mu\nu$ with $\rho\sigma$, symmetric under exchange of $\mu$ and $\nu$ separately as well as $\rho$ and $\sigma$.
This fixes 
$$\mathcal{D}_{\mu\nu,\rho\sigma} = \alpha D_{\mu\nu}D_{\rho\sigma} + \beta D_{\mu(\rho}D_{\sigma)\nu} + \gamma(D_{\mu\nu}p_{\rho}p_{\sigma} + D_{\rho\sigma}p_{\mu}p_{\nu}) + \delta ( D_{\mu(\rho}p_{\sigma)}p_{\nu} + D_{\nu(\rho}p_{\sigma)}p_{\mu}) + \epsilon p_\mu p_\nu p_\rho p_\sigma
$$
What's left is to find the coefficients.
A: Let's start with a basic example. For the case $m=0$ the operator (1) can be rewritten as
$$
\mathcal{O}^{\mu \nu, \alpha \beta} = \Box \left[ \frac{1}{2}(\eta^{\mu \alpha} \eta^{\nu \beta} + \eta^{\mu \beta} \eta^{\nu \alpha}) - \frac{1}{2}\eta^{\mu \nu}\eta^{\alpha \beta}  \right],
$$
after gauge freedom is fixed. In momentum space we have
\begin{align}
\mathcal{O}^{\mu \nu, \alpha \beta} & = -p^2 \left[ \frac{1}{2}(\eta^{\mu \alpha} \eta^{\nu \beta} + \eta^{\mu \beta} \eta^{\nu \alpha}) - \frac{1}{2}\eta^{\mu \nu}\eta^{\alpha \beta}  \right]
\\
& = -p^2 \left[ \eta^{\mu ( \alpha}\eta^{\beta) \nu} - \frac{1}{2}\eta^{\mu \nu}\eta^{\alpha \beta}  \right].
\end{align}
This operator also satisfies the symmetries (3) and again propagator should satisfy the equation (5). We should contain all the possible Lorentz invariant combinations of $\eta^{\mu \nu}$ (with suitably placed indices) when assuming the form of the propagator. We also note that propagator satisfies the same symmetries (3) with the operator $\mathcal{O}^{\mu \nu, \alpha \beta}$ and our metric tensor $\eta^{\mu \nu}$ is symmetric on its indices. Hence in our assumption for the propagator $\mathcal{D}_{\alpha \beta, \sigma \lambda}$ if we use (and we will indeed) the term $\eta_{\alpha \beta}\eta_{\sigma \lambda}$ for example, then we should ignore the terms $\eta_{\alpha \beta}\eta_{\lambda \sigma}, \; \eta_{\beta \alpha}\eta_{\sigma \lambda}, \;  \eta_{\beta \alpha}\eta_{\lambda \sigma}, \; \eta_{\sigma \lambda}\eta_{\alpha \beta}, \; \eta_{\sigma \lambda}\eta_{\beta \alpha}, \; \eta_{\lambda \sigma}\eta_{\alpha \beta}, \; \eta_{\lambda \sigma}\eta_{\beta \alpha}$ since they are all equivalent to each other and to $\eta_{\alpha \beta}\eta_{\sigma \lambda}$ by given symmetry considerations. Therefore, the correct form for the propagator should be 
$$
\mathcal{D}_{\alpha \beta , \sigma \lambda} = A(p)\eta_{\alpha ( \sigma}\eta_{\lambda) \beta} + B(p)\eta_{\alpha \beta} \eta_{\sigma \lambda},
$$
which by using (5) gives
$$
A(p) = -\frac{i}{p^2}, \quad B(p) = \frac{i}{p^2(D-2)}.
$$

Now for our question at hand, from equation (4) we understand that we should write down all Lorentz invariant combinations of $\eta^{\mu \nu}$ and $p^{\mu}$. However it is very lengthy and complicated to write them all and work with them. Hence we first write equation (4) in its longest form by explicitly showing all the indices, i.e; we get rid of all the parentheses we used to declare the symmetries on indices. Then by inspection we realize that with proper constants we can rewrite (4) by using
$$
\tag{6} P^{\mu \nu} \equiv \eta^{\mu \nu} + \frac{1}{m^2}p^{\mu}p^{\nu},
$$
and its combinations. Therefore in assuming the form of the propagator we can also use (6) instead of all the $\eta_{\mu \nu}$'s and $p_{\mu}$'s. Again by regarding the symmetry considerations, using $P_{\alpha \beta}P_{\sigma \lambda}$ for example, makes the terms $P_{\alpha \beta}P_{\lambda \sigma}, \; P_{\beta \alpha}P_{\sigma \lambda}, \;  P_{\beta \alpha}P_{\lambda \sigma}, \; P_{\sigma \lambda}P_{\alpha \beta}, \; P_{\sigma \lambda}P_{\beta \alpha}, \; P_{\lambda \sigma}P_{\alpha \beta}, \; P_{\lambda \sigma}P_{\beta \alpha}$ useless. Hence we assume the form
$$
\mathcal{D}_{\alpha \beta, \sigma \lambda} = A(p)P_{\alpha (\sigma}P_{\lambda) \beta} + B(p)P_{\alpha \beta}P_{\sigma \lambda},
$$
for the propagator. Using (5) gives 
$$
A(p) = -\frac{i}{p^2 + m^2}, \quad B(p) = \frac{i}{p^2 + m^2} \frac{1}{D-1},
$$
which is the correct result for the massive gravity propagator.

NOTE: In above calculations $D$ represents the dimension of our space-time and appears as a result of the multiplication $\eta^{\mu \nu}\eta_{\mu \nu}$.
