Moment of a function of Gaussian random variables: $\mathbb{E}[(a_{i}^{\top}AA^{\top}a_{j})^{q}]$

Let $$A$$ be an $$m\times k$$ matrix with iid $$\mathcal{N}(0,1)$$ entries and $$a_{i}$$ and $$a_{j}$$ be its $$i$$th and $$j$$th columns, respectively. I would like to compute the following quantity: $$$$\mathbb{E}[(a_{i}^{\top}AA^{\top}a_{j})^{q}],$$$$ for $$q\in\mathbb{N}$$.

One can expand $$a_{i}^{\top}AA^{\top}a_{j}$$ to get $$$$\|a_{i}\|_{2}^{2}a_{i}^{\top}a_{j}+\|a_{j}\|_{2}^{2}a_{i}^{\top}a_{j}+\sum_{t\in[k]\backslash\{i\}\cup\{j\}}(a_{i}^{\top}a_{t})(a_{j}^{\top}a_{t})$$$$ and explicitly compute the expectation for $$q=1, 2$$ etc. But how can I generalize for larger $$q$$?

[Posting this here since couldn't get answer on Math Stack Exchange and Cross Validated]

• I don't think it is reasonable to expect an exact answer in closed form; consider the simpler problem where instead of asking for moments of the matrix element $B_{ij}=a_i^T AA^T a_j$ you would ask for moments of ${\rm tr}\,B={\rm tr}\,W^2$ with Wishart matrix $W=AA^T$; for large $m,k$ you could use Marcenko-Pastur and make analytical progress, but exact results for any $m,k$ seem out of reach. – Carlo Beenakker Dec 6 '19 at 20:17