Clarification on showing covariance of density expressions under Lorentz boost

For a collection of point charges, the charge density is defined as

$$$$\rho(\textbf{x},t) = \sum_{k=1}^{N} q_k \delta^{(3)}\left[\textbf{x}-\textbf{x}_k(t)\right]$$$$

while the current density is defined as

$$$$\textbf{J}(\textbf{x},t) = \sum_{k=1}^{N} q_k \textbf{u}_k(t) \delta^{(3)}\left[\textbf{x}-\textbf{x}_k(t)\right]$$$$

To show that the expressions above are covariant, I assumed that I need to transform the above expressions under Lorentz boost and show that the continuity equation holds.

The transformation is as such

$$$$\rho'(\textbf{x}',t') = \sum_{k=1}^{N} q_k \delta^{(3)}\left[\textbf{x}'-\textbf{x}'_k(t')\right]-\frac{\gamma}{c^2}\textbf{v}\cdot\sum_{k=1}^{N} q_k \textbf{u}_k(t') \delta^{(3)}\left[\textbf{x}'-\textbf{x}'_k(t')\right]$$$$

$$$$\textbf{J}'_{\parallel}(\textbf{x}',t') = \gamma \sum_{k=1}^{N} q_k \textbf{u}_{k,\,\parallel}(t')\delta^{(3)}\left[\textbf{x}'-\textbf{x}'_k(t')\right]+\gamma \textbf{v} \sum_{k=1}^{N} q_k \delta^{(3)}\left[\textbf{x}'-\textbf{x}'_k(t')\right]$$$$

$$$$\textbf{J}'_{\perp}(\textbf{x}',t') = \sum_{k=1}^{N} q_k \textbf{u}_{k,\,\perp}(t')\delta^{(3)}\left[\textbf{x}'-\textbf{x}'_k(t')\right]$$$$

The continuity equation is given as such

$$$$\frac{\partial}{\partial t'}\rho'(\textbf{x}',t')=-\nabla'\cdot\textbf{J}'(\textbf{x}',t')$$$$

However, when I perform the explicit calculations, the steps were quite confusing and tedious. I could not show that the continuity equation holds too.

Is there a more efficient way to show the expressions for charge and current density are covariant under Lorentz boost with explicit calculations?

I think you are making your life difficult with considering collections of charges. Does one loose anything by establishing relations for a single point-charge? I don't think so.

Next, you have to be very careful about changing reference frames with 3-delta functions. Delta-functions are 'per-volume', so when you start messing about with reference frames one has to compensate the length contraction effects. The transformation rule would be $$\delta^{(3)}\to \gamma\delta^{(3)}$$, where $$\gamma$$ is the Lorentz factor, but I don't really like this approach since it misses the connection between space and time IMHO.

I think it is easier to first prove that four-current $$J^\mu=\left(c\rho,\,\mathbf{J}\right)^\mu$$ is a four-vector and then charge conservation in any reference frame becomes zero-divergence statement (follows from the anti-symmetry of the electromagnetic tensor and Covariant form of Ampere's law):

$$\partial_\mu J^\mu=\partial_{ct}\left(c\rho\right)+\boldsymbol{\nabla}.\mathbf{J}=0$$

Where $$c$$ is the speed of light

Showing that four-current of a single charge is a four-vector

As you said, the charge density of a single charge is:

$$\rho\left(t, \mathbf{r}\right)=Q\delta^{(3)}\left(\mathbf{r}-\mathbf{r}_Q\left(t\right)\right)$$

Where $$Q$$ is the charge, and $$\mathbf{r}_Q$$ is the position of the charge. What I can do next is convert the 3d delta to 4d delta ($$c$$ is the speed of light)

$$\rho\left(t, \mathbf{r}\right)=Q\int cdt'\,\delta\left(ct'-t\right)\delta^{(3)}\left(\mathbf{r}-\mathbf{r}_Q\left(t'\right)\right)=\frac{Q}{c}\int cdt'\,\,c\,\delta^{(4)}\left(x-x_Q\right)$$

Where $$x^\mu=\left(ct,\,\mathbf{r}\right)^\mu$$. The nice thing about this expression is that it now is based on 4-delta function therefore when you change reference frames all you have to multiply the $$\delta^{(4)}$$ by, is magnitude of the Jacobian of the coordinate transformation. If you stick to reflection, rotations and boosts that Jacobian magnitude will be 1, so you can treat $$\delta^{(4)}$$ as a 'normal function' (only works in Cartesian coordinates!).

Next, lets deal with integration variable $$t'$$, which is the temporal coorindate in your lab frame. This is inconvenient. It is easy to show that $$dt=\gamma d\tau$$, where $$\tau$$ is the proper time and $$\gamma$$ is the Lorentz factor. Therefore, swapping the integration variable we find:

$$c\rho\left(t, \mathbf{r}\right)=Q\int cd\tau\,\,\gamma c\,\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)$$

You can repeat similar machinations with your current density to arrive at:

$$\mathbf{J}\left(t, \mathbf{r}\right)=Q\int cd\tau\,\,\gamma \mathbf{v}\left(\tau\right)\,\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)$$

Finally note that in the lab-frame the four-velocity of the charge is $$u^\mu=\left(\gamma c,\,\gamma\mathbf{v}\right)^\mu$$, which is a four-vector. Therefore:

$$J^\mu = Q\int c d\tau\,\,u^\mu\left(\tau\right)\,\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)=\left(c\rho,\,\mathbf{J}\right)^\mu$$

is a four-vector (since integration with respect to proper time is frame-independent, so is definition of four-velocity, so is 4d delta function).

Divergence of four-current

It is interesting to try to take the divergence of the single-charge four-current

\begin{align} \partial_\mu J^\mu &= Q\int c d\tau\,\,u^\mu\left(\tau\right)\,\partial_\mu\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)=-Q\int c d\tau\,\,u^\mu\left(\tau\right)\,\frac{\partial}{\partial x_Q^\mu}\delta^{(4)}\left(x-x_Q\left(\tau\right)\right) \\ &=-Q\int c d\tau\,\,\frac{dx_q^\mu}{d\tau}\,\frac{\partial}{\partial x_Q^\mu}\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)=-cQ\int d\tau\,\,\frac{d}{d\tau}\delta^{(4)}\left(x-x_Q\left(\tau\right)\right)\\ &=0 \end{align}

Expression vanishes because the observer has to be essentially at the beginning or the end of time for integral not be be zero.