To be honest, I just learned about all this myself in the last months, so I am not sure whether this is actually correct.
Since you have this spherical symmetry, I think that you need spherical harmonics. They are orthogonal functions, think of them as a Fourier series on the surface of a sphere.
Your charge density $\sigma$ does not depend on $\phi$, therefore we can use the simpler Legendre-polynomials, where $m = 0$. Fist, we want to express the potential like so:
$$
\varphi(\vec x) = \frac1{\varepsilon_0} \sum_{l=0}^\infty \frac{1}{2l+1} \varrho_{l} \frac{1}{r^{l+1}} Y_{l,0} (\theta, \phi)
$$
The coefficients $\varrho_l$ are given by:
$$
\varrho_l = \int \mathrm d^3 x'\, r'^l \varrho(\vec x) Y^*_{l,0} (\theta', \phi')
$$
Now with
$$
Y^*_{l,0}(\theta', \phi') = \sqrt{\frac{2l+1}{4\pi}} P_l(\cos \theta')
$$
and
$$
P_0(x) = 1
,\quad
P_1(x) = x
$$
we can calculate the coefficitents.
But first, we need to convert the surface charge density $\sigma$ into a volume charge density. For that, we use the $\delta$-distribution:
$$
\varrho(\vec x) = \delta(r-R) \sigma_0 \cos(\theta)
$$
If you plug those into the $\varrho_l$, you will get $\varrho_0 = 0$ and $\varrho_1 = \sqrt{\frac{3}{4 \pi}} \frac{4}{3} \pi R^3 \sigma_0$. I hope this is correct.
Then we can put this into the first formula and get $\varphi$:
$$
\varphi(\vec x) = \frac{1}{\varepsilon_0} \sqrt{\frac{3}{4 \pi}} \frac{4}{3} \pi R^3 \frac{1}{r^2} \sigma_0 \cos \theta
$$
Since this is a pure dipol potential, the $1/r^2$ seems about right. And if you look at the dimensions, the $R^3 \sigma_0 / r^2$ have just the needed Charge/Length.
Spherical harmonics might be overkill, maybe there is a simpler method to do this.
In order to find the potential function, one usually employs the methods used here. Then, once you have this function for the potential:
$V(r)=\frac{2q}{4\pi\epsilon_0 a^2}(r-\frac{R^3}{r^2})cos(\theta)$
Find the electric field at the surface of the grounded sphere of radius R
$E(R)=\frac{\partial{V}}{\partial r}\bigg|_{r=R}=\frac{2q*cos(\theta)}{4\pi\epsilon_0 a^2}*3=\frac{6q*cos(\theta)}{4\pi\epsilon_0 a^2}$
$\frac{\sigma}{\epsilon_0}=E(R)
$
$\Rightarrow \sigma=\frac{6q*cos(\theta)}{4\pi a^2}$
That should answer your first question.
As for your second question, no, the potential inside each sphere will not change. Why?
Let's see what happens when a charged sphere is brought near another charged conducting sphere, or really any charge configuration changes the electric field in the environment of the charged conductiong sphere in question. The elctric field everywhere will change (except, of course inside the sphere) and there are lots of nifty tricks to figure out the new electric field, the method of image charges being the most prominent. Now, consider a regular conducting sphere with a charge q and a radius R sitting somewhere with no external fields. The potential at its surface and inside it $=\frac{q}{4\pi\epsilon_0R}$, a constant. That means that you can add this constant to the (new) potential function for the region outside without changing the electric field outside, because the gradient of a constant is $0$. Now, inside the sphere, $E=0$, which means the potential inside is a constant function (gradient of a constant is zero). You can set that constant to be whatever you want and $=\frac{q}{4\pi\epsilon_0R}$ looks pretty good, unless you want a discontinuity in the potential function.
Best Answer
This is a bit of an old question but I stumbled upon it and I believe it a great learning experience and quite a beautiful question.
What you got right (most of it)
I'll begin by pointing out that you are entirely right! your configuration creates a perfectly uniform field in the ball!
(I'll use ball to refer to the volume, while sphere would be used for the thin shell. I highly recommend using those definitions to avoid confusion)
As you pointed out the uniqueness theorem doesn't mean that your solution is unique. In fact there are infinite configurations of charges that would create uniform field. Unfortunately your solution is slightly problematic, because you described your charges as an nonuniform distribution of charge inside the volume of a perfect conductor.
Why is that a problem? Because inside of a perfectly conducting body there are no electric fields! (Any field the body could be subjected to would create surface currents that would immediately cancel the fields for the volume of the body.) Without electric fields, why would the charge ever arrange itself so nicely? As there are no way to create charge distribution inside perfect conductors we generally treat those solutions as "Non-physical", meaning mathematically correct but not real. Surface charges are on the other hand completely possible!
I will note, that your solution is on the right track. By taking the right limit and changing your notation a bit you can get the charge distribution as a surface charge.
So, now that we understand the problem, let's answer each of your questions
Is there a unique solution?
Or as you placed it "Is there some type of uniqueness theorem to determine this situation?" I'm going to take some liberties and assume the conductive sphere is neutrally charged, otherwise there are infinite solutions. (differing only by the charge of the sphere) I'll also assume all the charge is on the sphere and not on the ball. With those assumptions the solution is unique! Although I'm not aware of any named theorem to the problem, there is an easy way to see that there is only a single solution that has only surface charges.
We'll prove it thusly:
($u$ and $v$ are parameters that describe our surface. With a sphere we could choose for example $\theta$ and $\phi$. $A$ is an arbitrary constant and $\hat z$ is an arbitrary direction)
You can get the solution immediately here by making an argument about the energy of the new electric configuration, which must be zero. We'll make the same argument in the proof, albeit a bit more formally.
and so we proved that $\sigma$ is unique.
What is the right solution?
You are already on the right track! If we take the limit of $lim{a\rightarrow 0^+}$ and increase $rho$ accordingly so to keep $\rho \cdot a$ constant (meaning our polarity is the same as before) we can immediately get the surface charges.
Let's look at an arbitrary slice of the ball, at height H from the top, we'll consider the top half for now. (Up is $\hat a$ for simplicity) The radius of this slice is $\sqrt{R^2 - H^2}$ the area of this slice is therefore the square of that times $\pi$ meaning $A(H)= \pi \cdot (R^2 - H^2) $.
What is the total charge at a single slice? $Q (H) = \rho (A(H)-A(H+a))$. All the charge is distributed over the thin shell, since it's simple to prove that the charge is zero anywhere else. We immediately get that $\sigma (H) = \rho (A(H)-A(H+a))/{2 \cdot \pi \cdot \sqrt{R^2 - H^2}}$
$\sigma (H) = \rho (2H \cdot a)/{2 \cdot \pi \cdot \sqrt{R^2 - H^2}}$ We can immediately find $(\rho \cdot a)$ since we know the polarity of the final charge distribution. good luck