Derivation of gravitational dynamics using Lagrangian? 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.
 A: *

*It is important to mention that in this model the point particle masses $m_i$ serve as sources
$$ \rho({\bf x},t)=\sum_i m_i \delta^3({\bf x}-{\bf x}_i(t)) $$
for the gravitational potential field $\phi$.


*The full Lagrangian
$$L=T_1-V_2-V_3$$ for both the point particle positions ${\bf x}_i(t)$ and the field $\phi({\bf x},t)$ consists of

*

*a kinetic term $T_1=\frac{1}{2}\sum_i m_i \dot{\bf x}_i^2(t)$.


*an interaction term $V_2=\sum_i m_i \phi({\bf x}_i(t),t)$.


*a potential energy term $V_3=\int_{\mathbb{R}^3}\! d^3{\bf x}~{\cal V}_3({\bf x},t)$, where ${\cal V}_3=\frac{1}{8\pi G}(\nabla\phi)^2$.




*It is possible to integrate out the $\phi$ field to obtain a purely point-mechanical Lagrangian, cf. my related Phys.SE answer here.
A: The first Lagrangian is the Lagrangian for the motion of a point particle in  a Newtonian gravitational potential (or field), whereas the second is the Lagrangian for the gravitational potential (or field), the latter which eventually leads to the Newtonian field equation. So, actually there is no real relationship. The only connection  one might think of is that it is about Newtonian gravity. But the equation of motion (EOM) of a particle is fundamentally different from the field equation (FE). In the first (EOM) the motion of a single particle is searched for in a given known Newtonian gravitational potential (or field) --- a problem of mechanics, whereas in the second, an unknown Newtonian gravitational potential is searched for in a given mass density distribution --- a problem of field theory.
Also, the combined Lagrangian makes not much sense: in particular in the first phrase of the post $\phi(x)$ is considered as potential, whereas in the "combined" Lagrangian, it has the function of potential energy.
Lagrangians with different variables to be varied, actually exist.
But I wonder if a Lagrangian where, one varies in the first part of the coordinates and their time derivatives and in the second part of the field components, and their spacial derivatives, makes any sense?
By the way, the correct EOM (even in Newtonian gravity) should be
$$\ddot{x}_u = -\frac{\partial \phi}{\partial x}$$
So, according to Newton, the coincidental equality of inertial and gravitational mass, can have the mass factored out completely.
This can indeed be found by varying the coordinates and their time derivates in the first Lagrangian:
$$L = \frac{1}{2}m \dot{x}_u \dot{x}_u - m \phi(x_u)$$
where $\phi$ being the potential that is not equal to the potential energy $\phi \neq V$.
A: 
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?

No, it is not essential to use two different Lagrangians.
In my opinion it is just a conceptional step before merging
these two partial Lagrangians into a single common Lagrangian.

Second, is it possible to use one Lagrangian density that encompasses both.

Yes, you can derive the equations of motion for
$N$ particles (with masses $m_i$ at positions $\mathbf{x}_i(t)$, for $i=1, ... N$)
and for the gravitational field ($\Phi(\mathbf{x},t)$, for $\mathbf{x}\in\mathbb{R}^3$)
from a common Lagrangian.
However, we need to carefully distinguish between
the discrete positions $\mathbf{x}_i$ of the particles
and the continuous position $\mathbf{x}$ of the field.
The complete Lagrangian looks like this:
$$L[\mathbf{x}_i,\dot{\mathbf{x}}_i,\Phi,\nabla\Phi]=
  \sum_{i=1}^N \left(\frac{1}{2}m_i\dot{\mathbf{x}}_i(t)^2
  -m_i \Phi(\mathbf{x}_i(t),t)\right)
  -\frac{1}{8\pi G}\int d^3 x\ (\nabla\Phi(\mathbf{x},t))^2  \tag{1}$$
It has 3 components:

*

*a kinetic energy for each particle,

*an interaction energy between each particle and the gravitational field,

*the energy of the gravitational field.

We can define the density field made up by the particles as
$$\rho(\mathbf{x},t)=\sum_{i=1}^N m_i\ \delta^3(\mathbf{x}-\mathbf{x}_i(t))  \tag{2}$$
which has Dirac-like peaks where the particles are.
Using this density we can rewrite the Lagrangian (1) as:
$$L[\mathbf{x}_i,\dot{\mathbf{x}}_i,\Phi,\nabla\Phi]=
  \sum_{i=1}^N \frac{1}{2}m_i\dot{\mathbf{x}}_i(t)^2
  +\int d^3 x\underbrace{\left(
  -\rho(\mathbf{x},t)\Phi(\mathbf{x},t)-\frac{1}{8\pi G}(\nabla\Phi(\mathbf{x},t))^2
  \right)}_{\mathcal{L[\Phi,\nabla\Phi]}}  \tag{3}$$
From Lagrangian (1) or (3) we get the equations of motion.
From Lagrangian (1) the Euler-Lagrange equations
with respect to $\mathbf{x}_i(t)$ give Newton's law:
$$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\mathbf{x}}_i}\right)
  =\frac{\partial L}{\partial\mathbf{x}_i}$$
$$m_i\ddot{\mathbf{x}}_i=-m_i\nabla\Phi(\mathbf{x}_i)  \tag{4}$$
And from Lagrangian (3) the Euler-Lagrange equation
with respect to $\Phi(\mathbf{x},t)$
gives Poisson's equation:
$$\nabla\left(\frac{\partial{\mathcal{L}}}{\partial(\nabla\Phi)}\right)
  =\frac{\partial\mathcal{L}}{\partial\Phi}$$
$$-\frac{1}{4\pi G}\nabla^2\Phi=-\rho  \tag{5}$$
