The starred step is valid because (a) p and q have the same zeroth and second moments and (b) log(p) is a polynomial function of the components of x whose terms have total degrees 0 or 2.
You need to know only two things about a multivariate normal distribution with zero mean:
log(p) is a quadratic function of x=(x1,x2,…,xn) with no linear terms. Specifically, there are constants C and pij for which log(p(x))=C+n∑i,j=1pijxixj.
(Of course C and the pij can be written in terms of Σ, but this detail does not matter.)
Σ gives the second moments of the distribution. That is, Σij=Ep(xixj)=∫p(x)xixjdx.
We may use this information to work out an integral:
It breaks into the sum of two parts:
∫(q(x)−p(x))Cdx=C(∫q(x)dx−∫p(x)dx)=C(1−1)=0, because both q and p are probability density functions.
∫(q(x)−p(x))∑ni,j=1pijxixjdx=∑ni,j=1pij∫(q(x)−p(x))xixjdx=0 because each of the integrals on the right hand side, ∫q(x)xixjdx and ∫p(x)xixjdx, has the same value (to wit, Σij). This is what the remark “yield the same moments of the quadratic form” is intended to say.
The result follows immediately: since ∫(q(x)−p(x))log(p(x))dx=0, we conclude that ∫q(x)log(p(x))dx=∫p(x)log(p(x))dx.