Background: Given that equilibrium solutions of charges in a spherical shell formed from a conductor satisfy the Poisson problem $\nabla^2 \Phi = - \frac{\rho}{\epsilon_0}$ as a corollary of the differential form of Gauß's Law $\nabla \cdot E = \frac{\rho}{\epsilon_0}$, would it be natural to model the events proceeding the equilibrium with the diffusion equation?
To make my question more rigorous, we prescribe the system with initial charge density $\rho(r, \theta, \varphi)$ that is not uniform. Without using a Galerkin method, does there exist an analytical solution for $\Phi(r, \theta, \varphi, t)$?
Pardon me if my intuition is grossly over simplified; however, could the transition to an equilibrium state be viewed as a progression to positive eigenvalues? Taking the diffusion pde in $\mathbb{R}^3$ ($\partial_t u = k \nabla^2 u$) to be parabolic with degenerate eigenvalues, would not the equilibrium in the Poisson equation represent an elliptic PDE ($\nabla^2 \Phi = - \frac{\rho}{\epsilon_0}$)? Additionally, the steady state implies that $\partial_t u = 0$, which enforces the requirement of being harmonic on the solution.
Question: Would the following initial value problem suffice for modeling the redistribution of the charge density and voltage to reach equilibrium. My system is underconstrained, which leads me to believe that I am neglecting some physical relationship.
$$\partial_t \Phi = \frac{-k \rho(r, \theta, \varphi, t)}{\epsilon_0}$$
$$\rho(r, \theta, \varphi, 0) = \rho_0(r, \theta, \varphi) : \mathbb{R} \times [0, 2 \pi) \times [0, 2 \pi) \to \mathbb{R}$$
I am struggling to find the evolution of $\rho(r, \theta, \varphi, t)$ in time. I understand that areas to which the $\vec E$-field lines converge will have increasing charge density, and the opposite will hold for the case of diverging lines. Could I take the divergence of the following expression for the $\vec E$-field? Would some clever application of Gauß's law suffice instead of evaluating the following integral?
$$\vec E(r, \theta, \varphi) = \int_0^{\pi} \int_0^{2 \pi} \int_{r_0}^{r_1} \frac{dq}{4 \pi \epsilon_0 d((r^*, \theta^*, \varphi^*), (r, \theta, \varphi))^{\frac{3}{2}}} \vec d((r^*, \theta^*, \varphi^*), (r, \theta, \varphi)) dr^* d \theta^* d\varphi^*$$
$d$ is the chordal distance and $\vec d$ is the chordal displacement. They will correspond to the euclidean distance between the two points on the sphere written in Cartesian coordinates.
Meanwhile: I have found a solution for Laplace's PDE on the spherical shell. I am now attempting to use a Green's function to reconstruct a solution to the particular equation.