First, I suppose that the discontinuous $\epsilon_r$ is not a significant problem, because the discontinuous material properties can be expressed in terms of the Heaviside step function, and the Heaviside step function has its own derivative; this does not seem to be a difficulty.
For example, see here, where the magnetic permeability is discontinuous.
In any case, I will be writing about the electrostatic problem
where there are dielectric discontinuities.
The target equation is given as $\nabla\cdot\vec{D}=\rho$.
As usual, we include the scalar electric potential $\phi$ as the fundamental variable.
\begin{equation}
-\nabla\cdot\left[\epsilon(x,y,z)\nabla\phi(x,y,z)\right]-\rho(x,y,z)=0,\tag{1}
\end{equation}
where $\epsilon(x,y,z)$ is the coordinate dependent electric permittivity (positive everywhere)
and the charge densitiy distribution $\rho(x,y,z)$ is assumed to be given.
Let $\Omega$ be the target domain, and let Dirichlet and Neuman boundary conditions be specified at some outer boundary of $\Omega$.
The outer boundary surface of the domain is represented as $\partial\Omega$.
\begin{equation}
\begin{split}
\phi&=\text{prescribed value on }\partial\Omega_D \\
\frac{\partial \phi}{\partial n}&=0\text{ on }\partial\Omega_N\\
\partial\Omega&=\partial\Omega_D\cup\partial\Omega_N,\;\;\partial\Omega_D\cap\partial\Omega_N=\text{null}
\end{split}\tag{2}
\end{equation}
And we add the condition that at the outer boundary surface,
$\epsilon$ is a constant and its value is in vacuum (all dielectric materials are inside the domain).
\begin{equation}
\epsilon=\epsilon_0 \text{ on }\partial\Omega \tag{3}
\end{equation}
Let's start with the proof of uniqueness.
Suppose there are 2 different solutions $\phi_1$ and $\phi_2$ to (1) and (2) and let us define the difference as
\begin{equation}
\phi_D:=\phi_2-\phi_1\tag{4}
\end{equation}
Since $\phi_1$ and $\phi_2$ satisfy equation (1), the difference $\phi_D$ satisfies
\begin{equation}
\nabla\cdot\left[\epsilon(x,y,z)\nabla\phi_D(x,y,z)\right]=0.\tag{5}
\end{equation}
Consider solving this equation.
To introduce the weak convergence argument,
we have an arbitrary (i.e. all) weight function $\delta\phi$.
Think of $\delta\phi$ here as a one-digit quantity, not the product of two.
Note that $\delta\phi$ must vanish at the Dirichlet boundary:
\begin{equation}
\delta\phi=0\text{ on }\partial\Omega_D\tag{6}
\end{equation}
Multiply the weight function by equation (5) and integrate by volume in the domain.
\begin{equation}
\int_{\Omega}\delta\phi\nabla\cdot\left[\epsilon(x,y,z)\nabla\phi_D(x,y,z)\right]dV=0.\tag{7}
\end{equation}
Apply the partial integration to equation (7) as we always do.
\begin{equation}
-\int_{\Omega}\nabla\delta\phi\cdot\left[\epsilon(x,y,z)\nabla\phi_D\right]dV
+\int_{\partial\Omega}\epsilon_0\delta\phi\frac{\partial\phi_D}{\partial n}\cdot\vec{n}dS
=0.\tag{8}
\end{equation}
The second term of (8), which is the surface integration, is zero.
This is because of equation (6) and the Neumann boundary condition.
Since $\delta\phi$ is arbitrary, we can choose $\delta\phi=\phi_D$.
Substituting this into (8) we get,
\begin{equation}
\begin{split}
\int_{\Omega}\epsilon(x,y,z)(\nabla\phi_D)^2dV=0. \\
\therefore \phi_D=\phi_1-\phi_2=0
\end{split} \tag{9}
\end{equation}
This proves the uniqueness of $\phi$ in the case that dielectric materials are included.
Returning to your concern, that is,
the problem of a dielectric sphere in a uniform electric field.
It is when the dielectric constant is given as
\begin{equation}
\epsilon(x,y,z)=\epsilon_0\left[\epsilon_r1_{\text{IN}}(x,y,z)+1_{\text{OUT}}(x,y,z)\right]\tag{10},
\end{equation}
Here $1_{\text{IN}}(x,y,z)$ and $1_{\text{OUT}}(x,y,z)$ are simple extensions of Heaviside's step function,
which is a function of 1 inside the sphere and 0 outside, or 0 inside the sphere and 1 outside.
This case is included in the above described uniqueness theorem.
