Solution (source code)

= Solution

For a spherically symmetric <number density>, the radial <Laplacian> is $r^{-2}\partial_r(r^2\partial_rn)$. Applying <separation of variables> and choosing the separation constant $-k^2$ gives
$$
\frac{R''+2R'/r}{R}=\frac{T'}T-\lambda=-k^2,
$$
so the <ordinary differential equations> are
$$
\boxed{R''+\frac2rR'+k^2R=0,\qquad T'=(\lambda-k^2)T.}
$$
The radial problem is an <eigenvalue problem> for $-\Delta$ with <Dirichlet boundary conditions>, hence $k^2>0$: multiplying the radial equation by $r^2\overline R$ and integrating gives $k^2\int_0^ar^2|R|^2dr=\int_0^ar^2|R'|^2dr>0$ for a nonzero regular mode vanishing at $a$.

To solve it explicitly, set $R=f/r$, corresponding to the suggested power $p=-1$. The equation reduces to $f''+k^2f=0$, so
$$
R(r)=\frac{A\sin(kr)+B\cos(kr)}r.
$$
Finiteness at the origin forces $B=0$, and the <Dirichlet boundary condition> at $r=a$ requires $k=m\pi/a$ for an integer $m\geq1$. The absence of interior zeros excludes $m\geq2$, selecting the fundamental <normal mode>. If $n_c>0$ denotes the initial central density, the result is
$$
\boxed{n(r,t)=n_c\frac{\sin(\pi r/a)}{\pi r/a}\exp\left[\left(\lambda-\frac{\pi^2}{a^2}\right)t\right],\qquad 0\leq r\leq a,}
$$
where the sine ratio is defined as $1$ at the origin. It is strictly positive in the interior for every finite time. The arbitrary positive scale $n_c$ cannot be fixed without additional initial data.

The spatial factor is time independent, so the density and the total <neutron> population have the same exponential growth rate. Thus the <spherical reaction-diffusion growth threshold> is
$$
\boxed{\lambda>\frac{\pi^2}{a^2}.}
$$
At equality the fundamental mode is stationary; below it the mode decays. For fixed $\lambda>0$, the equivalent critical-radius condition is $a>\pi/\sqrt\lambda$. This is the competition between multiplication in the bulk and <diffusion> loss at the surface in the stated dimensionless model.