Dielectric sphere in a static field

I'll try to give the simplest solution I can think of to the classical problem of a dielectric sphere in a constant external field, see for instance  Landau & Lifshitz, vol. 8, chapter II, §8.

The sphere radius is $R$ and its dielectric constant is $\epsilon_i$, while that of the surrounding medium is $\epsilon_e$. The electric field at infinity is along the $z$ axis: $\mathbf{E}_0 = E_0\, \hat{z}$. Considering the symmetry of the object we will work in spherical coordinates $(r,\theta,\phi)$. Before writing down any equations, let us note the following points:

1. The system has rotational symmetry about the $z$ axis. Thus, the field has no component along $\hat{\phi}$.
2. The problem is antisymmetric with respect to $z$ (changing $E_0$ to $-E_0$ reverses the sign of the field everywhere).
3. The field scales by $E_0$ (doubling the applied field doubles the field at any point in the system).
4. The only length parameter is the sphere radius $R$.

The last two points express scaling relations, while the first two relate to the symmetry properties of the "complete" system: sphere + applied field and we can discuss them in terms of the multipole expansion of the potential created by a charge distribution:
$$\begin{equation} \Phi(r, \theta, \phi) = \frac{1}{4 \pi \epsilon_0} \sum_{\ell = 0}^{\infty} \frac{1}{r^{\ell + 1}} \sum_{m = - \ell}^{m = \ell} Q_{\ell m} Y_{\ell m}(\theta, \phi) \label{eq:finalpot} \end{equation}$$
Before doing that, we should notice that a constant field along $z$ is necessary for describing the potential outside the sphere (and is a trivial solution to Laplace's equation). However, such a field is not accounted for in the expansion, since it does not decay at infinity, and we must consider it separately. Aside from that, condition 1 amounts to forbidding any azimuthally modulated contribution, so that only terms with $m=0$ should be preserved. Condition 2 further limits the sum over $\Phi_{\ell}$ to odd values of $\ell$. For instance, a net charge $q$ at the origin (corresponding to $\ell = 0$) is symmetric, so that it cannot be engendered by the source field $\mathbf{E}_0$; thus, $q=0$.

Instead of the potential  $\Phi$, we will work directly with the field $\mathbf{E}$, written as the sum of the constant field along $z$ and a dipole field ($\ell = 1$) along the same direction. We will show that those terms are enough to fulfill the boundary conditions at the surface of the sphere; consequently, they describe the solution completely and terms of higher order than the dipole ($\ell = 3, 5, \ldots$) are not necessary.

Outside the sphere, the constant field is simply $\mathbf{E}_0$, to fulfill the boundary conditions at infinity. As for the dipole field, from the considerations above (1) it is parallel to $\hat{z}$ and proportional to $E_0$ (3); furthermore, we will rescale the distance $r$ by the sphere radius $R$ (4):
$$\begin{equation} \mathbf{E}_e(r,\theta) = \mathbf{E}_0 + A E_0 \left ( \frac{R}{r} \right )^3 \left [ 3 \cos (\theta) \hat{r} - \hat{z} \right ] \label{eq:Eext} \end{equation}$$
Inside the sphere the second term vanishes (otherwise, it would diverge at the origin), so we are left with a constant field:
$$\begin{equation} \mathbf{E}_i(r,\theta) = B \, \mathbf{E}_0 = B E_0\, \hat{z} \label{eq:Eint} \end{equation}$$
The constants $A$ and $B$ must be determined from the boundary conditions at the sphere surface:
$$\begin{equation} \begin{split} {E}_i^{\|}(R, \theta) &= {E}_e^{\|}(R, \theta) \qquad \mbox{and}\\ \epsilon_i{E}_i^{\bot}(R, \theta) &= \epsilon_e{E}_e^{\bot}(R, \theta) \end{split} \label{eq:BCsphere} \end{equation}$$

where subscripts $\|$ and $\bot$ denote the field components parallel and normal to the interface, respectively. Since the field has no component along $\hat{\phi}$ and the normal to the surface of the sphere is always $\hat{r}$, $\mathbf{E} = E^{\bot} \hat{r} + E^{\|} \hat{\theta}$. We therefore plug equations $\eqref{eq:Eext}$ and $\eqref{eq:Eint}$ in $\eqref{eq:BCsphere}$, using the identity $\hat{z} = \hat{r} \cos (\theta) - \hat{\theta} \sin (\theta)$. We are left with the system of linear equations:

$$\begin{equation} \begin{split} (1-A) &= B\\ \epsilon_e (1+2A) &= \epsilon_i B \end{split} \label{eq:BCcoeff} \end{equation}$$
which yields $\displaystyle A = \frac{\epsilon_r - 1}{\epsilon_r + 2}$ and $\displaystyle B = \frac{3}{\epsilon_r + 2}$, with $\epsilon_r = \epsilon_i / \epsilon_e$ the relative dielectric constant of the sphere with respect to the surrounding medium. Finally:
$$\begin{equation} \begin{split} \mathbf{E}_e(r,\theta) &= \mathbf{E}_0 + E_0 \frac{\epsilon_r - 1}{\epsilon_r + 2} \left ( \frac{R}{r} \right )^3 \left [ 3 \cos (\theta) \hat{r} - \hat{z} \right ]\\ \mathbf{E}_i(r,\theta) &= \frac{3}{\epsilon_r + 2} \mathbf{E}_0 \end{split} \label{eq:FinalE} \end{equation}$$
We see that the field does not depend on $\epsilon_i$ and $\epsilon_e$ separately, but only on their ratio $\epsilon_r$; furthermore, for $\epsilon_r = 1$ ($\epsilon_i = \epsilon_e$) the field $\mathbf{E}(r,\theta) = \mathbf{E}_0$ everywhere; the sphere is "transparent" and induces no perturbation to the applied field.