One way is to consider that Lorentz transformations apply more fundamentally to momentum/ energy, than to space/time. So, with a boost transformation along the $z$ axis, we will have :
$\begin {pmatrix} p'_z \\E' \end{pmatrix} = \gamma(v)\begin {pmatrix} 1 & -\frac{v}{c} \\-\frac{v}{c} &1\end{pmatrix}\begin {pmatrix} p_z \\E \end{pmatrix}$
It is not difficult to see that, if, $|\large \frac{p_zc}{E}|=1$, then $|\large \frac{p'_zc}{E'}|=1$
Now, by dimensional analysis, we have : $\frac{\vec Pc}{E} = \frac{\vec V}{c}$, where $\vec V$ has the dimension of a velocity. The most natural possibility is that $\vec V$ is the velocity of the particle (for instance, with $v=V_z$, you have $V'_z=0$). So, we see, that a particuleparticle with speed $|\vec V|=c$ in a galileanGalilean frame, has also $|\vec V'|=c$ in another galileanGalilean frame. We could also check that the quantity $E^2-\vec p^2c^2$ is conserved by Lorentz transformations, and call this quantity $m^2c^4$, where $m$ is the mass. So, particles who have $|\frac{\vec Pc}{E}| = |\frac{\vec V}{c}|=1$$\left|\frac{\vec Pc}{E}\right| = \left|\frac{\vec V}{c}\right|=1$, are massless particles.