For this problem, can anyone explain to me why when $x_1$ axis is aligned with the diagonal of the cube, the resulting inertia tensor will become diagonal?
How to interpret this result geometrically? I mean, is there some tricky way for some geometric shape that I don't have to compute the eigenvalue, eigenvectors to find the principal axises, but by observation, I can make a guess first.
Cross terms appear when the coordinate axes do not pass through the center of mass. That is if you start with a diagonal inertia matrix at the center of mass, when applying the parallel axis theorem cross terms will appear.
In vector form the parallel axis theorem is
$$ {\bf I} = {\bf I}_{cm} - m [{\bf r}\times] [{\bf r}\times] $$
where $[{\bf r}\times] = \begin{pmatrix}x\\y\\z\end{pmatrix} \times = \begin{bmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{bmatrix}$ is the cross product matrix operator.
So if we start with a diagonal inertia at the center of mass, when moved to a different point $(x,y,z)$ the inertia matrix is
$$ {\bf I} = \begin{bmatrix} I_x+m(y^2+z^2) & - m x y & -m x z \\ - m x y & I_y + m (x^2+z^2) & - m y z \\ - m x z & - m y z & I_z + m (x^2+y^2) \end{bmatrix} $$
So if any two of $x$ $y$ or $z$ are zero the result is still a diagonal matrix. This happens when one of the coordinate axis passes through the center of mass. In your case if the $x$ axis goes from the corner towards the center of mass (across diagonal) then $y=0$ and $z=0$ and the criteria is met.
So the question boils down to under which conditions the inertial matrix is diagonal at the center of mass. The answer to this has to do with symmetries. For example, when a particle at a positive $z$ coordinate is matched with an particle at a negative $z$ coordinate then when the cross terms get added up the result is zero. This can be seen from the structure of the inertia matrix (as shown above) if you think of this as a parallel axis theorem for a lot of small little point masses added together.
