$\newcommand\R{\mathbb R}$Letting $I:=\mathbf I$, $a:=\mathbf xi$, $b:=\mathbf xj$, $N:=(4M+I)^{1/2}$, $c_1:=a^\top Ma+b^\top Mb$, $m:=2N^{-1}M(a+b)$, $c_2:=m^\top m/2=2(a+b)^\top M(4M+I)^{-1}M(a+b)$, completing the squares, using the substitutions $z=Ny$ and $y-m=v$, 
and doing a bit of algebra, we see that the expectation in question is 
$$\frac{e^{c_2-c_1}}{\det(4M+I)^{1/2}}(N^{-2}+N^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top). \tag{1}\label{1}$$
Here one can get rid of $N$ and $m$ by noting that $N^{-2}=(4M+I)^{-1}$, so that $N^{-1}m=2(4M+I)^{-1}M(a+b)$ and hence $m^\top N^{-1}=2(a+b)^\top M(4M+I)^{-1}$. 

----

**Details:** The expectation in question is 
$$J:=(2\pi)^{-n/2}\int_{\R^n}dz\,(z-a)(z-b)^\top e^{g(z)}, \tag{2}\label{2}$$
where 
\begin{align}
g(z)&:=-(z-a)^\top M(z-a)-(z-b)^\top M(z-b)-z^\top z/2 \\ 
&=-z^\top(2M+I/2)z+2(a+b)^\top Mz-c_1 \\ 
&=-c_1-y^\top y/2+m^\top y \\ 
&=-c_1+c_2-v^\top v/2. \tag{3}\label{3} 
\end{align}
Also, 
\begin{align}
(z-a)(z-b)^\top&=(N^{-1}(v+m)-a)(N^{-1}(v+m)-b)^\top \\ 
&=h(v):=N^{-1}(v+m)(v+m)^\top N^{-1}+ab^\top \\ 
&\qquad\qquad-a(v+m)^\top N^{-1}-N^{-1}(v+m)b^\top \\ 
&=N^{-1}vv^\top N^{-1}+N^{-1}mm^\top N^{-1} \\ 
&+N^{-1}vm^\top N^{-1}+N^{-1}mv^\top N^{-1} \\ 
&+ab^\top -a(v^\top+m^\top) N^{-1}-N^{-1}(v+m)b^\top \\
\tag{4}\label{4}
\end{align}
and $dz=dy/\det N=dv/\det N=dv/\det(4M+I)^{1/2}$. 

So, using \eqref{2}, \eqref{3}, and \eqref{3}, and letting $V$ denote a standard normal random vector in $\R^n$, we see that 
\begin{align}
J&=\frac{e^{-c_1+c_2}}{\det(4M+I)^{1/2}}Eh(V). \tag{5}\label{5}
\end{align}
Noting finally that $EV=0$ and $EVV^\top=I$, from \eqref{5} we get the expression \eqref{1} for the expectation $J$ in question.