In the region between two spherical plates, we have,
$$\begin{array}{l}
\nabla \times E = 0\,\,\, (1)\\
\nabla .\left( {\varepsilon \left( \theta \right)E} \right) = 0\,\,\, (2)
\end{array}$$
The first equation leads to the definition of scalar potential $\Phi$, i.e. $E = - \nabla \Phi $. Therefore, from Eq. (2) we have,
$$\nabla .\left( {\varepsilon \left( \theta \right)\nabla \Phi } \right) = \nabla \varepsilon \left( \theta \right).\nabla \Phi + \varepsilon \left( \theta \right){\nabla ^2}\Phi = 0 \,\,\, (3)$$
Since the problem is invariant under $\varphi$-rotation, $\Phi=\Phi(r,\theta)$. Using the separation of variable method, $\Phi(r,\theta)=R(r)\Theta(\theta)$. By substituting this function into the Eq. (3), after some straightforward calculations (in spherical coordinate system), one can easily obtain the following equation,
$$-\frac{1}{{R(r)}}\frac{d}{{dr}}\left( {{r^2}\frac{dR(r)}{dr}} \right) = {\frac{1}{{\varepsilon \left( \theta \right)\Theta \left( \theta \right)}}\frac{{d\varepsilon \left( \theta \right)}}{{d\theta }}\frac{{d\Theta \left( \theta \right)}}{{d\theta }} + \frac{1}{{\sin \left( \theta \right)\Theta \left( \theta \right)}}\frac{d}{{d\theta }}\left( {\sin \left( \theta \right)\frac{{d\Theta \left( \theta \right)}}{{d\theta }}} \right)}
$$
One side of this equation is r-dependent and the other side is $\theta$-dependent. Therefore both side must be a constant value, namely $\lambda$,
$$\frac{d}{{dr}}\left( {{r^2}\frac{dR(r)}{dr}} \right) + \lambda R(r) = 0\,\,\, (4)$$
$$\lambda = \frac{1}{{\varepsilon \left( \theta \right)\Theta \left( \theta \right)}}\frac{{d\varepsilon \left( \theta \right)}}{{d\theta }}\frac{{d\Theta \left( \theta \right)}}{{d\theta }} + \frac{1}{{\sin \left( \theta \right)\Theta \left( \theta \right)}}\frac{d}{{d\theta }}\left( {\sin \left( \theta \right)\frac{{d\Theta \left( \theta \right)}}{{d\theta }}} \right)\,\,\, (5)$$
On the other hand, we know that a conductor has an equipotential surface , i.e. the functions $\Phi(r=r_1,\theta)$ and $\Phi(r=r_2,\theta)$ must be independent of $\theta$. Therefore, using $\Phi(r,\theta)=R(r)\Theta(\theta)$, we conclude that $\Theta(\theta)$ must be constant, namely one, and therefore
$$\Phi(r,\theta)=R(r)\Theta(\theta)=R(r)$$
is independent of $\theta$ for all values of $r$. This is the desired result.
In addition, in this case, $\frac{d\Theta(\theta)}{d\theta}=0$, from Eq. (5), we find $\lambda=0$, and therefore by using Eq. (4), we have
$$\frac{d}{{dr}}\left( {{r^2}\frac{{dR(r)}}{{dr}}} \right) = 0$$
with the solution,
$$R(r) = {c_1} + \frac{{{c_2}}}{r} \,\,\, (6)$$
where, $c_1$ and $c_2$ can be determined from the boundary conditions. (Here the total charge at $r=r_1$ is $Q$ and the total charge at $r=r_2$ is $-Q$).
By setting the reference point at infinity, i.e. by assuming $\Phi(r=\infty)=0$, we find $c_1=0$. For calculation of $c_2$ we use this fact that the total charge at $r=r_1$ is $Q$ (or, equivalently, the total charge at $r=r_2$ is $-Q$),
$$\int {\sigma \left( \theta \right)da} = \int {\sigma \left( \theta \right)r_1^2\sin } \left( \theta \right)d\theta d\varphi = Q \,\,\, (7)$$
Consider the general boundary condition
$$\left( {{\varepsilon _2}{E_2} - {\varepsilon _1}{E_1}} \right).{\hat n_{12}} = \sigma$$
where $\hat n_{12}$ is a unit vector from region 1 to region 2. In this problem the region 1 is inside the spherical plate at $r=r_1$ and the region 2 is its outside surface. Since the electric field inside a conductor is zero, we have
$$\sigma(\theta) = {\left. {\varepsilon \left( \theta \right){E_r}} \right|_{r = {r_1}}} = - \varepsilon \left( \theta \right){\left. {\frac{{d\Phi \left( r \right)}}{{dr}}} \right|_{r = {r_1}}} = - \varepsilon \left( \theta \right){\left. {\frac{{dR\left( r \right)}}{{dr}}} \right|_{r = {r_1}}} = \frac{{{c_2}}}{{{r_1}^2}}\varepsilon \left( \theta \right) \,\,\, (8)$$
Substituting Eq. (8) in Eq. (7), gives,
$${c_2} = \frac{Q}{{\int {\varepsilon \left( \theta \right)\sin } \left( \theta \right)d\theta d\varphi }} = \frac{Q}{{2\pi \int_0^\pi {\varepsilon \left( \theta \right)\sin \left( \theta \right)d\theta } }}\,\,\, (9)$$
which can be obtained easily.
Another Solution:
If, at beginning, we know that the electric potential is independent of $\theta$, we can use the Gauss's law,
$$\int {\varepsilon \left( \theta \right)E.\hat nda = Q} \,\,\, (10)$$
In this case since $\Phi$ is independent of $\theta$, $E = - \nabla \Phi = - \frac{{d\Phi (r)}}{{dr}}\hat r = e(r)\hat r$ is also independent of $\theta$ and radial. Therefore, using Eq. (10), we have,
$$e(r){r^2}\int {\varepsilon \left( \theta \right)\sin \left( \theta \right)d\theta d\varphi = Q} $$
and therefore (see Eq. (9)),
$$e(r) = \frac{Q}{{{r^2}\int {\varepsilon \left( \theta \right)\sin \left( \theta \right)d\theta d\varphi } }} = \frac{Q}{{2\pi {r^2}\int_0^\pi {\varepsilon \left( \theta \right)\sin \left( \theta \right)d\theta } }} = \frac{{{c_2}}}{{{r^2}}}\,\,\, (11)$$
Integrating the differential equation $- \frac{{d\Phi (r)}}{{dr}} = e(r) = \frac{{{c_2}}}{{{r^2}}}$, one obtain,
$$\Phi (r) = {c_1} + \frac{{{c_2}}}{r}$$
which is exactly the previous result ($c_1=0$ if we set the reference point at infinity).
