The standard textbook approach in Newtonian gravitational dynamics is to derive the particle dynamics using the particle Lagrangian:
$$L = T-V = \frac 1 2 m \dot x_u \dot x_u -m\phi(x_u)\tag{1}$$
With the corresponding Euler-Lagrange equation taking the variation w.r.t. $x$:
$$m\ddot x_u=\frac{ \partial\phi}{ \partial x_u}\tag{2}$$
or $$\textbf{F} = m \ \textbf{a};\tag{3}$$ and, derive the Poisson field equation using the Lagrangian density:
$$ \mathfrak L = -\rho \phi - \frac {1}{8\pi G} (\partial_u\phi)(\partial_u\phi)\tag{4}$$
which is substituted in the Euler-Lagrange equation taking the variation w.r.t. $\phi$:
$$\frac{\partial\mathfrak L}{\partial\phi}-\frac{\partial}{\partial x_u}(\frac{\partial\mathfrak L}{\partial\phi_u})=0\tag{5}$$
Yielding:
$$-\rho +\frac {1}{4\pi G}\nabla^2\phi=0 \ \ \ or \ \ \ \nabla^2\phi=4\pi G\rho.\tag{6}$$
I have several questions about this.
First, is it essential to use two different Lagrangians? What is the physical basis and meaning of having two Lagrangians? Is it a rigorous approach?
Second, is it possible to use one Lagrangian density that encompasses both. Something like:
$$\mathfrak L = -\rho(\phi-\frac 1 2 m \dot x_u \dot x_u)-\frac {1}{8\pi G} (\partial_u\phi)(\partial_u\phi)\tag{7}$$
and taking variation w.r.t. to $x$ and $\phi$. This sort of works if, in the variation w.r.t. $x$, the
$$\frac{\partial}{\partial x_u}\tag{8}$$
only acts on the first term in parenthesis, giving $F = ma$, but this does not seem right because both $\rho$ and $\partial_u\phi$ are also functions of $x$ and would add additional terms.