I am trying to solve a Heat transfer problem in a slightly ellipsoidal geometry using the Poission Equation, and Cauchy matching conditions on the boundary between the interior (which is a finite sphere) and exterior, where for the exterior only Laplace equation applies. But I seem to run into an inconsistency when performing the matching at this boundary. The problem setup is as follows (apologies if it is a bit lengthy, but I had to give sufficient context):
For the interior ($\rho$ is constant),
$ \nabla^2 \Psi(x,y,z) = \rho \ $
and a well-known solution in Rectangular coordinates is
$ \Psi(x,y,z) = \Psi_{0} - Ax^2-By^2-Cz^2$
For spherical symmetry, if we perform the usual coordinate transformations from Cartesian to Spherical coordinates:
$ x=r\sin\theta\cos\phi $
$ y=r\sin\theta\sin\phi $
$ z=r\cos\theta $
and also require that $ A=B=C $, we end up with a solution of the form
$ \Psi(x,y,z) = \Psi_{0} - Ar^2 $
which is what we would get if we directly solved the PDE in spherical polars. Similarly, the exterior solution (i.e. where $\rho=0$) is of the form
$ \Psi(x,y,z) = M - \frac{N}{r^2} $
and one can trivially match this at the spherical boundary. So far, so good.
But now instead of a sphere, consider a perturbed sphere, which is a tri-axial ellipsoid with small eccentricity. In this approximation, the Integration constants B and C can be approximated as $B=A+\epsilon$ and $C=A+\mu$, where the directional perturbations to the radius of the sphere ($\epsilon$ and $\mu$) are of a comparable magnitude. If we substitute these expressions into the above expression for the interior solution, keeping terms to linear order in $\epsilon$ and $\mu$, and then on the boundary of the interior we also substitute the spherical coordinate transformation, we end up with an expression of the form
$ \Psi(x,y,z)_{boundary} = \Psi_{0} - Ar^2 - \epsilon r^2\sin^2\theta sin^2\phi - \mu r^2 \cos^2\theta$
To make an easier identification with the solution of the exterior problem, both angular dependences can be written in terms of Spherical Harmonics, $Y_{lm}(\theta, \phi)$, where for the azimuthal dependence, $m=0,2$. So it would look something like (f, g, h, and l are numerical constants)
$ \Psi(x,y,z)_{boundary} = \Psi_{0} - Ar^2 - \epsilon r^2\ (fP_{0}(cos\theta) + gP_{2}(cos\theta)) (-2 + \exp(i2\phi) - \exp(-i2\phi)) - \mu r^2 (hP_{0}(cos\theta) + lP_{2}(cos\theta)) $
Now onto the exterior solution: The Laplace equation in Spherical polars can be solved in terms of inverse powers of r (the spherical part has already been shown above) and Spherical harmonics $Y_{lm}(\theta, \phi)$ containing Legendre Polynomials $P_{l}(\cos\theta)$, where there is a constraint on the suffix $m$ for the azimuthal dependence, namely that m cannot go beyond $+l$ or $-l$. As for choosing the suffixes $l$ and $m$, it all depends on the boundary limit of the interior solution. So using the last expression above to 'guide' the exterior solution, the perturbed part of the exterior solution must surely include contributions of the form
$ m(r)P_{0}(cos\theta) (-2 + \exp(i2\phi) - \exp(-i2\phi)) + n(r)P_{2}(cos\theta) (-2 + \exp(i2\phi) - \exp(-i2\phi)) $
But here is the problem that I seem to be running in to: According to the above selection rule on $m$ for the exterior solution, $m$ cannot go beyond $\pm l$. The second term ($P_{2}$) satifies the requirement, but the first term ($P_{0}$) clearly does not, and so cannot be included in the exterior solution. Yet, the boundary limit of the interior solution seems to demand it, otherwise the matching cannot be performed. What has gone wrong, and how is this resolved?