(Other) Rob's answer seems good to me, but let me offer another way of thinking.
As you approach the surface of the sphere very closely, the electric field should resemble more and more the electric field from an infinite plane of charge.
If you check Gauss's law (recalling that the field in the conductor is zero) you will see that if the surface charge density is $\sigma=Q/4\pi R^2$, then indeed the field at the surface is $\sigma/\epsilon_0$ as in the infinite charge of plane case.
Such a field is constant, the field lines are parallel and non-diverging, and the infinities associated with the field due to point charge do not arise.
Try doing the problem without the dielectric filling the sphere first as a warmup. Here's why that's useful:
Note that the definition of the displacement field is:$$\mathbf{D} = \varepsilon_0 \mathbf{E} + \mathbf{P},$$ where $\mathbf{P}$ is the dipole moment density. Taking the curl of both sides of the equation in the static state gives:
$$\nabla \times \mathbf{D} = \nabla \times \mathbf{P}.$$ Taking the divergence gives:
$$\nabla \cdot \mathbf{D} = \rho + \nabla \cdot \mathbf{P}.$$ We usually define the bound charge density as $\rho_{\mathrm{bound}} \equiv -\nabla \cdot \mathbf{P}$, and the remaining charge density is then called 'free', $\rho_{\mathrm{free}}\equiv \nabla \cdot \mathbf{D}$.
The divergence and curl equations above can be converted into the boundary conditions for the three vector fields across surfaces. Applying Gauss's law to the divergence equations gives:
$$\begin{align}\Delta \mathbf{P} \cdot \hat{n} & = \sigma_{\mathrm{bound}}, \\
\Delta \mathbf{D} \cdot \hat{n} & = \sigma_{\mathrm{free}},\ \mathrm{and} \\
\Delta \mathbf{E} \cdot \hat{n} & = \frac{\sigma_{\mathrm{bound}}}{\epsilon_0}, \end{align}$$ where $\hat{n}$ is a unit vector perpendicular to the surface and the $\Delta$s are understood to be a difference on one side of the surface from the other. These conditions are frequently expressed in terms of the components of the field perpendicular to the surface (eg $\mathbf{E}_\perp$), as implied by the dot product. In words, these mean that the change perpendicular components of the fields are dictated by their associated surface charge density.
Unlike the electric field, though, the displacement and polarization can have solenoidal components, as you can see from the curl equation above. Applying Stoke's theorem at the boundary surface for them gives:
$$\begin{align}\Delta \mathbf{D}\times \hat{n} &= \Delta \mathbf{P}\times \hat{n}, \ \mathrm{and} \\
\Delta \mathbf{E} \times \hat{n} &= 0.\end{align}$$
Like the conditions on the perpendicular components, this can be expressed as a condition on the components of the fields parallel to the surface (eg $\mathbf{E}_{||}$). In words, the parallel components of the electric field must be continuous, and the discontinuity in the parallel components of the polarization must be balanced by a discontinuity in the parallel components of the displacement field.
Note that I don't think you actually need to find $\mathbf{D}$ to solve this problem - it should be sufficient to solve the dielectricless problem, use that to infer the polarization in the sphere, and from that calculate how the net surface charge density is modified by the bound charge.
All of that said, you should find that the empty sphere is a pure dipole field on the inside and the outside. Keep in mind that a dipole field appropriate for $r\rightarrow 0$, ie 'inside', has a different structure than one for $r\rightarrow \infty$, 'outside'. Your textbook should have a section on solving such problems using Legendre polynomials, as mentioned in a comment by @PrasadMani. The multipole fields are separated by the fact that they behave under rotations in orthogonal fashions, and it shouldn't surprise you to find that adding a rotationally invariant dielectric sphere won't cause a mixing of multipole field types.
Best Answer
The total charge density, which is the sum of free charge density and bound charge density is constant all over the sphere, but the free and bound charge densities differ for upper and lower half of the sphere.
Your question about symmetry of total charge density can be answered easily, assuming that you know the symmetry of electric field. (your reasoning for it's spherical symmetry is correct):
($E_1$ means electric field in the first region and $E_2$ the field of the second region)
$$\mathbf{E}_1=\mathbf{E}_2 =\mathbf{E}_{out} \tag{as you stated}$$ $$\mathbf{E}_{out} . \hat r-\mathbf{E}_{in}.\hat r=\frac{\sigma}{\epsilon}_0 \,\,\,,\,\,\,\mathbf{E}_{in}=0 \to \sigma={\epsilon}_0 \mathbf{E}_{out} . \hat r$$
You can arrive at this result by calculating free and bound charge densities directly by solving for potential using Laplace equation too:
The problem is clearly azimuthal symmetric, and also according to your argument, it isn't a function of $\theta$ neither; i.e., the problem is 1-D and $V$ is a function of $r$.
Using the fact that potential is continuous at the boundary of the two regions and also vanishes at $\infty$, the only remaining term of the potential expansion in spherical coordinates will be:
$$V=\frac{Cq}{r}$$ $$ \mathbf{D}=-\epsilon \nabla V \to \cases{\mathbf{D}_1=\frac{\epsilon_1Cq}{r^2}\hat r \\ \mathbf{D}_2=\frac{\epsilon_2Cq}{r^2}\hat r}$$
Now, assuming a spherical surface surrounding the sphere and applying Gauss law for $\mathbf{D}$, we will find $C$:
$$\oint{\mathbf{D}.d\mathbf{S}}=Q_{free}\to C=\frac{1}{2\pi (\epsilon_1 + \epsilon_2)}$$
$$V=\frac{1}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r}$$
Now we find $\mathbf{D}_1$ , $\mathbf{D}_2$ and $\mathbf{E}$ :
$$\mathbf{E}=\frac{1}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r^2}\hat r$$
$$\mathbf{D}_1=\frac{\epsilon_1}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r^2}\hat r$$
$$\mathbf{D}_2=\frac{\epsilon_2}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r^2}\hat r$$
using the relation $\sigma_f=\mathbf{D}.\hat r $ we will have:
$$\to \cases{\sigma_{1f}=\frac{q\epsilon_1}{2\pi (\epsilon_1 + \epsilon_2)} \frac{1}{R^2} \\ \sigma_{2f}=\frac{q\epsilon_2}{2\pi (\epsilon_1 + \epsilon_2)} \frac{1}{R^2}}$$
From these two relations for $\sigma_{f}$s it is evident that the total free charge is equal to $q$, as expected: $$\sigma_{1f}\times 2\pi R^2+\sigma_{2f}\times 2\pi R^2=q$$
Now we find bound charge densities. First we should find polarization vectors $\mathbf{P}_1$ and $\mathbf{P}_2$:
$$\mathbf{P}=\mathbf{D}-\epsilon_0 \mathbf{E} \to \cases{\mathbf{P}_1= \frac{\epsilon_1- \epsilon_0}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r^2}\hat r \\ \mathbf{P}_2=\frac{\epsilon_2- \epsilon_0}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{r^2}\hat r }$$ and
$$\sigma_b =\mathbf{P}.\hat r\to \cases{\sigma_{b1}=-\frac{\epsilon_1- \epsilon_0}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{R^2} \\ \sigma_{b2}=-\frac{\epsilon_2- \epsilon_0}{2\pi (\epsilon_1 + \epsilon_2)} \frac{q}{R^2} }$$
Now you can see that on each half sphere the sum of bound and free charge densities is equal, as expected: $$\sigma_{b1}+\sigma_{f1}=\sigma_{b2}+\sigma_{f2}$$