I'm working on a problem that arises in some branches of physics. Let assume $$ \boldsymbol{I}_{k,D} = \frac{1}{\left(\sqrt{2\pi}\right)^{D}}\int_{\mathbb{R}^{D}} \boldsymbol{x}^{k}\text{exp}\left(-\left\|\boldsymbol{x}\right\|^2/2\right)\text{d}\boldsymbol{x} $$ where $\left\| \right\|$ is the euclidian norm and $\boldsymbol{x}$ is a tensor of order $k$ with $\left(\boldsymbol{x}^{k}\right)_{i_1, i_2, \dots, i_{k}} = \prod_{p=1}^{k}x_p$.
I want to know if it is possible to find a general formula to compute $I_{2k,D}$. For example, $$ I_{0,1} = 1, \ \boldsymbol{I}_{2, D} = \boldsymbol{I}_2, \ \boldsymbol{I}_{2k+1, D}=\boldsymbol{0} \text{ and } I_{2k, 1} = \frac{k!}{2^{k/2}\left(k/2\right)!} $$ where $\boldsymbol{I}_2$ is the identity matrix of size $2 \times 2$. What I've found so far is that for $\boldsymbol{x} = \left(x_i\right)_{i=1}^{D}$, we have $$\begin{aligned} \left(\boldsymbol{I}_{k,D}\right)_{i_1, i_2, \dots, i_{2k}} &= \frac{1}{\left(\sqrt{2\pi}\right)^{D}}\int_{\mathbb{R}^{D}} \prod_{k=1}^{D}x_i^{p_i}\text{exp}\left(-\left\|\boldsymbol{x}\right\|^2/2\right)\text{d}\boldsymbol{x}\\ &= \prod_{k=1}^{D}I_{p_i, 1} \end{aligned} $$ which happens because of the exponential properties. The integers $p_i$ depend on the indexes but $p_1 + p_2 + \dots + p_{D} = k$.
For example, if I choose $D = 2$ $$ \left(\boldsymbol{I}_{k,2}\right)_{i_1, i_2, \dots, i_{2k}} = I_{p_1, 1} I_{p_2, 1} $$ with $p_1+p_2 = k$ with $\left(p_1, p_2\right) \in \left\{0, 1, 2\right\}^2$. If $p_i$ is odd, then $I_{p_i, 1} = 0$ which means that the corresponding coefficient of $\boldsymbol{I}_{k,2}$ is null. The only non-zero component of $\boldsymbol{I}_{k,2}$ happens when $p_1$ and $p_2$ are even so here $\left(2,0\right)$ and $\left(0,2\right)$ which leads to the identity matrix for $k=2$. This can be generalized, the non-zero component for $k=4$ would be at indexes $\left(4, 0\right), \left(0,4\right), \left(2,2\right)$ but I'm struggling writing a general formula for $\boldsymbol{I}_{2k,D}$.
Any help would be greatly appreciated