Solution (source code)

= Solution

In orthonormal Cartesian coordinates, a second-rank <tensor> has components transforming under a change of axes $Q$ as $M'_{ij}=Q_{ik}Q_{j\ell}M_{k\ell}$, where $QQ^T=I$. Its contraction is invariant:
$$
M'_{ii}=Q_{ik}Q_{i\ell}M_{k\ell}=\delta_{k\ell}M_{k\ell}=M_{kk}.
$$
Thus \b[$M_{ii}$ is a scalar]. This is <tensor contraction> with the <Euclidean metric>. In general nonorthonormal coordinates the invariant expression is $g^{ij}M_{ij}$, not an unweighted sum of diagonal covariant components.

For the planar plate, every mass element has $z=0$. The <inertia tensor> is consequently
$$
M=\begin{pmatrix}
\int_Dy^2\rho\,dS&-\int_Dxy\rho\,dS&0\\
-\int_Dxy\rho\,dS&\int_Dx^2\rho\,dS&0\\
0&0&\int_D(x^2+y^2)\rho\,dS
\end{pmatrix}.
$$
It follows directly that $\mathbf e_z$ is an <eigenvector> and
$$
\boxed{M_\perp=\int_D(x^2+y^2)\rho\,dS.}
$$
The real <symmetric matrix> has three real <eigenvalues>. Its <trace> is $2\int_D(x^2+y^2)\rho\,dS=2M_\perp$, so $M_1+M_2+M_\perp=2M_\perp$. Hence the <perpendicular axis theorem> gives \b[$M_\perp=M_1+M_2$].

For a disc of constant areal <density> $\rho_0$, symmetry gives $\int_Dxy\,dS=0$, and polar integration gives $\int_Dx^2\,dS=\int_Dy^2\,dS=\pi a^4/4$. Therefore
$$
\boxed{M=\frac{\pi\rho_0a^4}{4}\operatorname{diag}(1,1,2)
=\frac{Ma^2}{4}\operatorname{diag}(1,1,2),\qquad M=\pi\rho_0a^2,}
$$
where the scalar $M$ in the last expression denotes the total mass rather than the <tensor>.