Issues with Sampling from Power-Law Potential DFs

[slizewsk]

🐛 Bug: Sampling Fails for PowerSphericalPotential

Attempting to sample from isotropic and constant anisotropy distribution functions (eddingtondf, constantbetadf) based on PowerSphericalPotential fails with a divide by zero error.

This issue appears both in the self-consistent case and when using a separate tracer density. While the potentials do diverge at very small r, the class should be designed to handle this appropriately.

Reproducible Example

For the self-consistent case:

from galpy.potential import PowerSphericalPotential
from galpy.df import eddingtondf

gamma = 0.5
sph_pot = PowerSphericalPotential(amp=1, alpha=gamma + 2)
iso_df = eddingtondf(pot=sph_pot, rmax=1e4)

samps = iso_df.sample(n=3000, return_orbit=True)  # Fails here

Full traceback: https://pastebin.com/5Z32qx0u

Expected Behavior

• The DF should allow bound orbit sampling for both isotropic and anisotropic forms.
• This issue does not occur for isotropicHernquistdf and constantbetaHernquistdf. See https://pastebin.com/XcfnDTtt for a working example using the Hernquist model to show the

System Details

• macOS-15.3.1-arm64-arm-64bit
• Python 3.11.8 | packaged by conda-forge | (main, Feb 16 2024, 20:49:36) [Clang 16.0.6 ]
• numpy 1.26.4
• scipy 1.13.1
• matplotlib 3.8.4
• galpy 1.10.2
• astropy 6.1.0

Additional Context

• I am testing PowerSphericalPotential for a Bayesian model fitting package using CmdStanPy.
• The potential is defined as: $\Phi(r) = -\frac{\Phi_0}{r^\gamma}$
• With a self-consistent density: $\rho(r) = \frac{\gamma (1-\gamma) \Phi_0}{4\pi G} r^{-\gamma-2}$
• In galpy, PowerSphericalPotential(alpha=γ+2) should correspond to this density.
• Sampling also fails for the KeplerPotential with the same error.
• For the NFWPotential, sampling works for the isotropic DF but fails for the constant-anisotropy DF. However, this produces a different error.
• I noticed the eddingtondf ._Emin returns -inf.. perhaps related?

Let me know if you'd like more details or a different test case.