How to determine charge density using Dirac deltas in advance? --- not after the fact Context
I have already asked one question regarding charge densities, Diracs, and Heavisides [0]. At the time of writing, that question remains open. More importantly, I still remain unclear regarding how to write charge densities with Diracs, and Heavisides. Here is an example to show that I lack clarity.
This question here revolves around how to define a surface charge. There is a large set up to the problem. Once set up, then comes the discussion of surface charge. I want to solve for the electric potential $V$ inside and outside a conductive ball of radius $R$ with a total charge $Q$. The problem can be solved with Gauss' Law, and it can be solved with spherical harmonics. I will do both.
Gauss' Law
Here, I use Gauss' Law, to find the electric field; and then I use integration and continuity of $V$ at the boundary to solve for $V$. I the potential at infinity as the zero referene.
\begin{align*}
\mathbf{E} {\left(\mathbf{r}\right)}
&= 
\begin{cases}
\mathbf{0},~&\text{for} 
~r< R
\\\\
  \frac{Q}{4\,\pi\,\epsilon_o \,r^2  }\,\mathbf{\hat{r}} ,~&\text{for} 
~r> R. 
\end{cases}
\end{align*}
\begin{align*}
V {\left(\mathbf{r}\right)}
&= 
\begin{cases}
  \frac{Q}{4\,\pi\,\epsilon_o \,R  }   
,~&\text{for} 
~r\leq R
\\\\
  \frac{Q}{4\,\pi\,\epsilon_o \,r  } ,~&\text{for} 
~r> R. 
\end{cases}
\end{align*}
Expansion in Spherical Harmonics
In a nut shell, since this problem has azimuthal symmetry, the solution can be written [1] in terms of Legendre polynomials of degree $\ell$, $P_\ell$, as
$$ V(r,\theta,\phi)
=
\begin{cases}
\sum_{\ell=0}^\infty A_\ell\, r^\ell\,P_\ell(\cos\theta), &(r\leq R)
\\
\sum_{\ell=0}^\infty \frac{A_\ell\,R^{2\,\ell+1}}{r^{\ell+1}} \,P_\ell(\cos\theta), &(r\geq R)
;
\end{cases} $$
where
$$  A_\ell = \frac{1}{2\,\varepsilon_o\,R^{\ell-1} }\,\int_0^\pi \sigma(\theta)  P_\ell(\cos\theta)\,\sin\theta\,d\theta . $$
Question
(1) What is $\sigma(\theta) $
(2) What am I misunderstanding about how to use Dirac delta distributions and Heaviside step function as I attempt to write charge densities in terms of the same?
Answer to Question 1
One the one hand, I argue, but not with self confidence and only because I know the answer in advance, that surface charge density is invariant with respect to $\theta$ (i.e., $\sigma(\theta) = \sigma_o$). In such case, according to the orthogonality of the Legendre polynomials [2]
\begin{align}
 A_\ell 
