Proof 4-current is indeed a 4-vector for a system of particles

$rho dV = dq$

$dq dx^{i} = rho dV dx^{i} = rho d^{4}x frac{dx^{i}}{dt}$

In L.H.S. $dq$ is invariant and $dx^{i}$ is a 4- vector hence R.H.S. should be a 4-vector too but we know $d^{4}x$ is invariant. Hence $rho frac{dx^{i}}{dt}$ is a 4-vector which is a 4-current.

The definition of 4-tuple (charge density, current density) is of course not enough to find how its values transform to another reference frame. There is a common idea that the continuity equation

$$ partial_mu j^mu = 0 $$ by itself implies $j^mu$ transforms as a four-vector. It does not, because even if there is a solution family (having one solution for every reference frame) $j_1^mu$ whose members transform via Lorentz transformation, there are other solution families that do not, for example $$ j_2^mu = j_1^mu + delta^{mu 0} C $$ where $C$ is a constant independent of reference frame. This solution has shifted charge density in all frames by a constant value and hence the transform that connects $j_2$ to $j_2'$ cannot be a linear transformation.

So, we need some other assumption / knowledge about this 4-tuple to derive its transformation properties.

One way is to use the knowledge that this 4-tuple obeys Maxwell's equations

$$ partial_{nu} F^{numu} = mu_0 j^mu $$ and the knowledge that matrix $F$, defined as 4x4 matrix of electric and magnetic components, transforms as $X^mu Y^nu$ via Lorentz transformation, where $X,Y$ are four-vectors. This is the historical way, because people noticed electric and magnetic field have to transform via Lorentz transformation (the simplest possible transformation) in order to preserve the form of Maxwell's equations in all inertial frames. Assuming that, the left-hand-side transforms as four-vector, and hence the right-hand side does as well. This is quite a general argument and it does not matter how charge or current is distributed in space. So it is valid for your system of point particles as well.

The other way is to assume electric charge density $rho$ and current density $mathbf j$ transform to other frames like particle number density $n$ and particle current density $nmathbf v$ for systems of conserved and invariant number of particles. This is natural if we imagine charge and current in some region of space to be formed by lots of point particles with the same velocity. That is, for any point of spacetime $P$ in the region, there is a frame of reference $S'$ where charged particles at spacetime point $P$ do not move and charge density there is the lowest for all frames of reference. Let this lowest value (proper density) be denoted $rho_0$, and the current 4-tuple there is

$$ J'(P) = [c rho_0, 0, 0, 0 ]~~~text{(no current, charges are at rest}). $$ In other frame $S$, moving with velocity $-mathbf V$ with respect to frame $S'$, velocity of an element of charged matter that is near $P$ is $mathbf V$ and current 4-tuple at the same spacetime point $P$ has values

$$ J(P) = [c gamma rho_0, gamma rho_0 mathbf V]. $$

Since $rho_0$ is invariant, we can put it as a prefactor and write $$ J(P) = rho_0 [cgamma, gamma mathbf V]. $$ which is invariant times four-velocity of the element.

This expression of $J(P)$ is valid for all possible velocities $mathbf V$. This expression is invariant value times a four-vector. Hence $J(P)$ transforms as four-vector. This argument is made for finite densities $rho_0$, but it can be shown the same idea works for delta distribution of charge and current as well.

Thus single particle current density 4-tuple is four-vector. Since Lorentz transformation is linear operation, any sum of quantities that transform as four-vector also transforms as four-vector. So also current density of a set of point particles transforms as four-vector.

This is a question that I wish more undergrads would ask! It’s a really important calculation, because it really requires you to “go back to the fundamentals” when you answer it.

First, let me remove an unnecessary complication you are adding: the Lorentz boost is linear and so if we solve this problem for one particle we will solve it by superposition for $N$ particles. You only need to handle one at a time. Furthermore we know that we can choose our coordinates arbitrarily so let us set our one particle moving in the xy-plane through the origin and then boost in the x-direction, that makes this essentially only a 2-dimensional problem. (If you do not like that multiply through everything by a $delta(z).$)

