# To prove the Lorentz invariance of density distribution functions for massless particles in phase space

One defines the density distribution function of a collection of $N$ particles in phase space as follows, $$f(\vec{x},\vec{p},t)=\sum_{i=1}^N\delta^{(3)}(\vec{x}-\vec{x}_i)\delta^{(3)}(\vec{p}-\vec{p}_i).$$ The claim is that the above distribution function is Lorentz invariant. To prove the claim one defines the following function, $$F(\vec{x},\vec{p},t)=\delta(p^2-m^2)\Theta(p^0)f(\vec{x},\vec{p},t).$$ If one is able to prove that $F$ is Lorentz invariant, then $f$ will also be Lorentz invariant. By manipulations one arrives at, $$F(\vec{x},\vec{p},t)=\frac{1}{2E_{\vec{p}}}\sum_{i=1}^N\delta^{(3)}(\vec{x}-\vec{x}_i)\delta^{(4)}(p-p_i).$$ As the delta function restricting the particle momenta on the shell has already converted the 3-dimensional delta function of momenta into a 4-dimensional delta functions, one only worries about the 3-delta function of position. To tackle this problem, integration over coordinate time is introduced, $$F(\vec{x},\vec{p},t)=\frac{1}{2E_{\vec{p}}}\sum_{i=1}^N\int dt_i\delta(t-t_i)\delta^{(3)}(\vec{x}-\vec{x}_i(t_i))\delta^{(4)}(p-p_i(t_i)).$$ Using the relation between proper time and coordinate time $d\tau_i=\frac{m}{E_{\vec{p}}}dt_i$, and suppressing the the index in proper time one arrives at the following equation, $$F(\vec{x},\vec{p},t)=\frac{1}{2m}\sum_{i=1}^N\int d\tau\delta^{(4)}(x-x_i(\tau))\delta^{(4)}(p-p_i(\tau)).$$ The above expression suggests $F$ is Lorentz invariant and therefore, so is $f$. But the above treatment only works for massive particles. My question is how to generalize this method to massless particles for which the mass and proper time both are inevitably zero?