I'd like to compute the following (Riemann) integral: $$\int_{\textbf{R}^n} \langle x,a \rangle ^2 e^{-\frac{1}{2} \|x\| ^2}$$ for some fixed $a\in \textbf{R}^n$. To get started I computed the case in $\textbf{R}^2$, letting $a=(p,q)$ and writing $x=(x,y)$ (pardon the abuse of notation.) We can write $\|(x,y)\|^2=x^2+y^2.$ So that \begin{align*} \int_{\textbf{R}^2} \langle x,a \rangle ^2 e^{-\frac{1}{2} \|x\| ^2}&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}(p^2x^2+2pqxy+q^2y^2)e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy\\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}p^2x^2e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy + \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}2pqxy e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy\\&+ \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}q^2y^2e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy. \end{align*} Applying Fubini's theorem and using the well known Gaussian integral, the first integral summand is $$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}p^2x^2e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy=p^2\int_{\infty}^{\infty}x^2e^{-\frac{1}{2}x^2}\,dx \int_{\infty}^{\infty}e^{-\frac{1}{2}y^2}\,dy=2\pi p^2.$$ By symmetry, $$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}q^2y^2e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy=2\pi q^2$$ If we factor out $2pq$ from the middle summand integral, and apply Fubini's theorem, we get $$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}2pqxy e^{-\frac{1}{2}x^2}e^{-\frac{1}{2}y^2}\,dxdy=2pq\int_{-\infty}^{\infty}xe^{-\frac{1}{2}x^2}\,dx\int_{-\infty}^{\infty}ye^{-\frac{1}{2}y^2}\,dy.$$ Both integrals are of odd functions taken over a domain symmetric about 0, so they both vanish. Then we get $$\int_{\textbf{R}^2} \langle x,a \rangle ^2 e^{-\frac{1}{2} \|x\| ^2}=2\pi(p^2+q^2).$$
This leads me to the following guess in the $\textbf{R}^n$ case: if $a=(a_1,a_2,\ldots,a_n)$, then $$\int_{\textbf{R}^n} \langle x,a \rangle ^2 e^{-\frac{1}{2} \|x\| ^2}=(2\pi)^{\frac{n}{2}}(a_1^2+\ldots+a_n^2).$$ The justification for this guess would similar to the result I got from the above computation: if we write $x=(x_1,\ldots,x_n),$ then we can split the integral by linearity according to the multinomial expansion of $\langle x,a \rangle ^2.$ In this expansion, all integrals of functions containing products of mixed $x_i$'s will vanish, since the maximal power of any term of mixed $x_i$'s will be 1 by the definition of squares of multinomials. This will leave only integrals of squares of one such $x_i$. Splitting each such integral by Fubini's theorem and evaluating will result in an $n$-fold product of $\sqrt{2\pi},$ multiplied by the corresponding square $a_i^2.$
Is this guess justified correctly? If not, any hints in the right direction would be appreciated. I have a suspicion this is not the slickest way to justify this.
First write $a = \alpha \mathbf{e}$ where $\alpha^2 = \|a\|^2= a_1^2 + \dots + a_n^2$ and $\mathbf{e}$ is a unit vector. The integral become $$ \int_{\mathbb{R}^n} \langle x,a \rangle ^2 e^{-\frac{1}{2} \|x\| ^2} = \alpha^2\int_{\mathbb{R}^n} \langle x,\mathbf{e} \rangle ^2 e^{-\frac{1}{2} \|x\| ^2} \, . $$ Next apply a rotation that maps $\mathbf{e}$ to $\mathbf{e_1}$, the first unit vector, and leaves the space that is perpendicular to $\mathbf{e}$ and $\mathbf{e_1}$ unchanged. This does not change the value of the integral, which now takes the form $$ \alpha^2\int_{\mathbb{R}^n} \langle x,\mathbf{e} \rangle ^2 e^{-\frac{1}{2} \|x\| ^2} = \alpha^2\int_{\mathbb{R}^n} x_1^2 e^{-\frac{1}{2} \|x\| ^2} $$ The integrand now is $x_1^2 \cdot e^{-x_1^2/2} \cdot e^{-x_2^2/2} \cdot \dots \cdot e^{-x_n^2/2}$. By Fubini's Theorem, therefore $$ \alpha^2\int_{\mathbb{R}^n} x_1^2 e^{-\frac{1}{2} \|x\| ^2} = \alpha^2 \left(2 \pi\right)^{n/2} $$
which is the expression that you derived.