In the context of evaluating the propagation of a flattened Gaussian beam, the following integral appears:
\begin{equation} \int (\mathbf x^T \mathbf F \mathbf x)^n \exp \left [ - \mathbf x^T \mathbf G \mathbf x + \mathbf x^T \mathbf h \right] d\mathbf x \end{equation} where x and h are 2 by 1 matrices and the uppercase matrices are 2 by 2
A solution to this integral is implicitly used in this paper:
https://www.sciencedirect.com/science/article/pii/S0030402604701394
that is used to find the propagation integral for a flattened Gaussian beam, however the written solution seems to be wrong.
The only thing reported is that the following integral was used: \begin{equation} \int_0^\infty x^{2n} \exp(-a^2x^2)\cos(xy)dx = (-1)^n \pi^{1/2} 2^{-(2n+1)} a^{-(2n+1)} \exp \bigg (-\frac{y^2}{4a^2} \bigg ) H_{2n} \bigg (\frac{y}{2a} \bigg ) \end{equation}
The closest integral in the exchange I've been able to find is: Evaluation of Multidimensional Integral
Thank you,
Alex
Please note that if you have the expression for $J(\mathbf{G},\mathbf{h})=\int \exp \left [ - \mathbf x^T \mathbf G \mathbf x + \mathbf x^T \mathbf h \right] d\mathbf x$, then you can calculate any "Gaussian times polynomial" by judicially differentiating (possibly several times) by $h_i$.
Alternatively, consider the expression $$J(\mathbf{G}+\lambda \mathbf{F},\mathbf h)=\int \exp \left [ - \mathbf x^T (\mathbf G + \lambda \mathbf F) \mathbf x + \mathbf x^T \mathbf h \right] d\mathbf x$$ and differentiate it by $\lambda$ ($n$ times) and then take $\lambda \rightarrow 0$