🐛 fix(potential): finite gradient & hessian at r=0 for radial potentials

[nstarman]

gradient and hessian returned NaN at exactly r = 0 for essentially every
radial potential. The cause is that a radius built as sqrt(x² + y² + z²) has an
infinite derivative where its argument vanishes, so the chain rule produces
0 · inf = NaN for any potential written in terms of r — HernquistPotential,
PlummerPotential, NFWPotential, LogarithmicPotential, and the rest.

The fix

A single helper, galax.potential._src.utils.safe_sqrt, offsets the radicand by
the smallest normal float:

tiny = jnp.finfo(jnp.promote_types(q2.dtype, float)).tiny
return jnp.sqrt(q2 + tiny)

A thin safe_vector_norm wraps safe_sqrt(sum(square(x))) so call sites read as
a norm. r_spherical uses it, covering the twelve modules that build a spherical
radius; the ellipsoidal radii in TriaxialHernquistPotential,
Vogelsberger08TriaxialNFWPotential, and LeeSutoTriaxialNFWPotential call
safe_sqrt directly, since their quadratic form is not a plain norm.

The offset is ~1e-308, a no-op at any meaningful magnitude. Values agree with
jnp.linalg.vector_norm to within a rounding ulp — not bitwise: under
jax.jit, XLA fuses sum(square(x)) differently from vector_norm's scaled
algorithm, and ~11% of sampled elements differ by 1 ulp (≤3e-16 relative).
safe_vector_norm also overflows above |x| ~ 1e154, which vector_norm's
scaling avoids and no physical position reaches.

It has to be applied to the primal rather than only to the tangent. I tried the
usual jnp.where "double-where" and a custom_jvp that keeps the primal exact;
both make the gradient finite but silently zero every Hessian at the origin,
and the custom_jvp also reports zero force at the centre of a Kepler point
mass. Recovering the correct Hessian of a cored profile needs Φ'(r)/r
evaluated at a genuinely non-zero r.

Second cause found along the way

Making the radius differentiable exposed a log(1 + x) underflow: at
x ~ 1e-155, log(1 + x) flushes to exactly 0, so NFWPotential returned
-0 at the origin instead of -Gm/r_s. Switched to log1p in the NFW, Burkert,
Lee–Suto, and Vogelsberger profiles — a genuine precision fix at small r
independently of this PR.

Results at r = 0

	before	after
PlummerPotential, IsochronePotential	NaN	exact analytic gradient and Hessian (GM/b³)
HernquistPotential, TriaxialHernquistPotential, LogarithmicPotential, StoneOstriker15Potential, LM10Potential	NaN	gradient 0 — the net force at the centre of a spherically symmetric system, by symmetry
NFWPotential, BurkertPotential, MilkyWayPotential, MilkyWayPotential2022, gNFWPotential	NaN value	correct value -Gm/r_s

Cuspy profiles report a very large finite Hessian (~1e18) rather than NaN, which
reflects the genuine divergence of the tidal field at a 1/r cusp.

Known limitation, pinned by a test

NFWPotential and BurkertPotential still have a NaN gradient at the
origin, from a second and independent cause: their profile's r-derivative is
itself 0/0 there (d/dr[log1p(x)/x] = [x/(1+x) − log1p(x)]/x²). A
differentiable radius cannot help; those formulas need a small-x series. Their
values are now correct. test_profile_formula_singularity_is_a_separate_issue
pins this so that fixing the profiles trips the test.

Two genuinely singular potentials also change character: KeplerPotential and
JaffePotential diverge at r = 0, and their value there goes from -inf to a
large finite number. Kepler's gradient stays NaN. Happy to special-case these if
you'd rather they keep returning -inf.

Tests

tests/unit/potential/test_origin.py, covering safe_sqrt value-exactness and
derivative finiteness, the exact analytic Hessian for cored profiles, zero
gradient for cuspy ones, and the known limitation above. Full suite passes; the
two gNFWPotential doctests that asserted (nan, nan) at the origin were
updated to the now-correct values.

CI note

tests/functional/test_mockstreamgenerator.py::test_second_deriv fails on this
branch, but it also fails identically on the base commit d1a031bc — I ran
it in a separate worktree to check:

	computed	reference	max rel. diff
base d1a031bc	-978142	-978157	0.00098106
this branch	-978137	-978157	0.00101161

The committed pytest-arraydiff reference is stale by ~1e-3 on darwin/arm64.
This change moves the value by a further ~5e-6 relative: the 1-ulp norm
differences above, amplified through second-order autodiff of a chaotic
10,000-step integration. test_first_deriv, same mechanism and same
rtol=1e-7, passes. I have deliberately not regenerated the reference
file.

The mypy pre-commit hook also reports 50 pre-existing errors in 23 files, an
identical count with and without this change; none are in the files touched here.

Full suite otherwise: 6047 passed, 182 skipped, 14 xfailed.