Moments of the trace of orthogonal matrices Let $O_n$ be the (real) orthogonal group of $n$ by $n$ matrices.
I am interested in the following sequence which showed up in a calculation I was doing
$$a_k = \int_{O_n} (\text{Tr } X)^k dX$$
where the integral is taken with respect to the normalized Haar measure on $O_n.$  I noticed some work of Dianconis et al. covers a more general set of moments, but essentially the formulas I saw only seem to hold for small $k.$ For the case of $O_2$ the sequence looks like $a_{2k} = \frac{1}{2}\binom{2k}{k},$ which can be calculated by hand since $O_2$ has some special structure.
I was hoping  there would be some general known formula (since, if I recall correctly, this calculates the multiplicity the trivial representation in a tensor power the standard representation of $O_n$), or that it might be in the OEIS, but I don't know exactly what to look for, and don't know enough small terms to search.
 A: This answer is a follow-up to the other answers (particularly to Carlo Beenakker's answer from September 4).
First off, Carlo's conjecture is indeed true. That is:

Theorem. If $n \geq k$ then
  $$a_{2k}=\pi^{-1/2}2^{k}\Gamma(k+1/2).\quad [*]$$

Idea of proof: Recall that the Brauer algebra describes the set of invariants of $2k$-fold tensor copies of orthogonal matrices, which the following operator projects onto:
$$
\int_{O_n}X^{\otimes 2k} dX. \quad [\dagger]
$$
If we let $\{e_i\}$ be the standard basis of $\mathbb{R}^n$, then after you unwrap everything, you find that $[\dagger]$ is the orthogonal projection onto the span of the $\pi^{-1/2}2^{k}\Gamma(k+1/2)$ different vectors of the form
$$
W_\pi\sum_{i_1,i_2,\ldots,i_k=1}^n e_{i_1} \otimes e_{i_1} \otimes e_{i_2} \otimes e_{i_2} \otimes \cdots \otimes e_{i_k} \otimes e_{i_k}, \quad [**]
$$
where $W_\pi$ is a unitary operator that permutes the $2k$ different tensor factors according to an arbitrary permutation $\pi$ (this is essentially Theorem 3.1 in [1]). The reason that there are exactly $\pi^{-1/2}2^{k}\Gamma(k+1/2)$ vectors of the form $[**]$ is that this is the number of perfect matchings of $2k$ objects (which in turn is equal to the dimension of the Brauer algebra, as expected).
So far, this shows that $a_{2k} \leq \pi^{-1/2}2^{k}\Gamma(k+1/2)$ always (regardless of $n$ and $k$). To show that the other inequality holds when $n \geq k$, it suffices to show that the vectors $[**]$ are linearly independent. I will just refer to Theorem 3.4 of [1] for this claim, but it's not difficult to prove directly using standard linear algebra tools.

This still leaves the question of what $a_{2k}$ equals when $n < k$. Here are some random observations and conjectures based on the numerical evidence found by Carlo and that I've dug up elsewhere:

Conjecture. If $n = k-1$ then $a_{2k} = \pi^{-1/2}2^{k}\Gamma(k+1/2) - \frac{(2k)!}{k!(k+1)!}$. This is the number of perfect matchings of $2k$ objects that have at least one "crossing" (in the same sense that the Catalan numbers give the number of perfect matchings of $2k$ objects with no "crossings").
Conjecture. If $n = k-2$ then $a_{2k}$ seems to be the number of perfect matchings of $2k$ objects with at least three crossings (there is a known formula for this quantity, but it's messy so I won't write it here).

The other answers gave formulas for when $n = 2$ or $n = 3$. It turns out that there is also already a known formula for the $n = 5$ case (see OEIS A095922: the even-indexed terms are the values of $a_{2k}$ that we want). That OEIS entry cites the book [2], but I haven't been able to track down a copy of it. Numerical evidence also suggests the following:

Conjecture. If $n = 4$ then
  $$ a_{2k} = \frac{1}{2}\frac{(2k)!}{k!(k+1)!}\left(\frac{(2k)!}{k!(k+1)!}+1\right).$$


References:


*

*G. Lehrer and R. Zhang, The second fundamental theorem of invariant theory for the orthogonal group, Ann. of Math. 176:2031-2054 (2012).

*A. Mihailovs, A Combinatorial Approach to Representations of Lie Groups and Algebras, Springer-Verlag New York (2004).

A: Below I'll mention a classical result (due to Alan T. James (1964), * ). For the moments to be nonzero, $k$ must be even, so following Carlo, I'll also write $a_{2k}$ to denote the following quantity:
\begin{equation*}
  a_{2k}(X) := E_X[ \mbox{tr}(AX)^{2k}] = \sum_{\tau \in P_k} \frac{2^k \cdot k!\, d_\tau}{(2k)!C_\tau(I_n)}C_\tau(AA^T)
\end{equation*}
where $n\times n$ orthogonal $X$ is uniformly distributed, $C_\tau(\cdot)$ denotes the (normalized) Zonal symmetric polynomial indexed by partition $\tau$ ($P_k$ is the set of partitions of integer $k$), and it is known that
\begin{equation*}
 d_\tau = \frac{\chi_{[2\tau]}(1)2^kk!}{(2k)!} = \frac{2^k k!\prod_{i< j}^{\ell(\tau)} (2\tau_i -2\tau_j -i + j)}
 {\prod_{i=1}^{\ell(\tau)} (2\tau_i + \ell(\tau) - i)!},
\end{equation*}
where $\chi_{[2\tau]}(1)$ is the dimension of the representation of $[2\tau]=[2\tau_1,\ldots,2\tau_{\ell(\tau)}]$ of the symmetric group on $2k$ symbols, and $\ell(\tau)$ is the length (i.e., number of parts) of the partition $\tau$.
This quantity was obtained by James (1964) using representation theory (as is now evident from Qiaochu's more general answer).
As noted, in the OP, $A=I_n$, so actually the desired moment is (for $k\le n$)
\begin{equation*}
a_{2k}(I_n) = \frac{2^k k!}{(2k)!}\sum_{\tau \in P_k} d_\tau,
\end{equation*}
which should hopefully simplify further, as conjectured by Carlo.

The reference (that on page 40 includes the abovecited result) from which I took the above material: Zonal Polynomials by Akimichi Takemura (see e.g., google books link) --- though note that this book uses $n$ where I've written $k$ (to be consistent with the OP).
A: I am writing more than a year after the question was posted, only to spell out some more details regarding the calculation implicit in the solution presented by Suvrit, and to clarify the dependence on $n$ in that solution.
Write $[{\rm Tr}(AX)]^{2k}=p_{1^{2k}}(AX)$, using power sum symmetric functions. Next, expand into Schur functions, $p_{1^{2k}}(AX)=\sum_{\lambda\vdash 2k}d_\lambda s_\lambda(AX)$. Notice that in this expansion we must have $\ell(\lambda)\le n$, otherwise the Schur function is identically zero.
Now use the definition of zonal polynomials (average of irreducible characters of $U(n)$ over the subgroup $O(n)$),
$$ \int_{O(n)}s_\lambda(AX)dX=\frac{Z_\tau(AA^T)}{Z_\tau(1^n)},$$ which holds if $\lambda=2\tau$ (the integral vanishes otherwise) to get
$$ \int_{O(N)}[{\rm Tr}(X)]^{2k}dX=\sum_{\tau\vdash k}d_{2\tau}, $$ with the condition $\ell(\tau)\le n$ (which is the only dependence of the answer on $n$).
A: Here is a general comment. Let $G$ be a compact group and let $V$ be a (finite-dimensional, continuous, complex) representation of $G$. This data determines a locally finite directed graph, the representation graph $\Gamma$ of $G$ and $V$, as follows. The vertices of $\Gamma$ are the irreducible representations of $G$, and the number of edges from an irreducible representation $U_1$ to an irreducible representation $U_2$ is the multiplicity of $U_2$ in $U_1 \otimes V$.
Let $\chi_V$ denote the character of $G$.

Claim: The moment
$$\int_G \chi_V(g)^k \, dg$$
is the number of closed walks of length $k$ from the trivial representation to itself in $\Gamma$.

The idea of the proof is to look at the relationship between the decomposition of $V^{\otimes k}$ into irreducibles and the decomposition of $V^{\otimes (k+1)}$ into irreducibles. The application is that if $G$ is sufficiently nice then we can try to figure out what $\Gamma$ looks like and work from there.
Let's do $G = \text{O}(3)$ ($V$ the complexification of the standard representation) this way to demonstrate how the method works. First observe that we can immediately reduce to $G = \text{SO}(3)$, since the integral over $\text{O}(3)$ splits up into an integral over $\text{SO}(3)$ and over the other coset, and since the latter is $-1$ times the former (here we are using the fact that $3$ is odd) we get that
$$\int_{\text{O}(3)} \chi_V(g)^{2k+1} \, dg = 0$$
and
$$\int_{\text{O}(3)} \chi_V(g)^{2k} \, dg = \int_{\text{SO}(3)} \chi_V(g)^{2k} \, dg.$$
Recall that the irreducible representations of $\text{SO}(3)$ are precisely the odd-dimensional irreducible representations $V_1, V_3, V_5 \dots$ of its Lie algebra $\mathfrak{so}(3) \cong \mathfrak{su}(2)$, with $V = V_3$.

Claim: For $n \ge 3$, we have $V_n \otimes V_3 \cong V_{n-2} \oplus V_n \oplus V_{n+2}$.

This should be familiar from the representation theory of $\mathfrak{su}(2)$, and it gets us a very explicit description of the representation graph in this case: the vertices are labeled by the odd positive integers $1, 3, 5 \dots$ and for every vertex there are three edges, one from each vertex to its left neighbor, to itself, and to its right neighbor. (The special case $V_1$ is easy to handle since $V_1 \otimes V_3 \cong V_3$.)
As a corollary, closed walks from $1$ to itself of length $k$ are very close to being counted by the Motzkin numbers $M_k$ (A001006), which would be the correct answer if the representation graph had an additional edge from $1$ to itself. Some computations reveal that the actual sequence we get is A005043, which I've never seen before but which OEIS says are called the Riordan numbers. The sequence begins
$$1, 0, 1, 1, 3, 6, 15, 36, 91 \dots$$
so to get the original sequence of moments we wanted, we start with this sequence and replace every odd term with $0$, giving
$$1, 0, 1, 0, 3, 0, 15, 0, 91 \dots$$
This sequence appears in the OEIS with the zeroes removed as A099251.
For general compact connected Lie groups the representation graph should look a bit like the intersection of a Weyl chamber with the weight lattice, or something like that. I think the combinatorics gets much hairier as soon as the rank is higher than $1$ though.
A: 
Pastur and Vasilchuk have extended the result of Diaconis and Evans for $a_{2k}$ from $2k\leq n/2$ to $2k\leq n-1$:
$$a_{2k}=\pi^{-1/2}2^{k}\Gamma(k+1/2)\;\;\text{for}\;\;2k\leq n-1\quad\quad[*]$$

As suggested by Liviu Nicolaescu, for small $n$ you might directly integrate over the probability distribution of the eigenvalues in $O_n$, which you can find here. (This is the socalled "circular real ensemble" of random matrices.) 
$\bullet$ For $n=2$ one reproduces the result 
$$a_{2k}=\frac{1}{\pi}\int_0^{\pi/2}dx\,(2\cos x)^{2k}=\tfrac{1}{2}(2k)!(k!)^{-2}$$
you quoted above, which agrees with $[*]$ for $k=1$ and $k=2$. (For $k=3$ it gives $10$ instead of $15$.)
$\bullet$ For $n=3$ one obtains
$$a_{2k}=\frac{1}{2\pi}\int_0^{2\pi}dx\,(1-\cos x)(1+2\cos x)^{2k}$$
which evaluates to
$$a_2=1,\;\; a_4=3,\;\; a_6=15,\;\; a_8=91,\;\;a_{10}=603,\;\;a_{12}=4213$$
in agreement with the series of Qiaochu Yuan's answer. The formula $[*]$ for $a_{2k}$ holds for $k\leq 3$ (for $k=4$ it gives $105$ instead of $91$).
$\bullet$ For $n=4$ one finds
$$a_{2k}=\frac{1}{8\pi^2}\int_0^{2\pi}dx\int_0^{2\pi}dy\,(\cos x-\cos y)^2(2\cos x+2\cos y)^{2k}$$
$$\quad\quad+\frac{1}{2\pi}\int_0^{2\pi}dz\,(1-\cos^2 z)(2\cos z)^{2k}$$
which evaluates to
$$a_2=1,\;\;a_4=3,\;\;a_6=15,\;\;a_8=105,\;\;a_{10}=903,\;\;a_{12}=8778.$$
Now the formula $[*]$ for $a_{2k}$ holds for $k\leq 4$ (for $k=5$ it gives $945$ instead of $903$).
$\bullet$ Continuing with $n=5$, we have
$$a_{2k}=\frac{1}{2\pi^2}\int_0^{2\pi}dx\int_0^{2\pi}dy\,(1-\cos x)(1-\cos y)(\cos x-\cos y)^2(1+2\cos x+2\cos y)^{2k}$$
which evaluates to
$$a_2=1,\;\;a_4=3,\;\;a_6=15,\;\;a_8=105,\;\;a_{10}=945,\;\;a_{12}=10263.$$
The formula $[*]$ for $a_{2k}$ holds for $k\leq 5$ (for $k=6$ it gives $10395$ instead of $10263$).

I am now tempted to $$\text{conjecture}\quad\quad a_{2k}=\pi^{-1/2}2^{k}\Gamma(k+1/2)=(2k-1)!!\;\;\text{for}\;\;k\leq n\quad\quad[**]$$
  thus doubling the range of validity of $[\ast]$. Is it true?

