Solution (source code)

= Solution

Here $I_{ij}$ denotes the unreduced <second mass moment tensor>, distinct from its trace-free <mass quadrupole moment>. In coordinates referred to the chosen origin,
$$
\boxed{I_{ij}(t)=\int_{\mathcal V}T^{00}(t,\mathbf y)y_i y_j\,d^3y.}
$$
This expression uses $c=1$. In SI units the leading nonrelativistic mass density is $T^{00}/c^2$. For slowly moving <point masses>, $T^{00}\simeq\sum_A m_{(A)}\delta^{(3)}(\mathbf y-\mathbf y_{(A)}(t))$, so
$$
\boxed{I_{ij}(t)=\sum_{A=1}^N m_{(A)}y_{(A)i}(t)y_{(A)j}(t),\qquad
Q_{ij}=I_{ij}-\frac13\delta_{ij}I_{kk}.}
$$
Kinetic corrections to the energy density are higher order in the velocity. The <second mass moment tensor> is not the mechanical <moment of inertia> tensor, which instead has components $\delta_{ij}I_{kk}-I_{ij}$.