potential ${\phi_{\rm N}}$ is replaced by the field equation $$\nabla\cdot\left[\mu\left({\Vert\nabla\phi\Vert\over{a_0}}\right)
                 \nabla\phi\right] = 4\pi G \rho,
\label{eqMOND}$$ where ${a_0}\simeq 1.2 \times 10^{-10} {\rm m\, s^{-2}}$ is a characteristic acceleration, $\Vert ...\Vert$ is the standard Euclidean norm, $\phi$ is the MOND gravitational potential produced by the density distribution $\rho$, and in finite mass systems $\nabla\phi\to 0$ for $\Vert{{\bf x}}\Vert\to\infty$. The MOND gravitational field ${{\bf g}}$ experienced by a test particle is $${{\bf g}}=-\nabla\phi,
\label{eqgv}$$ and the function $\mu$ is such that $$\mu(y)\sim\cases{y&for $y\ll 1$,\cr 1&for $y\gg 1$;}
\label{eqmulim}$$ throughout this paper we use $$\mu (y)={y\over\sqrt{1+y^2}}.
\label{eqmu}$$ In the so-called ‘deep MOND regime’ (hereafter dMOND), describing low-acceleration systems ($\Vert\nabla\phi\Vert \ll{a_0}$), $\mu(y)=y$ and so equation (\[eqMOND\]) simplifies to $$\nabla\cdot\left({\Vert\nabla\phi\Vert}\nabla\phi\right) = 4\pi G {a_0}\rho.
\label{eqdMOND}$$ The source term in equation (\[eqMOND\]) can be eliminated by using equation (\[eqPoisson\]), giving $$\mu\left(\Vert\nabla\phi\Vert\over{a_0}\right)\nabla\phi=\nabla{\phi_{\rm N}}+{{\bf S}},
\label{eqcurl}$$ where ${{\bf S}}={\rm curl\,}{{\bf h}}$ is a solenoidal field dependent on $\rho$ and in general different from zero. When ${{\bf S}}=0$ equation (\[eqcurl\]) reduces to Milgrom’s (1983) formulation and can be solved explicitly. Such reduction is possible for configurations with spherical, cylindrical or planar symmetry, which are special cases of a more general family of stratifications (BM; Brada & Milgrom 1995). Though 