Rigorous formulation of electromagnetic field theory for a system of moving charges (non relativistic) using distribution theory

Question

Suppose we have a system of $n$ charged particles with trajectories given by a paths $x_j:\mathbb{R}\to\mathbb{R}^3$ then the Maxwell equations for this system are given first by defining the charge and current densities $$\begin{align} \rho(t,x)&=\sum\limits_{j=1}^ne_j\delta(x-x_j(t))\\ \mathbf{J}(t,x)&=\sum\limits_{j=1}^ne_j\dot{x}_j(t)\delta(x-x_j(t)) \end{align}$$ With $\delta$ being the usual Dirac delta distribution. And then for vector fields $\mathbf{E},\mathbf{B}:\mathbb{R}\times\mathbb{R}^3\to \mathbb{R}^3$ we impose that they satisfy $$\begin{align} \nabla\cdot\textbf{E}(t,x)&=4\pi\rho(t,x)\\ \nabla\cdot\textbf{B}(t,x)&=0\\ \nabla\times\textbf{E}(t,x)&=-\frac{1}{c}\frac{\partial \textbf{B}(t,x)}{\partial t}\\ \nabla\times\textbf{B}(t,x)&= \frac{1}{c}\left( 4 \pi\mathbf{J}(t,x)+\frac{\partial\mathbf{E}(t,x)}{\partial t}\right) \end{align}$$ In field theory, it is much more common to work with the potentials, these are constructed by first using the second equation listed here and Helmholtz theorem to say that a vector potential $\textbf{A}$ must exist. Then by manipulating the third equation we you get an irrotational vector field and this lets you construct a scalar potential $\phi$. With this a non relativistic Lagrangian is constructed by considering the kinetic energy of the particles, the energy from the interaction of the particles with the fields and the energy of the fields itself as $$\begin{align} \nonumber \mathcal{L}\left(x,\dot{x},\phi,D\phi,\mathbf{A},D\mathbf{A}\right)&=\sum\limits_{j=1}^n\frac{m_j}{2}|\dot{x}_j|^2-\int\left(\rho\phi-\frac{1}{c}\mathbf{J}\cdot\mathbf{A}-\frac{1}{8\pi}(|\mathbf{E}|^2-|\mathbf{B}|^2)\right)dx\\ \nonumber &=\sum\limits_{j=1}^n\left(\frac{m_j}{2}|\dot{x}_j|^2-e_j\phi(t,x_j(t))+\frac{e_j}{c}\dot{x(t)}\cdot\mathbf{A}(t,x_j(t))\right)\\ &\;+\frac{1}{8\pi}\int\left(|-\nabla\phi(t,x)-\frac{1}{c}\frac{\partial}{\partial t}\mathbf{A}(t,x)|^2-|\nabla\times\mathbf{A}|^2\right)dx \end{align}$$ Now I know that from how I have constructed the theory, plenty of mathematical problems arise. I do understand, or at least I am familiar with distribution theory but I'm not too sure how to use it to rigorously justify all these steps that are present in many physics texts. Though I think the questions reduce to

1. The charge and current densities are time dependent and so they are curves in the space of tempered distributions. What does it mean here for the curve to be smooth? I understand the distributional derivative is well defined, but what does it mean to take the time derivative of the densities? I suppose it reduces to constraints in the paths of the charges which need to have some sort of second derivative (In a traditional or weak sense maybe).
2. The fields, at least, must be time dependent vector-distributions so the question 1 applies, but also Helmholtz decomposition theorem must apply to construct the potentials. I am not aware of any distributional versión of this theorem so how is it justified?
3. The Lagrangian uses the formal properties of the Dirac delta distribution by integrating the densities against the potentials, however the Dirac delta is only defined for Schwartz functions and in general the potentials are not in this space, so is there a way to define the energy do to the interaction between the fields and the charges in a general way?
4. Also, the energy due to the fields appear as an integral that it is well known not to converge for even for a single static charge so what does this Lagrangian actually mean? And while the Lagrangian density is well defined, in order to apply the Euler-Lagrange equations again I don't know if there is a distributional formulation of these types of techniques or if the fields must at least belong to some Sobolev space (not sure considering the case of a single static charge)