Solution (source code)

= Solution

Let a change of orthonormal coordinates be represented by a <rotation matrix> $R$. Since both <angular momentum> and <angular velocity> are vectors,
$$
L'_i=R_{ip}L_p,
\qquad
\omega'_j=R_{jq}\omega_q.
$$
Using $L_i=I_{ij}\omega_j$ in both frames gives
$$
I'_{ij}R_{jq}\omega_q=R_{ip}I_{pq}\omega_q
$$
for every vector $\boldsymbol\omega$. Multiplication by $R^{-1}=R^T$ yields
$$
\boxed{I'_{ij}=R_{ip}R_{jq}I_{pq}},
$$
which is exactly the transformation law for a rank-two <tensor>. Thus $I_{ij}$ is a rank-two tensor.

Put $\mathbf r=\mathbf x-\mathbf a$. The <vector triple product> gives
$$
\mathbf r\times(\boldsymbol\omega\times\mathbf r)
=r^2\boldsymbol\omega
-\mathbf r(\mathbf r\cdot\boldsymbol\omega).
$$
Therefore
$$
\boxed{
I_{ij}(\mathbf a)
=\rho\int_B
\left[(x_k-a_k)(x_k-a_k)\delta_{ij}
-(x_i-a_i)(x_j-a_j)\right]dV}.
$$
At the centre of mass,
$$
\boxed{
I_{ij}(\mathbf0)
=\rho\int_B(x_kx_k\delta_{ij}-x_ix_j)\,dV}.
$$
Expanding the first formula, the terms linear in $x_i$ vanish because the centre of mass is the origin, while $\rho\int_BdV=M$. Hence the <parallel axis theorem> in tensor form is
$$
\boxed{
I_{ij}(\mathbf a)
=I_{ij}(\mathbf0)
+M(a_ka_k\delta_{ij}-a_ia_j)}.
$$