A point source in a isotropic homogeneous sphere, continuously produces particles (like charges) and generates currents flowing in the sphere.
A point sink conversely absorbs particles.
When the source and the sink produces and absorbs particles at the same rate, a steady current field is then formed.
In the genral case, an equal amount of sources and sinks are distributed in the sphere, and form a steady current field $\vec J(\vec r)$.
From conservation of flow we know $$ \nabla\cdot\vec J(\vec r)=-\frac{\partial \rho}{\partial t} =F(\vec r),\quad 0\leq r\leq R, $$ where $\vec J(\vec r)$ is the vector field to be solved, $F(\vec r)$ a given scalar field, $R$ the sphere's radius. Because particles cannot run out of the sphere, the current must turn back before approaching the spherical surface. Therefore the boundary condition reads, $$ \left.{\vec J}\cdot{\vec n}\right|_{r=R}=0, $$ where $\vec n$ is the normal vector of the surface.
Is this problem well-defined? How to solve it?
A point source can be mathematically described as $I\delta(\vec r)$ since $$ \int\mathrm d^3\vec rI\delta(\vec r)=I=\int \mathrm d\vec S\cdot\vec J(\vec r)=\int\mathrm d^3\vec r\nabla\cdot\vec J(\vec r). $$ And a dipole (a point source and a point sink of the same current, infinitesimally close to each other) can be described as $$ I\delta(\vec r+\vec d/2)-I\delta(\vec r-\vec d/2)=I\vec d\cdot\nabla\delta(\vec r)=\vec p\cdot\nabla\delta(\vec r), $$ where $\vec d$ is the vector from the source to the sink, and $I\to\infty,\vec d\to0$ while $\vec p=I\vec d$ remains a finite constant value, called a dipole.
The field $F(\vec r)$ can be seen as a summation of distributed dipoles. For brevity, let $F(\vec r)=\vec p\cdot\nabla\delta(\vec r)$, so $$ \nabla\cdot\vec J(\vec r)=\vec p\cdot\nabla\delta(\vec r),\quad 0\leq r\leq R,\\ \left.\vec J(\vec r)\cdot\vec n\right|_{r=R}=0. $$
I tried Fourier transformatin and get $$ \mathrm e^{-\mathrm i\vec k\cdot\vec R}\vec J(\vec R)+\mathrm i\vec k \cdot\vec j(\vec k)=\mathrm i\vec k\cdot\vec p, $$ where $\vec j(\vec k)$ is the transformed form. I can't go any further.
I thought I should find the vector eigenfunction $\vec u$ of the operator $\nabla$, which means $$ \nabla\cdot\vec u(\vec r)=\vec k\cdot\vec u(\vec r), $$ where $\vec k$ is a constant vector, the eigenvalue. I need three vector eigenfunctions for a complete basis.