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))\right) -\frac{1}{8\pi G}\int d^3 x\ (\nabla\Phi(\mathbf{x},t))^2 \tag{1}$$$$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$$\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}$$