Energy density of dust in f(R) gravity

Question

Usually during the calculation of the Einstein-Hilbert action, we consider the integral of the variation of the Ricci tensor to be zero as the variation of it, at the boundary was assumed to be zero, as:

$$-\frac{1}{16\kappa{\pi}}\int{d^{4}x[\delta\sqrt{|g|}(R+\Lambda)+\sqrt{|g|}(\delta{g^{\mu{\nu}}})R_{\mu{\nu}}+\sqrt{|g|}{g^{\mu{\nu}}}(\delta{R_{\mu{\nu}}})-8\pi{\kappa}}\sqrt{|g|}\delta{g^{\mu{\nu}}}T_{\mu{\nu}}]$$ 
$\kappa=\frac{G}{c^{4}}$

Here:
$$g^{\mu{\nu}}\delta{R_{\mu{\nu}}}=D_{\mu}\left[g^{\alpha{\beta}}\Gamma^{\mu}_{\alpha{\beta}}-g^{\alpha{\mu}}\Gamma^{\beta}_{\alpha{\beta}}\right]=D_{\mu}\delta{U^{\mu}}$$
$$\int_{\mathcal{M}}{d^{4}x\sqrt{|g|}D_{\mu}\delta{U^{\mu}}=0}$$

Now, if we are to continue calculating this integral with the vanishing boundary, we obtain the following expression for the variation of the Ricci tensor:
$$g^{\mu{\nu}}\delta{R_{\mu{\nu}}}=g_{\mu{\nu}}\Box\delta{g^{\mu{\nu}}}-D_{\mu}D_{\nu}(\delta{g^{\mu{\nu}}})$$
Using the above result in the action, we obtain:
$$R_{\mu{\nu}}-\frac{1}{2}g_{\mu{\nu}}R-\frac{1}{2}g_{\mu{\nu}}\Lambda+\left[g_{\mu{\nu}}\Box-D_{\mu}D_{\nu}\right]=8\pi{\kappa}T_{\mu{\nu}}$$
Now, under the vacuum conditions, $\Lambda=0$, and  $T_{\mu{\nu}}=0$, we obtain a value for the Ricci scalar as:
$$R+\delta^{\mu}_{\mu}\Box=\frac{1}{2}\delta^{\mu}_{\mu}R+g^{\mu{\nu}}D_{\mu}D_{{\nu}}$$
$$R=3\Box$$
Now, in the case of Dust (in co-moving coordinates), $T_{\mu{\nu}}=\rho{g_{0\mu}g_{0\nu}}$, we obtain the following:
$$R-\frac{1}{2}\delta^{\mu}_{\mu}R=8\pi{\kappa}T+2\Lambda-4\Box$$
This yields:
$$T_{\mu{\nu}}g^{\mu{\nu}}=T=\rho=\frac{\Box-2\Lambda}{8\pi{\kappa}}$$
But we know that the density in the dust model has the following form:
$$\rho=\sum_{q=1}^{N}\gamma{m_{q}}\delta^{(3)}[\vec{x}-\vec{z_{q}}(S)]\frac{1}{\sqrt{|g|}}$$
This implies:
$$\Lambda=\frac{1}{2\sqrt{|g|}}\left[\Box{\sqrt{|det(g_{\mu{\nu}})|}}-8\pi{\kappa}\sum_{q=1}^{N}\gamma{m_{q}}\delta^{(3)}[\vec{x}-\vec{z_{q}}(S)]\right]$$

If we replace $1/G$ with a scalar field $\phi$, which varies from one place to another, then we obtain the following form:
Note: the D'alembertian for a scalar field is: $\Box=\frac{1}{\sqrt{|g|}}\partial_{\mu}(\sqrt{|g|}\partial^{\mu})$, hence, $\kappa=G/c^{4}=1/(c^{4}\phi)$.
From Brans-Dicke theory:$$\Box{\phi}=\frac{8\pi{T}}{(3+2\omega)}$$
where $\omega$ is the dimensionless Dicke coupling constant.
This yields: >$$\Lambda\approx\frac{4\pi{T}}{\phi(3+2\omega)}$$

Under strong coupling, i.e. $\omega>-1.5$, $\Lambda$ is positive ($\Lambda\approx\frac{2\pi{\rho}}{3\phi}$ when $\omega=1.5$). Under weak coupling, i.e. $\omega<-1.5$, $\Lambda$ is negative. Is my derivation correct? What can be the possible implications of this? I am aware that in f(R) gravity, in Einstein field equations when linearized on a de Sitter background gives rise to wave equations for the de Sitter covariant field $\phi_{\mu{\nu}}$: $$(\Box-2\Lambda)\phi\approx{8\pi{\kappa}T_{\mu{\nu}}}$$

When $\Lambda$ is negative, we can observe that the equation takes the form of Klein-Gordon equations for a massive spin-zero field.

Comments

Somehow it is not clear to me what do you mean by expressions such as $R = 3 \Box$. On the left-hand-side I see a scalar and on the right a differential operator. Are you assuming a space of test functions or tensors? What is it and what is its meaning?

---

I am not sure about that part too, but it seems to make sense when I use it in Einstein field equations when linearized on a de Sitter background, which gives rise to wave equations for the de Sitter covariant field ϕμν. ​Also the d'alembertian would change its form for a scalar field as when I replace (1/G) with a scalar field ϕ​.