&= 
\frac{\sigma_o}{2\,\varepsilon_o\,R^{\ell-1} }\,\int_0^\pi P_\ell(\cos\theta)\,\sin\theta\,d\theta . 
\\
&= 
\begin{cases}
0 & (\ell\neq 0)
\\
\frac{\sigma_o\,R }{ \varepsilon_o }.
\end{cases} 
\end{align}
Thus,
$$ V(r,\theta,\phi)
=
\begin{cases}
\frac{\sigma_o\,R }{ \varepsilon_o } , &(r\leq R)
\\
\frac{\sigma_o\,R }{ \varepsilon_o }\,\frac{ R }{r }, &(r\geq R)
;
\end{cases} $$
It appears that for the answer here to be self-consistent with the answer from Gauss' law, that
$$ \sigma(\theta) = \frac{Q}{4\,\pi\,R^2}.\tag{1}$$
This seems like a sensible description of the surface charge. Ulitmately,
$$ V(r,\theta,\phi)
=
\begin{cases}
\frac{ Q   }{ 4\,\pi\,R \,\varepsilon_o } , &(r\leq R)
\\
\frac{ Q   }{ 4\,\pi \,\varepsilon_o\,r }, &(r\geq R)
.
\end{cases} $$
Well, that appears to have went well.
One the other hand, I argue that surface charge density varies with respect to $\theta$. In such case, according to the orthogonality of the
I write that
\begin{align}
Q 
&=
\int_{\mathbb{R}^3} \rho(r,\theta,\phi)\,d\tau
\\
&=
\int_{\phi=0}^{2\,\pi}\,\int_{\theta=0}^{\pi}\,\int_{r=0}^\infty \rho(r,\theta,\phi)\,r^2\,\sin\theta\,dr\,d\theta\,d\phi.
\end{align}
Since the charge only exists at $r=R$, I write $\rho(r,\theta,\phi) = \sigma(r,\theta,\phi)\,\delta(r-R)$, then
\begin{align}
Q 
&=
\int_{\phi=0}^{2\,\pi}\,\int_{\theta=0}^{\pi}\,\int_{r=0}^\infty \sigma( \theta,\phi)\,\delta(r-R)\,r^2\,\sin\theta\,dr\,d\theta\,d\phi.
\\
&=
\int_{\phi=0}^{2\,\pi}\,\int_{\theta=0}^{\pi} \sigma(\theta,\phi) \, R^2\,\sin\theta\,dr\,d\theta\,d\phi .
\end{align}
I argue, but not self-convincingly, that since the differential area element, $da$ is $da =R^2\,\sin\theta\,d\theta\,d\phi$ and since the surface charge density should be constant with respect to $da$, that $ \sigma( \theta,\phi) =\sigma_o$. Then,
\begin{align}
Q 
&=
\int_{\phi=0}^{2\,\pi}\,\int_{\theta=0}^{\pi} \sigma_o \, R^2\,\sin\theta\,dr\,d\theta\,d\phi .
\\
&=
4\,\pi \sigma_o \, R^2.
\end{align}
Again, in oder to be consistent, it seems that
$$\sigma = \frac{Q}{4\,\pi\,R^2}. \tag{2}$$
In such case, I will get the same answer as before for $V(r,\theta,\phi)$.
My true question
Is there some way to more formal method to determine how to use Dirac delta distributions and Heavside step functions? Both ways that I performed here are round about. In Equation 1 and Equation 2, I determine the charge density after the fact. I do not like that, as I believe that I will fail to properly apply Dirac delta functions and Heaviside step functions on more complicated surfaces---when I can't use an after-the-fact argument.
Bibliography
[0] Dirac delta, Heaviside step, and volume charge density
[1] Griffiths, "Introdution to Electrodynamics" 2nd Edition p. 142-143.
[2] Wikipedia contributors. "Legendre polynomials." Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 21 Feb. 2021. Web. 8 Mar. 2021.
 A: 
One the one hand, I argue, but not with self confidence and only because I know the answer in advance, that surface charge density is invariant with respect to $\theta$.

