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.