Determining vector fields in a Spherical Capacitor with Quadrant Dielectrics

I am working with a spherical capacitor consisting of an inner radius R1 and an outer radius R2. The region between the plates of the capacitor is filled with four dielectric materials, each shaped like a spherical wedge with an angle of π/2. These materials, occupying equal quadrants, extend from the inner radius R1 to the outer radius R2. Each dielectric is linear, homogeneous, isotropic, and has permittivities ε1, ε2, ε3, and ε4, respectively.

How do I determine the vectors D (electric displacement), E (electric field), and P (polarization) between the plates?

I would start by calculating the electric field because it only depends on the total charge, so by using Gauss's theorem and due to the spherical symmetry of the problem, it is easily obtained (I think)

$$E = \frac{Q}{4\pi \epsilon_0 R^2}$$

From here I suppose that D and P are straightforward, but I'm not sure.

Spherical coordinates are used here; let $$r$$ be radial coordinate, $$\theta$$ the polar angle and $$\phi$$ the azimuthal angle. Let $$V$$ denote the electric scalar potential; the unit of $$V$$ is volte. Although the coordinate dependence of $$V$$ is $$V(r,\theta,\phi)$$, let us restrict ourselves to the case where $$V$$ only depends on $$r$$ for a while.
Assume that the electric scalar potential is given as $$\begin{equation*} V(r)=k_1\frac{1}{r}+k_2\tag{1} \end{equation*}$$ where $$k_1$$ and $$k_2$$ are constants, given by, \begin{align*} k_1=-\frac{R_2R_1(V_2-V_1)}{R_2-R_1} \\ k_2=\frac{R_2V_2-R_1V_1}{R_2-R_1} \tag{2} \end{align*} Equations (1) and (2) satisfy the boundary conditions upon metal surfaces: \begin{align*} V(R_1)=V_1 \\ V(R_2)=V_2.\tag{3} \end{align*} The electric field in this case has only r-component and given as $$\begin{equation*} E_r(r)=-\mathrm{grad}V(r)=-\frac{k_1}{r^2}\tag{4} \end{equation*}$$
Up to this point, we have assumed that the physical quantity depends only on $$r$$. From this point on, we will change the equation to depend on other coordinate variables as well. Electric flux density ($$\mathbf{D}$$), has r-component: \begin{align*} D_r(r,\theta,\phi)&=\epsilon(\theta,\phi) E_r(r)=\epsilon(\theta,\phi)\left(-\frac{\partial V(r)}{\partial r}\right)=-\epsilon(\theta,\phi)\frac{k_1}{r^2} \\ &=\epsilon(\theta,\phi)\frac{R_2R_1(V_2-V_1)}{R_2-R_1}\frac{1}{r^2}\tag{5} \end{align*} As a special case, let us focus on the case $$r = R_1$$. Electric flux density becomes \begin{align*} D_r(r=R_1,\theta,\phi) =\epsilon(\theta,\phi)\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}.\tag{6} \end{align*} Total electric charge on $$r=R_1$$ metal surface is \begin{align*} Q=\int_{S}D_r(r=R_1,\theta,\phi)dS =\int_{S}\left\{\epsilon(\theta,\phi)\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}\right\}dS.\tag{7} \end{align*} where $$dS$$ is the small area at $$r=R_1$$ sphere (in unit of $$m^2$$).
If in the particular case that permitivity is uniform, i.e.,$$\epsilon(\theta,\phi)=\epsilon_1$$, then \begin{align*} Q=\epsilon_i\frac{R_2(V_2-V_1)}{R_1(R_2-R_1)}4\pi R_1^2 =\epsilon_i\frac{4\pi R_1R_2}{R_2-R_1}(V_2-V_1)\tag{8} \end{align*} From this equation, the capacity matrix is found to be: $$\begin{equation*} C=\epsilon_i\frac{4\pi R_1R_2}{R_2-R_1}.\tag{9} \end{equation*}$$ For another case that the questioner is focusing on, the capacity is the following equation, $$\begin{equation*} C=\frac{\epsilon_1+\epsilon_2+\epsilon_3+\epsilon_4}{4}\frac{4\pi R_1R_2}{R_2-R_1}.\tag{10} \end{equation*}$$
Now let me justify that equations (1), (2), (4), (5) and (6) are the solution for the case where the permitivity depends on $$\theta,\phi$$: $$\begin{equation*} \epsilon=\epsilon(\theta,\phi)\;\;\;\text{if }R_1 The general starting electro-static field equation is as follows. \begin{align*} \mathrm{div}\mathbf{D}=\rho,\tag{12} \\ \mathbf{D}=\epsilon\mathbf{E},\tag{13} \\ \mathbf{E}=-\mathrm{grad}V.\tag{14} \end{align*} The symbols in (12) to (14) are according to the standard conventions. The electric free charge density $$\rho$$ is $$0$$ in the insulating materials and $$\rho \neq 0$$ at the metal surfaces. Thus within insulationg region, (12) becomes \begin{align*} \mathrm{div}\left(\epsilon\mathbf{E}\right) &=\mathrm{grad}\epsilon\cdot\mathbf{E}+\epsilon\mathrm{div}\mathbf{E} \\ &=\frac{1}{r}\frac{\partial\epsilon}{\partial\theta}E_{\theta} +\frac{1}{r\sin\theta}\frac{\partial\epsilon}{\partial\phi}E_{\phi} +\epsilon\left\{ \frac{1}{r^2}\frac{\partial(r^2E_r)}{\partial r} +\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(E_{\theta}\sin\theta) +\frac{1}{r\sin\theta}\frac{\partial E_{\phi}}{\partial\phi}\right\} \\ &=\frac{1}{r}\frac{\partial\epsilon}{\partial\theta}\left(-\frac{1}{r}\frac{\partial V}{\partial\theta}\right) +\frac{1}{r\sin\theta}\frac{\partial\epsilon}{\partial\phi}\left(-\frac{1}{r\sin\theta}\frac{\partial V}{\partial\phi}\right) \\ &+\epsilon\left\{ \frac{1}{r^2}\frac{\partial\left(r^2\left(-\frac{\partial V}{\partial r}\right)\right)}{\partial r} +\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\left(\left(-\frac{1}{r}\frac{\partial V}{\partial\theta}\right)\sin\theta\right) +\frac{1}{r\sin\theta}\frac{\partial}{\partial\phi}\left(-\frac{1}{r\sin\theta}\frac{\partial V}{\partial\phi}\right)\right\}\\ &=0\tag{15} \end{align*} Equation (15) is the equation to be solved, although it is complex, with V as the unknown. Starting from equation (15), trying to find V is difficult, but substituting (1) into (15) shows that it is one of the solutions. On the other hand, someone has proved the uniqueness of the solution, so we can conclude that (1) is the solution of (15).