Solution (source code)

= Solution

Put $V=|v_0|>0$, take $q>0$ without loss of generality, and write $s=R^2+z^2/q^2$. The <Newtonian gravitational potential> is singular at the origin; all pointwise claims below concern $s>0$. The zero of the <Newtonian gravitational potential> is the one in the paper, with an implicit fixed length unit inside the <logarithm>.

The <Poisson equation for Newtonian gravity> in <cylindrical coordinates> gives
$$
4\pi G\rho=\frac1R\partial_R(R\Phi_R)+\Phi_{zz},\qquad \Phi_R=\frac{V^2R}{s},\qquad \Phi_z=\frac{V^2z}{q^2s}.
$$
The two contributions to the <Laplacian> are
$$
\frac1R\partial_R(R\Phi_R)=\frac{2V^2}{s}-\frac{2V^2R^2}{s^2},\qquad \Phi_{zz}=\frac{V^2}{q^2s}-\frac{2V^2z^2}{q^4s^2}.
$$
Consequently the <mass density> of the <axisymmetric logarithmic gravitational potential> is
$$
\boxed{\rho(R,z)=\frac{V^2}{4\pi Gq^2}\frac{R^2+(2-q^{-2})z^2}{(R^2+q^{-2}z^2)^2}.}
$$
The coefficient of $R^2$ is positive. The coefficient of $z^2$ determines the <density positivity for an axisymmetric logarithmic potential>: \b[strict positivity at every noncentral point requires $q^2>1/2$ and $v_0\ne0$]. Nonnegativity permits $q^2=1/2$, where the <mass density> vanishes on the noncentral symmetry axis. For $q^2<1/2$ it is negative near that axis. Either sign of a nonzero real $q$ gives the same model, and $q=0$ is undefined. Setting $v_0=0$ yields vacuum away from the singular origin, not strictly positive <mass density>. The $r^{-2}$ central cusp is locally integrable and contains no hidden point <mass>: the <Newtonian gravitational field> flux through a shrinking sphere is of order its radius.

An equatorial <circular orbit> requires inward <gravitational acceleration> equal to $v_c^2/R$. Thus the <circular speed> is
$$
v_c^2=R\Phi_R(R,0)=V^2,\qquad \boxed{v_c(R)=V\quad(R>0).}
$$
This is a <flat galaxy rotation curve>. The singular center does not define an additional circular orbit at $R=0$.

\Image[/past-exam-of-the-mathematics-course-of-the-university-of-cambridge/2012/iii/paper-59-rotation-curve.png]
{title=Flat equatorial rotation curve of the singular logarithmic galaxy; the central point is excluded}

The <specific orbital energy> $E=\tfrac12(v_R^2+v_\phi^2+v_z^2)+\Phi$ is conserved because the <Newtonian gravitational potential> is time independent. The axial <specific angular momentum> $L_z=Rv_\phi$ is conserved because the <Newtonian gravitational potential> is axisymmetric. Hence the <Jeans theorem>, or directly $dF(E,L_z^2)/dt=0$ in the <Collisionless Boltzmann equation>, permits a stationary <two-integral galactic distribution function>. Using $L_z^2$ makes this <galactic distribution function> even under reversal of the azimuthal <velocity>, so it describes a nonstreaming population. A general axisymmetric population can have an odd part in $L_z$, or dependence on a third <integral of motion>; neither is required for this construction.

At each position define $\rho\langle g\rangle=\int gF\,d^3v$. The <two-integral galactic distribution function> is invariant under swapping $v_R$ and $v_z$, since only their sum of squares occurs in $E$. It is also even in each <velocity> component. Swapping integration variables proves <meridional velocity isotropy of a two-integral distribution>, while reflecting one variable makes each mixed integrand odd. Therefore
$$
\boxed{\langle v_R^2\rangle=\langle v_z^2\rangle,\qquad \langle v_Rv_z\rangle=\langle v_Rv_\phi\rangle=\langle v_\phi v_z\rangle=0.}
$$
All mean <velocities> vanish, so these are also the corresponding statements for the <velocity ellipsoid> and its centered <covariance> tensor. They do not require $\langle v_\phi^2\rangle$ to equal either meridional second moment.

To normalize the proposed <logarithmic-potential two-integral distribution>, first separate the <mass density> into terms with the same spatial factors as the two exponentials:
$$
\rho=\frac{V^2}{4\pi Gq^2}\left[\frac{2q^2-1}{s}+\frac{2(1-q^2)R^2}{s^2}\right].
$$
Since $e^{-2\Phi/V^2}=s^{-1}$ and $e^{-4\Phi/V^2}=s^{-2}$, the <Gaussian integral> and its derivative give
$$
\int e^{-v^2/V^2}\,d^3v=\pi^{3/2}V^3,\qquad \int v_\phi^2e^{-2v^2/V^2}\,d^3v=\frac{\pi^{3/2}V^5}{8\sqrt2}.
$$
Thus velocity integration of the <galactic distribution function> yields $AR^2s^{-2}\pi^{3/2}V^5/(8\sqrt2)+Bs^{-1}\pi^{3/2}V^3$. Matching the independent spatial factors determines
$$
\boxed{A=\frac{4\sqrt2(1-q^2)}{\pi^{5/2}Gq^2V^3},\qquad B=\frac{2q^2-1}{4\pi^{5/2}Gq^2V}.}
$$
This verifies both the stationary <Collisionless Boltzmann equation> and the required <Poisson equation for Newtonian gravity>. A shift $\Phi\mapsto\Phi+C$ changes the normalizations to $Ae^{4C/V^2}$ and $Be^{2C/V^2}$; the physical <mass density> is unchanged.

There is an important distinction between positive <mass density> and a nonnegative <galactic distribution function>. For $1/2\le q^2\le1$, both coefficients above are nonnegative. For a prolate model $q^2>1$, $A<0$ and one must test the sum rather than reject the model just because one term is negative. On accessible phase space,
$$
L_z^2e^{-2E/V^2}=\frac{R^2}{s}v_\phi^2e^{-v^2/V^2}\le\frac{V^2}{e},
$$
with equality in the equatorial plane at $v_R=v_z=0$, $v_\phi^2=V^2$. The <nonnegativity bound for a prolate logarithmic-potential distribution> is consequently
$$
\boxed{q^2\le\frac{16\sqrt2-e}{16\sqrt2-2e}\quad(q^2>1).}
$$
For larger $q$, the displayed algebraic expression still integrates to the positive <mass density>, but is negative at those equatorial <circular orbit> phase points and is not a physical <galactic distribution function>. This extra restriction is distinct from the earlier density-only answer.