Lagrangian density for Lorentz force of continuous charge distribution in external field? It's frequently an exercise to derive the Lorentz force law for a particle with charge $q$ in an external electromagnetic field given by the following Lagrangian:
$$L = -mc^2\sqrt{1-\frac{\dot{r}^2}{c^2}} - q \phi + q \mathbf{\dot{r}} \cdot \mathbf{A}$$
Which leads to the relativistic Lorentz force law:
$$\mathbf{F} = \frac{d}{dt} \Bigg(\frac{m \mathbf{\dot{r}}}{\sqrt{1-\frac{\dot{r}^2}{c^2}}} \Bigg) = q(\mathbf{E} + \mathbf{\dot{r}} \times \mathbf{B})$$
For continuous distributions, we have:
$$\mathbf{f} = \rho \mathbf{E} + \mathbf{J} \times \mathbf{B}$$
I'm trying to find the corresponding Lagrangian density which results in this force. I know if the charge distribution is treated as the source you can use the standard Lagrangian Density for electromagnetism but doing so will not give you the Lorentz Force Equation. However in my specific case I'm ignoring the self-field of the charge distribution, the fields are purely external, I don't need the electromagnetism Lagrangian density for my problem. Naively one might replace all instances of $m$ with a mass density term, $\rho_m$, and all instances of $q$ with a charge density term, $\rho_q$, where $\mathbf{J} = \rho_q \mathbf{\dot{r}}$. 
$$\mathcal{L} = -\rho_mc^2\sqrt{1-\frac{v(\mathbf{r},t)^2}{c^2}} - \rho_q \phi + \rho_q \mathbf{v}(\mathbf{r},t) \cdot \mathbf{A} = -\rho_mc^2\sqrt{1-\frac{v(\mathbf{r},t)^2}{c^2}} - \rho_q \phi + \mathbf{J} \cdot \mathbf{A}$$
However the densities are also a function of the coordinates and moreover the mass densities and charge densities are related to each other in some unknown way. If we assume all of our particles are electrons then we can scale the densities by the electron mass and charge. If we take the variation of this Lagrangian density with respect to the charge density then we get the following:
$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\rho}} + \frac{d}{dx}\frac{\partial \mathcal{L}}{\partial \rho_x} + \frac{d}{dy}\frac{\partial \mathcal{L}}{\partial \rho_y} + \frac{d}{dz}\frac{\partial \mathcal{L}}{\partial \rho_z}= \frac{\partial \mathcal{L}}{\partial \rho}$$
The LHS is clearly zero since we have no dependence on derivatives of the density in our Lagrangian density. The RHS just gives us:
$$0 = -\frac{mc^2}{e}\sqrt{1-\frac{\dot{r}^2}{c^2}} - \phi + \mathbf{v}(\mathbf{r},t) \cdot \mathbf{A}$$
So clearly this Lagrangian density is not correct or we should not be varying with respect to the density. Similarly one can take the variation with respect to the velocity field but this also does not result in the correct equation. I feel like I'm having a fundamental misunderstanding here but I can't find any reference that works through this. What is the correct Lagrangian Density? What is the correct quantity to vary the action?
 A: Comments to the question (v3):


*

*OP is essentially asking about the Lagrangian field-theoretic formulation of a relativistic fluid in an external electromagnetic background $A_{\mu}$. 

*Fluid dynamics have both a Lagrangian and an Eulerian picture. (Note that the word Lagrangian is used in two different meanings.) In the relativistic context, there is also the issue of a manifestly Lorentz-invariant formulation.

*Here is the simplest Lagrangian Lagrangian relativistic formulation (with Lorentz symmetry but without manifest Lorentz symmetry). This essentially boils down to replacing discrete sums in point mechanics with continuous integrals in field theory. We put the speed of light $c=1$ to one for simplicity. The 3-position field ${\bf r}:\mathbb{R}^3\times \mathbb{R}\to\mathbb{R}^3$ depends on a continuous labelling variable ${\bf a}\in\mathbb{R}^3$  and time $t\in\mathbb{R}$.
The action becomes 
$$ S[{\bf r}] ~=~ \left. \int \! dt~ d^3a~ {\cal L}({\bf r}({\bf a},t),{\bf v}({\bf a},t),{\bf a} ,t)  \right|_{{\bf v}=\dot{\bf r}} \quad,\tag{1} $$
where the Lagrangian density is
$$ {\cal L}({\bf r}({\bf a},t),{\bf v}({\bf a},t),{\bf a} ,t) 
~=~ -\frac{\mu({\bf a})}{\gamma({\bf v}({\bf a},t))} 
-\rho({\bf a})~\phi({\bf r}({\bf a},t),t)+\rho({\bf a})~{\bf v}({\bf a},t)\cdot{\bf A}({\bf r}({\bf a},t),t) ,\tag{2}$$
and where 
$$\gamma({\bf v})~:=~\frac{1}{\sqrt{1-{\bf v}^2}}\tag{3}$$
is the gamma factor. The rest mass $m(\Omega)$  and the charge $Q(\Omega)$ in a region $\Omega\subseteq \mathbb{R}^3$ of labelling 3-space is given by
$$ m(\Omega) ~=~ \int_{\Omega} \!  d^3a~\mu({\bf a})\quad\text{and}\quad Q(\Omega) ~=~ \int_{\Omega} \!  d^3a~\rho({\bf a}),\tag{4}$$
respectively.

*The Eulerian formulation is more complicated (already in the non-relativistic case) due to labelling gauge-symmetry, cf. e.g. this Phys.SE post and references therein. If time permits, I might write down the Eulerian formulation explicitly in a future update.
