Why is the gravitational multipole expansion taken with the origin at the centre of mass instead of the centre of gravity? Wikipedia's article on the multipole expansion of the gravitational potential expands it as
$$
 V(\mathbf{x}) = - \frac{GM}{|\mathbf{x}|} - \frac{G}{|\mathbf{x}|} \int \left(\frac{r}{|\mathbf{x}|}\right)^2 \frac {3 \cos^2 \theta - 1}{2} dm(\mathbf{r}) + \cdots
$$
with respect to a system of reference with the origin at the centre of mass $\mathbf r_\mathrm{COM}=\int\mathbf r\: \mathrm dm(\mathbf{r})$ of the mass distribution.
I don't understand this choice of system of reference. Why is the origin taken at the centre of mass? To me it seems to make more sense to put the origin at the centre of gravity, where the gravitational potential is 0. In the center of mass, it may not be!
 A: The multipole expansion in its most general form reads
$$
V(\mathbf r)=\sum_{l=0}^\infty \sum_{m=-l}^lQ_{lm}\frac{Y_{lm}(\theta,\phi)}{r^{l+1}},
$$
where the $Q_{lm}$ are the multipole moments of the system and the $Y_{lm}$ are spherical harmonics, and in this form it is applicable to bodies of any shape, and located at any point in space, so long as the body is contained in a finite spherical region of radius $R$ and the test point $\mathbf r$ is located outside this spherical region.
The Wikipedia article you linked to is concerned with the first few terms in this expansion, which can be rewritten as
$$
V(\mathbf r)=-\frac{GM}{r}+\frac{\mathbf d\cdot\mathbf r}{r^3}+\frac{p_2(x,y,z)}{r^5},
$$
where $p_2$ is a traceless homogeneous quadratic polynomial, and the dipole moment of the system is defined as
$$
\mathbf d=G\int\mathbf r\:\mathrm dm(\mathbf r). \tag 1
$$
These first three terms go down with the distance to the origin $r$ as $1/r$, $1/r^2$ and $1/r^3$ respectively, which means that the relative importance of the second, dipole term goes as $1/r$ with respect to the first, monopole term.
This changes, however, if you set up your origin at the centre of mass of the distribution, since then the dipole moment $(1)$ vanishes. If you do this, then you can completely ignore the dipole term, and the subleading term of the expansion becomes the quadrupole term, which goes as $1/r^2$ of the monopole term, which makes the truncated series much more accurate.
That's about it, really.

In addition, though, you should note that the centre of gravity as you've linked to cannot be defined for a body in isolation: it is a property of a mass distribution $\mathrm dm(\mathbf r)$ and an externally-produced gravitational field $\mathbf g(\mathbf r)$ acting on the distribution. In these conditions you can find the total force
$$
\mathbf F=\int\mathbf g(\mathbf r)\mathrm dm(\mathbf r)
$$
and the total torque with respect to an arbitrary origin $\mathbf r_0$,
$$
\boldsymbol\tau=\int(\mathbf r-\mathbf r_0)\times\mathbf g(\mathbf r)\mathrm dm(\mathbf r),
$$
and then you can define the centre of gravity as the $\mathbf r_0$ at which $\boldsymbol \tau$ vanishes. If you don't have the external field $\mathbf g(\mathbf r)$, though, this doesn't make any sense: the centre of gravity is the point through which a given field configuration can be thought of as acting, and it depends on the field configuration.
A: Unlike the centre of masss the centre of gravity is not an intrinsic property of a mass distribution. The centre of gravity only becomes distinct from the centre of mass when there is an external gravitational field that is non-uniform i.e. varies significantly across the distribution of masses. If the external gravitational field is higher on one side of the mass distribution than the other it effectively pulls the centre of mass towards the higher gravity side.
The multipole calculation you describe involves no external gravitational field so the centre of gravity coincides with the centre of mass.
A: I think you are puzzled because you are not distinguishing the gravitational potential produced by the system and the one acting on the system.
The multipole expansion you link to, is performed on the expression for the potential produced by the system on a distant point. In this case the system has a mass center, but no gravity center can be defined for it without having an external gravitational field acting on it. Also the resulting field will be always along the line passing by the center of mass, that is why is convenient here.
But when you have an external gravitational field acting over some system, you can define a gravity center for it, as the point which would not "feel" any torque regardless of the system orientation in the external field, or just using the formula you linked to. But in this case it might not coincide with the mass center because the field might be inhomogenous throughout the system, so different parts can feel different forces.
Summing up, the choice is based on the symmetry of the potential (the field produced by the system always points towards its mass center) and because no gravity center can be defined w/o external field.
