Can anybody suggest how I can compute the second moment (or the whole moment generating function) of the cosine of two gaussian random vectors $x,y$, each distributed as $\mathcal N (0,\Sigma)$, independent of each other? IE, moment for the following random variable
$$\frac{\langle x, y\rangle}{\|x\|\|y\|}$$
The closest question is Moment generating function of the inner product of two gaussian random vectors which derives MGF for the inner product. There's also this answer from mathoverflow which links this question to distribution of eigenvalues of sample covariance matrices, but I don't immediately see how to use those to compute the second moment.
I suspect that second moment scales in proportion to half-norm of eigenvalues of $\Sigma$ since I get this result through algebraic manipulation for 2 dimensions, and also for 3 dimensions from guess-and-check. For eigenvalues $a,b,c$ adding up to 1, second moment is:
$$(\sqrt{a}+\sqrt{b}+\sqrt{c})^{-2}$$
Using the following for numerical check
val1[a_, b_, c_] := (a + b + c)/(Sqrt[a] + Sqrt[b] + Sqrt[c])^2
val2[a_, b_, c_] := Block[{},
x := {x1, x2, x3};
y := {y1, y2, y3};
normal := MultinormalDistribution[{0, 0, 0}, ( {
{a, 0, 0},
{0, b, 0},
{0, 0, c}
} )];
vars := {x \[Distributed] normal, y \[Distributed] normal};
NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]]
val1[1.5,2.5,3.5] - val2[1.5,2.5,3.5]
Checking the formula for 4 variables (within numerical bounds):
val1[a_, b_, c_,
d_] := (a + b + c + d)/(Sqrt[a] + Sqrt[b] + Sqrt[c] + Sqrt[d])^2
val2[a_, b_, c_, d_] := Block[{},
x := {x1, x2, x3, x4};
y := {y1, y2, y3, y4};
normal :=
MultinormalDistribution[{0, 0, 0,
0}, {{a, 0, 0, 0}, {0, b, 0, 0}, {0, 0, c, 0}, {0, 0, 0, d}}];
vars := {x \[Distributed] normal, y \[Distributed] normal};
NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]]
val1[0.5, 1.5, 2.5, 3.5] - val2[0.5, 1.5, 2.5, 3.5]