given by equation (\[D3MG\]) in the appendix \[polyderiv\]. It should be also noted that equations (\[growth1st\]), (\[eq:D2sources\]) and f(\[D3MG\]) can be either solved by inverting the linear operator $\left(\mathcal{\hat{T}}-A(k)\right)$ using its Green function, as done in [@Matsubara:2015ipa; @Aviles:2017aor], or by numerically solving the corresponding differential equations. Even though both methods give results that agree with each other well, we chose to proceed with the latter because it is computationally faster. The differential equations were solved using a $5^{th}$ order Runge-Kutta scheme, implemented in $\it{Mathematica}$ [@Mathematica].

We finish this section by showing the particular expressions for the mass terms in (\[interaction\]) and the sources in (\[eq:sources\]) for the $f(R)$ and $n$DGP models we study. For the $f(R)$ Hu-Sawicki model, the Brans-Dicke scalar is simply $\omega_{BD}=0$, while the mass terms are given by the expansion [@PhysRevD.79.123512]: $$\label{eq:fRmasses}
M_{n} = \frac{d^{n}\bar{R}(f_{R_0})}{df_{R_0}^{n}},$$ which, using (\[fRHu\]), gives [@Aviles:2017aor]: $$\begin{aligned}\label{eq:fRmasses2}
M_{1}(a) &= \frac{3 H_0^2}{2 |f_{R_0}|}\frac{\left(\Omega_m a^{-3}+4 \Omega_{\Lambda}\right)^3}{\left(\Omega_m+4 \Omega_{\Lambda}\right)^2},\\
M_{2}(a) &= \frac{9 H_0^2}{4 |f_{R_0}|}\frac{\left(\Omega_m a^{-3}+4 \Omega_{\Lambda}\right)^5}{\left(\Omega_m+4 \Omega_{\Lambda}\right)^4},\\
M_{3}(a) &= \frac{45 H_0^2}{8 |f_{R_0}|}\frac{\left(\Omega_m a^{-3}+4 \Omega_{\Lambda}\right)^7}{\left(\Omega_m+4 \Omega_{\Lambda}\right)^6}.
\end{aligned}$$ In the case of the $n$DGP braneworld model, a similar procedure, informed by (\[eq:poisson\_dgp\]), gives the relevant expressions [@PhysRevD.79.123512; @Aviles:2018saf]: $$\begin{aligned}\label{eq:nDGPmasses}
&M_{1}(a) = 0,\\
&M_{2}(\bold{k}_1,\bold{k}_2,a) = 2\frac{n^2}{H_0^2 a^4}\left(k_1^2k_2^2-(\bold{k_1}\cdot\bold{k_2})^2\right),\\
&M_{3}(\bold{k}_1,\bold{k}_2,\bold{k}_3,a) = 18\frac{n^2 