randomVec = #*Normalize@RandomReal[{-1, 1}, 3] & (*[taken from here][1]*)
n = 2; (*particles count*)
u = Table[randomVec@1, {n}]
matQ = Table[
   1/(2.*n)*
    Sum[3*Part[u, i, a]*Part[u, i, b] - KroneckerDelta[a, b], {i, 1,
      n}], {a, 3}, {b, 3}];

SeedRandom[123];
n = 1000000;
u = RandomPoint[Sphere[{0, 0, 0}], n];

f = u [Function] NDSolve`FEM`MapThreadDot[
 ArrayReshape[u, Append[Dimensions[u], 1]],
 ArrayReshape[u, Insert[Dimensions[u], 1, 2]]
 ];

Q = 0.5 (ConstantArray[3./n, n].f[u] - N@IdentityMatrix[3]);
{λ, U} = Eigensystem[Q];
order = Ordering[λ];
λ = λ[[order]];
U = U[[order]];
p = λ[[-1]]

U[[-1]]