Let me do as is my custom and define the relativistic time coordinate $w = ct$ as well, with the charged particle's velocity vector $mathbf u$ transformed into some vector $veceta = (eta_x, eta_y) = mathbf u/c.$ In this language we would write $$rho(w, x, y) = q~delta(x - eta_x w)delta(y - eta_y w),mathbf J(w, x, y) = veceta~crho(w, x, y).$$ We have our 2D setup in the unprimed coordinates, great.

Now comes the part where we “go back to the fundamentals” and ask, what does the transformed charge density $rho'(w', x', y')$ really mean? We have to go back to the definition, it is the charge contained in the box formed by the product of open intervals $(x', x' + delta x') times (y', y' + delta y')$ at the time $w'$, divided by $delta x'~delta y'$. We can then transform that primed-stationary box back to the unprimed coordinates and we see that it is moving across space, and it is length-contracted. The length contraction turns out to be kind of gorgeously irrelevant here due to the fact that we're dealing with a point charge, as is the fact that “simultaneous” in the primed frame is not “simultaneous” across this square in the unprimed frame: we anticipate that the charge “should be” a Lorentz-scalar so we just need to find the time when the box hits our point charge.

So the box is a set of points $(x, y) in (gamma~(x' + beta w'), gamma~(x' + beta w') + delta x'/gamma) times (y', y' + delta y')$ at time $w = gamma~(w' + beta x')$. And so the point particle described by $rho$ is inside this volume only when $$x' + beta w' = eta_x (w' + beta x') y' = eta_y ~gamma~(w' + beta x'). $$ This in turn implies that the point charge is seen traveling in the primed coordinates at the point $$x' = frac{eta_x -beta}{1 - eta_x beta}~w', y' = frac{gamma^{-1}~eta_y} {1 - eta_x beta}~ w',$$ and so if it is to still be a point charge with magnitude $q$ (that is, if we assume that charge is a Lorentz scalar) then the proper expression for $rho', mathbf J'$ must be just, $$rho'(w', x', y') = q~deltaleft(x' - frac{eta_x -beta}{1 - eta_x beta}~w'right) ~deltaleft(y' - frac{gamma^{-1}~eta_y} {1 - eta_x beta}~ w'right),J_x' = frac{eta_x - beta}{1 - eta_x beta} ~crho', J_y' = frac{gamma^{-1}~eta_y} {1 - eta_x beta}~ crho'.$$ Now compare this with the naive transformation of the fields we have, $$ begin{align} rho(w', x', y') &= q~deltabig(gamma (x' + beta w') - eta_x gamma (w' + beta x')big) deltabig(y' - eta_y gamma (w' + beta x')big) &=frac{q}{gamma~(1-betaeta_x)} deltaleft(x' - frac{eta_x - beta}{1 - betaeta_x}right)deltaleft(y' - eta_y frac{gamma^{-1}~eta_y} {1 - eta_x beta}~ w'right) &= frac{rho'}{gamma~(1 - beta eta_x)}, end{align}$$and note how this special pre-factor comes in because $delta(a x) = delta(x)/|a|.$ Once you see this pre-factor, I claim you have solved the entire problem. At long last you find that in fact, $$ c rho' = gamma ~(c rho - beta J_x), J_x' = gamma~ (J_x - beta c rho), J_y' = J_y,$$ and therefore the single particle current density has transformed like a 4-vector under boosts, with transformation matrix $$begin{bmatrix}gamma&-gammabeta&0-gammabeta&gamma&0&0&1end{bmatrix}.$$ Combine this with the fact that $mathbf J$ transforms like an ordinary vector under rotations, and you have that arbitrary single-particle $(crho, mathbf J)$s are 4-vectors, and thus their superpositions are, too.

For the $k-$particle the charge density and the charge current density are respectively begin{align} rho_kleft(mathbf{x},tright) & boldsymbol{=} e_k delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right] tag{01a}label{01a} mathbf{j}_kleft(mathbf{x},tright) & boldsymbol{=} e_k delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right]dfrac{mathrm d mathbf{x}_kleft(tright)}{mathrm d t}boldsymbol{=} e_k delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right]mathbf{u}_kleft(tright) tag{01b}label{01b} end{align} where begin{equation} mathbf{u}_kleft(tright)boldsymbol{equiv}dfrac{mathrm d mathbf{x}_kleft(tright)}{mathrm d t} tag{02}label{02} end{equation} is the velocity 3-vector of the $k-$particle.

