Evaluation of Gaussian multivariable integral 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

*

*$\mathbf x$ and $\mathbf h$ are $2\times 1$ matrices and

*the uppercase matrices $\mathbf F$ and $\mathbf G$ are $2 \times 2$.

A solution to this integral is implicitly used in this paper, where it 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 the one in this Q&A.
Any hint on how to proceed to solve the first expression will be greatly appreciated.
Thank you,
Alex
(Question also posted here)
 A: Multiplying by $\epsilon^{n} /n! $ and summing over $n$, we obtain the generating function
\begin{eqnarray}
M &=& \sum_{n=0}^{\infty } \frac{\epsilon^{n} }{n! }
\int d^2 x\ (\mathbf x^T \mathbf F \mathbf x)^n \exp \left [ - \mathbf x^T  \mathbf G \mathbf x + \mathbf x^T \mathbf h \right] \\
&=& \int d^2 x\ \exp \left [ - \mathbf x^T (\mathbf G - \epsilon \mathbf F) \mathbf x + \mathbf x^T \mathbf h \right]
\end{eqnarray}
Assuming $\mathbf G - \epsilon \mathbf F$ is symmetric, we can diagonalize
$$
O (\mathbf G - \epsilon \mathbf F) O^T = D = \mbox{diag} (d_1 ,d_2 )
$$
with orthogonal $O$ (otherwise, one has to be a bit more careful with the transformation). Substituting $\mathbf y=O\mathbf x$ and $\mathbf k=O\mathbf h$, we have
\begin{eqnarray}
M &=& \int d^2 y\ \exp \left [ - \mathbf y^T D \mathbf y + \mathbf y^T \mathbf k \right] \\
&=& \frac{\pi }{\sqrt{d_1 d_2 } }\exp \left[ \frac{k_1^2 }{4d_1 } +  \frac{k_2^2 }{4d_2 } \right]
\end{eqnarray}
where $d_i $, $k_i $ can be given explicitly in terms of the elements of $\mathbf F $, $\mathbf G$, and $\mathbf h$, as well as $\epsilon $, and the original problem is reduced to obtaining the coefficient of $\epsilon^{n} $ in the expansion of $M$. This is now a purely mechanical calculation, but it is not clear to me whether it reduces to anything reasonably succint without further assumptions on $\mathbf F $, $\mathbf G$, and $\mathbf h$.
