Compute the moment generating function of $Y = X_1X_2 + X_1X_3 + X_2X_3$

Define $$A=\begin{bmatrix} 0&1&1\\ 1&0&1\\ 1&1&0 \end{bmatrix},\quad X=\begin{bmatrix} X_1\\ X_2\\ X_3 \end{bmatrix}$$ the other choices of $A$ would give the correct $Y$, but I need it to be symmetric later. Then $$ Y=\frac{1}{2}X^TAX $$ The pdf for 3d standard normal distribution $$ p(x)=(2\pi)^{-3/2}\exp\left[-\frac{1}{2}x^Tx\right] $$ The moment generating function for y is then $$ \mathbb{E}(e^{\lambda y})=\int d^3x(2\pi)^{-3/2}\exp\left[\frac{\lambda}{2}x^TAx\right]\exp\left[-\frac{1}{2}x^Tx\right]=(2\pi)^{-3/2}\int d^3x\exp\left[-\frac{1}{2}x^T(I-\lambda A)x\right] $$ The above integral is solved by the following for a real symmetric matrix $M$ $$ \int d^3x\exp\left[-\frac{1}{2}x^TMx\right]=\sqrt{\frac{(2\pi)^3}{\det M}} $$ Inserting we find $$ \mathbb{E}(e^{\lambda y})=\frac{1}{\sqrt{\det(I-\lambda A)}}=\frac{1}{\sqrt{-2\lambda^3-3\lambda^2+1}} $$


There may be some smarter ways to solve this, but repeated applications of the tower rule $$\mathbf E [X] = \mathbf E[\mathbf E[X|Y]]$$ will give you the result.

Take the definition of the moment generating function $$M_Y(t) = \mathbf E[\mathrm e^{t Y}]$$ and, in the first step, condition on $X_2$ and $X_3$; you get that $$\begin{aligned} M_Y(t) &= \mathbf E\big[ \mathbf E [ \mathrm e^{t Y} | X_2, X_3]\big]\\ &= \mathbf E\big[ \mathbf E [ \mathrm e^{t X_1(X_2 + X_3)} | X_2, X_3] \mathrm e^{t X_2X_3}\big] \end{aligned} $$ Note that, conditioned on $X_2$ and $X_3$, the random variable $X_1(X_2+X_3)$ is $N\big(0,(X_2+X_3)^2\big)$, so (using the definition of the MGF of a normal random variable) $$ \mathbf E [ \mathrm e^{t X_1(X_2 + X_3)} | X_2, X_3] = \mathrm e^{\frac{1}{2}(X_2 + X_3)^2t^2}$$ So we have that $$M_Y(t) = \mathbf{E} \big[ \mathrm e^{\frac{1}{2}(X_2 + X_3)^2t^2 + tX_2 X_3}\big].$$

Apply the same trick again to integrate out $X_2$: $$\begin{aligned} M_Y(t) &= \mathbf{E} \big[ \mathbf E[\mathrm e^{\frac{1}{2}(X_2 + X_3)^2t^2 + tX_2 X_3}|X_3]\big]\\ &=\mathbf{E} \big[ \mathbf E[\mathrm e^{\frac{1}{2}X_2^2 t^2 + X_2 X_3t^2 + tX_2 X_3}|X_3]\mathrm e^{\frac{1}{2} X_3^2t^2}\big]. \end{aligned}$$ The conditional expectation can be computed writing out the integral and completing the square at the exponent; you may need to prove the following intermediate result

Let $X\sim N(0, 1)$, then $$ \mathbf{E} [ \mathrm e^{ \frac{1}2 a X^2 + b X}] = \frac{1}{\sqrt{1-a}}\mathrm e^{\frac{1}{2}\frac{b^2}{1-a}}$$

Doing the same "complete the square" once more to integrate out $X_3$ should give you the final answer.