Diagonal matrices have obvious advantages in speeding up numerical computations, and Wolfgang Bangerth's answer is a good explanation of how to calculate a diagonal mass matrix, but it doesn't answer the OP's question "why does this work" in the sense of "why is it a good approximation to the physics you are modelling".
Conceptually, you can separate the response of an element into three parts: translational motion of a rigid body, rigid rotation about the element center of mass, and the deformation of the element.
The basic function of the element mass matrix is to represent the element KE as a quadratic form (i.e. $\frac 1 2 v^T M v$ where $v$ are the nodal velocities).
As the element size decreases, the contribution to the KE from rigid body rotation decreases faster than the contribution from translation, (for a solid element with a typical linear size of $a$, the mass is proportional to $a^3$ but the moments of inertia are proportional to $a^5$) and the contribution from the element deformation is negligible (at least for problems with small elastic strains).
Therefore, you only really need a "good" approximation to the rigid body parts of the motion, i.e. 6 DOFs, and in fact a good approximation to only the KE from rigid body translation, i.e. 3 DOFs, will converge as the element size is reduced.
The diagonal terms of the element matrix contain more than enough independent parameters to represent those 3 or 6 KE terms with sufficient accuracy. In fact for higher order elements, you can use mass diagonal mass matrices where the diagonal terms for the mid-side nodes are zero.
Note that this is a completely different situation from the element potential energy, where the contributions from rigid body translation and rotation are zero, and the only thing that matters is representing the strain energy corresponding to the element deformation. A diagonal stiffness matrix would therefore not be a feasible idea!