\beta(a)}{H_0^2 a^6 4 \pi G\bar{\rho}_m}\Biggl((\bold{k_1}\cdot\bold{k_2})k_1^2k_3^2+2(\bold{k_1}\cdot\bold{k_2})^2k_3^2\\
& -2(\bold{k_1}\cdot\bold{k_2})(\bold{k_1}\cdot\bold{k_3})(\bold{k_2}\cdot\bold{k_3})-(\bold{k_1}\cdot\bold{k_2})(\bold{k_1}\cdot\bold{k_3})^2\Biggr), \\
\end{aligned}$$ with $\beta$ defined in (\[eq:dgpbeta\]). It is interesting to notice that, even though the interaction term in (\[eq:poisson\_dgp\]) contains only second order derivatives, in (\[eq:nDGPmasses\]) a $3^{rd}$ order mass contribution is now present, that arises when transforming the Eulerian derivatives to the Lagrangian space through (\[Jacobian\]).

### 2-point statistics for biased tracers in GR {#sec:twopointGR}

The perturbative theory of galaxy clustering [@Desjacques:2016bnm], which aims to describe the statistics of biased tracers in the quasi-linear regime, consists of a perturbative description for the evolution of the underlying dark matter density field, combined with an analytical description for the bias parameters at each given order. In the case of cold dark matter, the calculation of the 2-point statistics, even in LPT, is a more straightforward process, since one just needs to plug the ${\bf q}$-space Lagrangian overdensity, mapped to the Eulerian frame through (\[delJac\]), into the common expressions for the autocorrelation function: $$\label{xi}
\xi(r) = \langle \delta_m(\bold{x})\delta_m(\bold{x+r})\rangle$$ and its Fourier space counterpart, the matter power spectrum $$\label{Pk}
\left(2\pi\right)^3\delta_D(\bold{k}+\bold{k'})P(k) = \langle \tilde{\delta}_m(\bold{k})\tilde{\delta}_m(\bold{k'})\rangle.$$ When studying biased tracers, like for example dark matter halos, we need an analytical model to describe their statistical prevalance with respect 