For some purpose, I have to calculate the well-known linearization of the Einstein-Hilbert action in the de Sitter background. I encountered a problem: assuming the de Sitter metric, my resulting action has a Fierz-Pauli mass term with different coefficients. Indeed, let $g\leadsto g+h$ be a change in metric, where $g$ is the background de Sitter metric, and $h$ is the perturbation. Then, using this StackExchange answer: \begin{align} &\int d^4x \sqrt{-g}R \nonumber \\ \leadsto& \int d^4x \sqrt{-g}\left(1+\frac{1}{2}h+\frac{1}{8}(h^2-2h^\rho_\lambda h^\lambda_\rho)+\cdots \right)\left[ \vphantom{\frac{0}{0}}R +(-h^{\alpha \beta}R_{\alpha \beta}+\nabla_\beta \nabla_\alpha h^{\alpha \beta}-\square h)\right. \nonumber \\ &+\left( 2h^\gamma_\alpha h^{\alpha \beta}R_{\beta \gamma}-h^{\alpha \beta}h^{\gamma\lambda}R_{\alpha \gamma \beta \lambda}+h^{\alpha \beta}\nabla_\beta \nabla_\alpha h +\frac{1}{4} h \square h + h^{\alpha \beta} \nabla_\alpha \nabla_\gamma h^\gamma_\beta \right. \nonumber \\ &\left.\left.- h \nabla^\beta \nabla_\gamma h^\gamma_\beta-2h^{\alpha \beta}\nabla_\gamma \nabla_\beta h^\gamma_\alpha+\frac{1}{2} h_{\alpha \gamma}\nabla_\beta \nabla^\gamma h^{\alpha \beta}+\frac{1}{4} h_{\alpha \beta} \square h^{\alpha \beta}\right)+\cdots \right]\tag{1} \end{align} But: \begin{align} -2h^{\alpha \beta} \nabla_\gamma \nabla_\beta h^\gamma_\alpha &= -2h^{\alpha \beta} \nabla_\beta \nabla_\gamma h^\gamma_\alpha-2h^{\alpha \beta}R_{\rho \beta}h^\rho_\alpha+2h^{\alpha \beta}R^\rho_{\ \alpha \gamma \beta}h^{\gamma}_\rho\tag{2} \\ \frac{1}{2} h_{\alpha \gamma}\nabla_\beta \nabla^\gamma h^{\alpha \beta} &= \frac{1}{2} h_{\alpha \gamma} \nabla^\gamma \nabla_\beta h^{\alpha \beta}+\frac{1}{2}h^\gamma_\alpha R^\alpha_{\ \rho \beta \gamma} h^{\rho \beta}+\frac{1}{2}h^\gamma_\alpha R_{\rho \gamma} h^{\alpha \rho}\tag{3} \end{align} Thus we have, after some gathering of terms: \begin{align} &\int d^4x \sqrt{-g}R \nonumber \\ \leadsto&S_\text{EH}+S_\text{FP}+\int d^4x \sqrt{-g} \left( -h_{\mu \nu} G_{\mu \nu}+\frac{1}{8}(h^2-2h^\rho_\lambda h^\lambda_\rho)R -\frac{1}{2}(hh^{\mu \nu}-h^\mu_\alpha h^{\alpha \nu})R_{\mu \nu} \right. \nonumber \\ &\left.+\frac{3}{2}h^{\alpha \beta}h^{\gamma\lambda}R_{\alpha \gamma \beta \lambda} \right)+\cdots\tag{4} \end{align} My problem arises from $-\frac{1}{2}(hh^{\mu \nu}-h^\mu_\alpha h^{\alpha \nu})R_{\mu \nu}$. Indeed in de Sitter background, $R_{\mu \nu}=\frac{1}{4}g_{\mu \nu}R$, which means that $\frac{1}{8}(h^2-2h^\rho_\lambda h^\lambda_\rho)R$ and $\frac{1}{2}(hh^{\mu \nu}-h^\mu_\alpha h^{\alpha \nu})R_{\mu \nu}$ will not compensate (one is zero, but the other isn't). But if they don't, this means that I will have a Fierz-Pauli theory with two different masses in absolute value, which is known to be inconsistent [1].
I can't for my life find where I am wrong with my derivation...
So my question: are there any resources where this calculation has been done step by step, so I can compare and see where is my error? If someone sees where is my mistake, it would be great too.
[1] https://arxiv.org/abs/2102.10813 eq 2.3-2.5
