Solution (source code)

= Solution

Put $\Psi=-\Phi=GM(r^p+a^p)^{-1/p}$. The spherical <Poisson equation> gives the <hypervirial model> density
$$
\boxed{
\rho(r)=\frac{(p+1)Ma^p}{4\pi}
\frac{r^{p-2}}{(r^p+a^p)^{2+1/p}}}
=\frac{(p+1)Ma^p}{4\pi(GM)^{2p+1}}
r^{p-2}\Psi^{2p+1}.
$$

For the proposed <hypervirial distribution function>, write $\mathcal E=-E=\Psi-v^2/2$ and $n=(3p+1)/2$. Direct velocity integration gives
$$
\begin{aligned}
\rho
&=A r^{p-2}
\int_0^{\sqrt{2\Psi}}v^p
(\Psi-v^2/2)^n\,dv
\int_0^{2\pi}d\varphi
\int_0^\pi\sin^{p-1}\alpha\,d\alpha\\
&=A\,2^{(p+1)/2}\pi^{3/2}
\frac{\Gamma(p/2)\Gamma(3(p+1)/2)}
{\Gamma(2p+2)}
r^{p-2}\Psi^{2p+1}.
\end{aligned}
$$
Matching this expression to the density fixes
$$
\boxed{
A=\frac{(p+1)a^p\Gamma(2p+2)}
{2^{(p+5)/2}\pi^{5/2}G^{2p+1}M^{2p}
\Gamma(p/2)\Gamma(3(p+1)/2)}}.
$$
This coefficient is positive for $0<p\leq2$, so the distribution function self-consistently generates the stated density and potential.