The charge current density 4-dimensional vector is begin{equation} mathbf{J}_kboldsymbol{=} left(vphantom{tfrac{a}{b}}crho_k,mathbf{j}_kright)boldsymbol{=} e_k delta^3left(vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kright)left(vphantom{tfrac{a}{b}}c,mathbf{u}_kright) tag{03}label{03} end{equation} hence begin{equation} boxed{:: mathbf{J}_kboldsymbol{=} rho_{0k}mathbf{U}_k vphantom{dfrac{a}{b}}::} tag{04}label{04} end{equation} where begin{equation} boxed{::mathbf{U}_k boldsymbol{=}gamma_kleft(vphantom{tfrac{a}{b}}c,mathbf{u}_kright)vphantom{dfrac{a}{b}}::} tag{05}label{05} end{equation} is the velocity 4-vector of the $k-$particle, begin{equation} boxed{::rho_{0k}left(mathbf{x},tright) boldsymbol{=}e_k dfrac{delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right]}{gamma_k}boldsymbol{=}e_kleft(1boldsymbol{-}leftVertdfrac{mathrm d mathbf{x}_k}{c,mathrm d t}rightVert^2right)^{tfrac12} delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right]::} tag{06}label{06} end{equation} is a scalar coefficient, function of the space-time coordinates $left(mathbf{x},tright)$ and begin{equation} boxed{::gamma_k boldsymbol{=}left(1boldsymbol{-}dfrac{u_k^2}{c^2}right)^{boldsymbol{-}tfrac12} boldsymbol{=}left(1boldsymbol{-}leftVertdfrac{mathrm d mathbf{x}_k}{c,mathrm d t}rightVert^2right)^{boldsymbol{-}tfrac12}::} tag{07}label{07} end{equation} is the $gamma-$factor relevant to the velocity of the $k-$particle.

Now, if we prove that the charge current density 4-dimensional vector $mathbf{J}_k$, equation eqref{04}, for one particle is a Lorentz 4-vector then the total charge current density 4-dimensional vector
begin{equation} mathbf{J}_{texttt{total}}boldsymbol{=}sumlimits_{k}mathbf{J}_k tag{08}label{08} end{equation} will be also a Lorentz 4-vector.

Note that the vector $mathbf{J}_k$ is a scalar multiple of a Lorentz 4-vector, the scalar being the function $rho_{0k}left(mathbf{x},tright)$ of equation eqref{06} and the Lorentz 4-vector being the velocity 4-vector $mathbf{U}_k$ of equation eqref{05}. But it could be proved that

Proposition 1 : If a 4-dimensional vector is the scalar multiple of a Lorentz 4-vector, the scalar being a Lorentz invariant scalar function, then this 4-dimensional vector is also a Lorentz 4-vector.

Above Proposition 1 could be proved on the ground of the following Proposition 2

Proposition 2 : Let $X_{mu}$ be the components of an arbitrary covariant tensor of the first order. If the inner product $A^{mu}X_{mu}$ is an invariant, then $A^{mu}$ are the components of a contravariant tensor of the first order.

For a proof of Proposition 2 see my answer here

Proving an object is a 4-vector given its inner product with a 4-vector is a scalar

In the following we'll prove that the function $rho_{0k}left(mathbf{x},tright)$ of equation eqref{06} is a Lorentz invariant scalar function. It is the proper charge density of the $k-$particle as seen from the rest frame of this same particle. Note that the point charge $e_k$, a by hypothesis Lorentz invariant scalar, is a factor of $rho_{0k}left(mathbf{x},tright)$ so the Lorentz invariance of the latter is restricted to the proof that begin{equation} boxed{::left(1boldsymbol{-}leftVertdfrac{mathrm d mathbf{x}_k}{c,mathrm d t}rightVert^2right)^{tfrac12} delta^3left[vphantom{tfrac{a}{b}}mathbf{x}boldsymbol{-}mathbf{x}_kleft(tright)right]boldsymbol{=}texttt{Lorentz invariant scalar}::} tag{09}label{09} end{equation} Note that equation eqref{09} is a mathematical property of the Lorentz transformation of the Dirac function independent of the theory of Classical Electrodynamics. This invariant was so strange to me when I found it 4.5 years ago, see my answer here

How is 4-current a 4-vector?