The gravitational potential of ellipsoid In the literature (Kirchhoff G. - Mechanic (1897), Lecture 18 or Lamb, H. - Hydrodynamics (1879)) one can find the following analytical closed form expression for the gravitational potential of homogeneous ellipsoid of unit density, whose surface is given by 
\begin{equation}
 \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1 \;.
\end{equation}
Gravitational point for internal points is
\begin{equation}
 \Omega=\pi abc\int_0^\infty\left(1-\frac{x^2}{a^2+\lambda}-\frac{y^2}{b^2+\lambda}-\frac{z^2}{c^2+\lambda}\right)\frac{d\lambda}{\Delta}
\end{equation}
and for external points
\begin{equation}
 \Omega=\pi abc\int_u^\infty\left(1-\frac{x^2}{a^2+\lambda}-\frac{y^2}{b^2+\lambda}-\frac{z^2}{c^2+\lambda}\right)\frac{d\lambda}{\Delta} \;,
\end{equation}
where
\begin{equation}
 \Delta=\sqrt{(a^2+\lambda)(b^2+\lambda)(c^2+\lambda)}
\end{equation}
and $u$ is the positive root of equation
\begin{equation}
 \frac{x^2}{a^2+u}+\frac{y^2}{b^2+u}+\frac{z^2}{c^2+u}=1 \;.
\end{equation}
The expressions in these formulas appear similar to confocal ellipsoidal coordinates.
How can these formulas be derived? (perhaps something more readable than the original papers)
Can they be derived in terms of ellipsoidal harmonics?
 A: This is a very late answer. There is the book
Ellipsoidal Figures of Equilibrium
by the god himself in this field S. Chandrasekhar. Chapter 3 is devoted fully to understanding the gravitational potentials of ellipsoids. Theorems 3 and 9 are what you are looking for.
You also asked

Can they be derived in terms of ellipsoidal harmonics?

Chandrasekhar does not derive the equations in terms of ellipsoidal harmonics. In fact, he states that very early on in the introduction (section 16). Instead he employs spherical polar coordinates and proceeds by establishing a series of lemmas on the moments of the mass distribution. This amounts to considering integrals of the form
$$I(u) = a_1 a_2 a_3 \int_u^{\infty} \frac{du}{\Delta}; \qquad A_i(u) = a_1 a_2 a_3 \int_u^{\infty} \frac{du}{\Delta (a_i^2 + u)}$$
where $\Delta^2=(a_1^2+u) (a_2^2+u) (a_3^2+u)$ and $a_i$ are the semi-major axes of the ellipsoid. Then come the two theorems you need
Theorem 3: At a point $x_i$ interior to the ellipsoid, the potential is
$$\Phi = \pi  G \rho \Big[I(0) - \sum_{i=1}^3 A_i(0) x_i^2 \Big]$$
Theorem 9: At a point $x_i$ exterior to the ellipsoid, the potential is
$$\Phi = \pi  G \rho \Big[I(\lambda) - \sum_{i=1}^3 A_i(\lambda) x_i^2 \Big]$$
where $\lambda$ is the (algebraically) largest root of
$$\sum_{i=1}^3 \frac{x_i^2}{a_i^2 + \lambda} = 1$$
A: There is a nice and clear derivation in 'The theory of potential', W. D. MacMillan, Dover Publications, New York, 1958. 
In paragraph 32 you can find the derivation of the potential at an interior point, and in paragraph 35 an elegant derivation of the potential at an exterior point. 
The latter uses the fact that the gradient of the potential at an exterior point can be related to that at an interior point by simple considerations. Given that the potential gradient at an interior point is known, the one at an exterior point can thus be obtained. Finally, the derivation shows how to obtain the potential that produces that gradient. 
I do not know about the derivation in terms of ellipsoidal harmonics. 
