Let $X\in\mathbb{R}^{d}$ have independent, mean zero entries that are strictly subgaussian with parameter $\sigma^{2}$ and let $A=ab^{\top}$ for some $a,b\in\mathbb{R}^{d}$. I am interested in computing \begin{equation} \mathbb{E}[(X^{\top}AX)^{q}], \end{equation} for $q\in\mathbb{N}$. For example, when $q=1$, it is easy to see that this is just $\sigma^{2}\text{Tr}(A)=\sigma^{2}(a^{\top}b)$.
This work [1] for example computes the first four moments when $X$ is Gaussian and $A$ is symmetric [equation 2.32]. Overall, it is a combination of polynomials in the trace of $A$ (or its higher powers).
Any ideas on how this would extend to subgaussian $X$ and nonsymmetric rank-one $A=ab^{\top}$?
[1] Issie Scarowsky, Quadratic forms in normal variables. 1973