Let $\phi$ be a scalar field and $g_{\mu \nu} = \eta_{\mu \nu}+h_{\mu \nu}/M_p$ where $M_p$ is the Planck mass (so we assume we deal with perturbations). Let $\Lambda_2,\Lambda_3$ be energy scales such that $\Lambda_2 \gg \Lambda_3$. These are defined by $\Lambda^2_2 = M_p H_0$ and $\Lambda_3^3 = M_p H_0^2$. The Horndeski action is:

$$S = \int d^4 x \sqrt{-g} \sum^5_{i=2} \mathcal{L}_i,$$

where

\begin{align} \mathcal{L_2}&=\Lambda_2^4 G_2,\nonumber \\ \mathcal{L_3}&=\Lambda_2^4 G_3 [\Phi], \nonumber \\ \mathcal{L_4}&= M_p^2 G_4 R + \Lambda_2^4 G_{4,X}([\Phi]^2 - [\Phi^2]),\nonumber \\ \mathcal{L_5}&= M_p^2 G_5 G_{\mu \nu}\Phi^{\mu \nu} - \frac{1}{6}\Lambda_2^4 G_{5,X}([\Phi]^3 - 3[\Phi][\Phi^2] + 2[\Phi^3]), \nonumber \end{align}

where $G_2,G_3,G_4,G_5$ are functions of $\phi$ and $X = -\frac{1}{2}\nabla^\mu \phi \nabla_\mu \phi /\Lambda_2^4$, $\Phi^{\mu}_{ \ \nu}:= \nabla^\mu \nabla_\nu \phi/\Lambda_3^3$ and square brackets indicate the trace, e.g. $[\Phi^2] = \nabla^\mu \nabla_\nu \phi \nabla^\nu \nabla_\mu \phi/\Lambda_3^6$ and $,$ denote partial derivatives. We assume that $\langle \phi \rangle = 0$ on the background, bars will indicate evaluation at the background.

I am trying to derive the expression for the graviton propagator in this theory but it does not work out well. My idea was to identify the Lagrangian for $2$ gravitons (where $h = h^\mu_\mu$):

\begin{align} \mathcal{L}_{hh} &= \Lambda^4_2 \bar{G}_2\sqrt{-g} + \bar{G}_4 M_{\mathrm{pl}}^2 R \sqrt{-g}\\ &\approx \bar{G}_2 H_0^2 \Big(1+\frac{1}{2}h + \frac{1}{8}h^2 - \frac{1}{4}h_{\mu \nu}h^{\mu \nu}\Big) + \bar{G}_4 \Big(-\frac{1}{2}(\partial_\sigma h)(\partial_\mu h^{\mu \sigma}) + \frac{1}{2}(\partial_\nu h)(\partial^\nu h) + \frac{1}{2}(\partial_\nu h^{\mu \nu})(\partial^\sigma h_{\sigma \mu})+ \frac{1}{2}(\partial^\sigma h_{\sigma \nu})(\partial_\mu h^{\mu \nu}) - \frac{1}{2}(\partial_\nu h)(\partial_\mu h^{\mu \nu}) - \frac{1}{2}(\partial_\beta h^{\mu \nu})(\partial^\beta h_{\mu \nu})\Big). \end{align}

In ordinary GR one would add the gauge fixing term $-\bar{G}_4 (\partial_\nu h^{\mu \nu} - \frac{1}{2}\partial^\mu h)^2$ so that the last term becomes $\bar{G}_4(\frac{1}{4} (\partial_\nu h)^2 - \frac{1}{2} (\partial_\beta h^{\mu \nu})^2)$. However, in this case in the paper arXiv:1904.05874 they mention that the graviton propagator $\mathcal{P}^{\nu \beta}_{\mu \alpha}$ is found from ($\delta$ is the generalised Kronecker delta):

$$\frac{\bar{G}_4}{2}\mathcal{P}^{\nu \beta}_{\mu \alpha}(p) \delta^{\alpha \rho \mu^\prime}_{\beta \sigma \nu^\prime} p_\rho p^\sigma = -i \delta^{\mu^\prime}_\mu \delta^{\nu^\prime}_\nu.$$

My questions are: How do they get rid of the $\bar{G}_2$ term and what gauge fixing term has taken to arrive at this result for the definition of the propagator?