Your solution using Gauss' law implicitly assumes spherical symmetry.  If $\sigma(\theta)$ varies with $\theta$, then $\mathbf E \neq \frac{Q}{4\pi \epsilon_0 r^2} \hat r$.
In any case, if your surface charge density on the surface of the sphere is $\sigma(\theta)$, then the spatial charge distribution (in 3D) is given by $\rho(\mathbf r) = \sigma(\theta) \delta(r-R) = \frac{Q}{4\pi R^2} \delta (r-R)$. One can see immediately that this works, because
$$Q_{inside} = \int _0^{r_0} r^2 \mathrm dr \int_0^{2\pi} \mathrm d\phi \int_{-1}^1 \mathrm d(\cos(\theta)) \ \left[\frac{Q}{4\pi R^2} \delta(r-R)\right] = \begin{cases}0 & r_0<R\\Q & r_0 > R\end{cases}$$
On a more complicated surface, perhaps defined by $R=R(\theta, \phi)$, then if the surface charge distribution is given by $\sigma(\theta,\phi)$, the volume charge distribution will be given by $\rho(r,\theta,\phi) = \sigma(\theta,\phi) \delta\big(r - R(\theta,\phi)\big)$.
A: The question here asks about doing this in advance. To determine the matter in advance requires knowledge of the scale factors of the various curvilinear coordinates. For example, the spherical coordinate system    has the scale factors $1$,  $r$, and $r\sin\theta$ for coordinate $r$, $\theta$, and $\varphi$, respectively.
I am not going to do an exotic example here. But let's take a look at a problem with spherical symmetry. Let's say a half-dome shell of radius $R$ in the upper $z$ plane. The goal is determine
$$
\rho = \rho(r,\theta, \varphi).$$
For each coordinate we have two choices: either a Dirac Delta  or a Heaviside step. So,  in $\rho $ we have to have a factor that is exclusively
$$ \frac{\delta(r-r_o)}{1}\quad \text{or}\quad \left[H(r-r_i )  -  H(r_f-r ) \right],$$
a factor that is exclusively
$$ \frac{\delta(\theta-\theta_o)}{r}\quad \text{or}\quad \left[H(\theta-\theta_i )-  H(\theta_f-\theta ) \right],$$
and a factor that is exclusively
$$ \frac{\delta(\varphi-\varphi_o)}{r\,\sin\theta}\quad \text{or}\quad \left[\frac{\delta(\varphi-\varphi_i)}{1} -  H(\varphi_f-\varphi ) \right].$$
In the case of a half-dome shell in the positive $z$ plane, we then have
$$ 
\rho = \rho(r,\theta,\varphi) = \rho_o\, 
\left[
\frac{\delta(r-R)}{1}
\right]
\left[
H(0-\theta ) 
 -  H{\left(\frac{\pi}{2}-\theta \right)} 
\right]
\left[
\frac{\delta(\varphi-0)}{1} -  H(2\,\pi-\varphi ) 
\right]\,
 $$
In this example $$\rho_o = \frac{Q}{2\,\pi\,R^2}$$
Let's put it all together.
\begin{align*}
Q 
&= 
\int_{\mathbf{R}^3} \rho(r^\prime ,\theta^\prime,\varphi^\prime)
 dV^\prime 
\\
&= 
\int_{0}^\infty 
\int_{0}^{\pi} 
\int_{0}^{2\,\pi\infty} 
\frac{Q}{2\,\pi\,R^2}  \,\left[ \frac{\delta(r- R )}{1}  \right]  \,
\left[ 
H(\theta - 0  )
-
H{\left(\frac{\pi}{2} - \theta \right)}
\right] 
\left[  
H(\phi - 0  )
-
H(2\,\pi - \phi   )
   \right]{r^\prime}^2\,\sin\theta^\prime  \,dr^\prime
 \,d\theta^\prime \,d\varphi^\prime
\\
&= 
\int_{0}^{\infty} 
\frac{Q}{2\,\pi\,R^2}  \,\left[ \frac{\delta(r^\prime- R )}{1}  \right]  \,
 {r^\prime}^2\, dr^\prime
\\
&\times 
\int_{0}^{\pi} 
\left[ 
H(\theta^\prime  - 0  )
-
H{\left(\frac{\pi}{2} - \theta^\prime  \right)}
\right] 
 \sin\theta^\prime    \,d\theta^\prime  
\\
&\times 
\int_{0}^{2\,\pi } 
\left[  
H(\phi^\prime - 0  )
-
H(2\,\pi - \phi^\prime   )
   \right]  \,d\varphi^\prime
\\
&= 
\frac{Q}{2\,\pi\,R^2} 
\left(R^2\right)
\left(1\right)
\left(2\,\pi\right)
\end{align*}
Thus, we find the consistent result, which is that
$$Q
 =   Q    \,.   $$
