A Laplace's equation describes a scalar field, say, density $\rho(r,z)$ distribution in a 2D axisymmetric coordinate system:
$$\nabla^2 \rho=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial \rho}{\partial r}\right)+\frac{\partial^2 \rho}{\partial z^2}=0,\quad (r,z) \in \Omega$$
where $r$ and $z$ are the radial and axial coordinates, $\Omega$ is a spherical domain with an elliptical source of species in its center.
If one imposes a zero density gradient in the direction of the unit normal, $\mathbf{n}=l^{-1}(r,z)$, of the outer boundary of $\Omega$, which is a bigger and bigger sphere surface $l=(r^2+z^2)^{1/2}$ at infinity with:
\begin{equation} \lim_{l \to +\infty}\frac{\partial \rho}{\partial \mathbf{n}} =\lim_{l \to +\infty}\frac{r\rho_r+z\rho_z}{l}=0, \quad \mbox{with} \quad l=(r^2+z^2)^{1/2}. \end{equation}
I could say that the density is constant in a local region at infinity based only on the condition of zero density gradient in normal at infinity and the fact that there is no density gradient tangential to the sphere surface.
My question is: Could I further deduce that the density $\rho$ is uniform everywhere in the domain $\Omega$ up to (and including) the border between $\Omega$ and the source? I came up with this guess according to the nature of Laplace's equation: within the domain, the value (here is density) is never above the max. on the boundary and never below the min. on the boundary. Can I make the above conclusion just according to the boundary value on one of the boundaries as in my case? Thanks for any comments.
Imposing the homogeneous Neumann boundary condition on one of two boundary components does not allow you to conclude that the function is constant. It all depends on what happens on the other boundary component ("elliptical source"). For example, in the spherical shell $r<|x|<R$ one can solve the Laplace equation with boundary conditions $$\frac{\partial u}{\partial n}|_{ |x|=R} = 0, \quad \frac{\partial u}{\partial n}|_{ |x|=r} = g(x) $$ for any function $g$ such that $\int_{|x|=r}g(x) = 0$. Indeed, the total flow across the boundary is zero for any such $g$, which makes these boundary conditions compatible with $\Delta u=0$.