fix: Apparent unphysical behaviour in the cNFWMCRLudlowSph class

[mwiet]

Description

The cNFWMCRLudlowSph returns zero convergence_2d_from and an over-amplified deflections_yx_2d_from for  dimensionless core fractions f_c greather than 0.1 The convergence is also negative in this regime.

The cored-NFW family in autogalaxy.profiles.mass.dark.cnfw / cnfw_mcr suffers three related numerical issues at the core fractions (f_c = r_core / r_s ≈ 0.1 – 0.3) most relevant to self-interacting dark matter (SIDM) lensing applications:

1. convergence_2d_from returns identically zero everywhere on a uniform image-plane grid, for any combination of mass_at_200, redshift, and f_c we tested.
2. deflections_yx_2d_from scales non-monotonically with f_c: the deflection magnitude blows up by roughly an order of magnitude near f_c ≈ 0.2, and the implied kappa_s becomes negative for f_c ≳ 0.5. Compared to an equivalent NFWTruncatedSph (same mass_at_200, autolens-shipped MCR concentration), cNFWMCRLudlowSph with f_c = 0.2 produces deflections ~20× larger at radii ~r_s.
3. The deflection vector field flips sign (points INWARD insteadof outward) for f_c ≳ 0.22-0.30, well before kappa_s becomes negative. A cored halo in this regime acts as a negative-mass perturber: it deflects light toward its centre. The previous draft of this report only checked deflection magnitudes via np.linalg.norm and missed this — the inward-flip happens silently within the "positive kappa_s" range, so any user checking only kappa_s > 0 to validate a cored halo is unaware of the qualitative rendering error.

Steps to Reproduce 1

import numpy as np
import autolens as al
import autogalaxy as ag

mass = 1.0e9
# Four probe points 1 arcsec from a halo at the origin, on the +/-x/y axes.
probes = [(0.0, 1.0), (0.0, -1.0), (1.0, 0.0), (-1.0, 0.0)]
probe_grid = al.Grid2DIrregular(values=probes)

for f_c in (0.01, 0.10, 0.20, 0.30, 0.50):
    cnfw = ag.profiles.mass.dark.cnfw_mcr.cNFWMCRLudlowSph(
        centre=(0.0, 0.0), mass_at_200=mass, f_c=f_c,
        redshift_object=0.5, redshift_source=1.0,
    )
    defl = np.asarray(cnfw.deflections_yx_2d_from(grid=probe_grid))
    print(f"f_c={f_c}  kappa_s={cnfw.kappa_s:+.3e}")
    for (y, x), (ay, ax) in zip(probes, defl):
        outward = (ax * x >= 0) and (ay * y >= 0)
        print(f"  probe ({y:+.1f}, {x:+.1f}):  α=({ay:+.3e}, {ax:+.3e})  "
              f"{'OUTWARD (OK)' if outward else 'INWARD (sign-flipped)'}")

Expected Behaviour 1

For f_c ∈ (0, 1):

• convergence_2d_from should produce a finite, positive, smoothly varying field that integrates to the projected mass within the evaluation aperture.
• deflections_yx_2d_from magnitudes should vary continuously with f_c and approach the NFW value as f_c → 0. Beyond f_c ≈ 0.1 they currently jump by 10× before falling back to NFW-comparable values at f_c ≈ 0.5, which is not physically plausible.
• kappa_s should remain positive across the full f_c ∈ (0, 1) range.

Actual Behaviour 1

f_c=0.01  kappa_s=+8.132e-03   OUTWARD at every probe
f_c=0.10  kappa_s=+1.754e-02   OUTWARD at every probe
f_c=0.20  kappa_s=+4.125e-01   OUTWARD at every probe  (but magnitude 20× NFW)
f_c=0.30  kappa_s=-2.433e-02   INWARD at every probe   ← qualitatively wrong
f_c=0.50  kappa_s=-8.388e-03   INWARD at every probe   ← qualitatively wrong

The sign flip is uniform across the field (every probe point on every axis flips at the same f_c threshold), which rules out a local numerical artifact and points at an analytic-formula sign error in one of the branches of the cNFW deflection. The threshold is somewhere in f_c ∈ (0.20, 0.30).

Steps to reproduce 2

import numpy as np
import autogalaxy as ag

mass = 1.0e9
nfw = ag.mp.NFWTruncatedSph.from_m200_concentration(
    centre=(0.0, 0.0),
    m200_solar_mass=mass,
    concentration=12.0,
    redshift_halo=0.1,
    redshift_source=1.0,
    cosmology=ag.cosmology.Planck15(),
    truncation_factor=100.0,
)
grid = ag.Grid2D.uniform(shape_native=(101, 101), pixel_scales=0.1)

# 1. convergence_2d_from on cNFW returns identically zero
for f_c in (0.01, 0.1, 0.2, 0.5):
    cnfw = ag.profiles.mass.dark.cnfw_mcr.cNFWMCRLudlowSph(
        centre=(0.0, 0.0),
        mass_at_200=mass,
        f_c=f_c,
        redshift_object=0.1,
        redshift_source=1.0,
    )
    kappa = np.asarray(cnfw.convergence_2d_from(grid=grid).native)
    print(f"f_c={f_c:>4}: kappa_s={cnfw.kappa_s: .3e}  sum(kappa_grid)={kappa.sum(): .3e}")

# 2. Deflection magnitude scales non-monotonically with f_c
pt = ag.Grid2D.no_mask([[(1.0, 0.0)]], pixel_scales=1.0)
nfw_alpha = np.linalg.norm(np.asarray(nfw.deflections_yx_2d_from(grid=pt)))
print(f"NFW   |alpha|(r=1\")={nfw_alpha:.3e}")
for f_c in (0.01, 0.1, 0.2, 0.5):
    cnfw = ag.profiles.mass.dark.cnfw_mcr.cNFWMCRLudlowSph(
        centre=(0.0, 0.0), mass_at_200=mass, f_c=f_c,
        redshift_object=0.1, redshift_source=1.0,
    )
    a = np.linalg.norm(np.asarray(cnfw.deflections_yx_2d_from(grid=pt)))
    print(f"cNFW f_c={f_c:>4}: |alpha|(r=1\")={a:.3e}")

Expected Behaviour 2

For f_c ∈ (0, 1):

• convergence_2d_from should produce a finite, positive, smoothly varying field that integrates to the projected mass within the evaluation aperture.
• deflections_yx_2d_from magnitudes should vary continuously with f_c and approach the NFW value as f_c → 0. Beyond f_c ≈ 0.1 they currently jump by 10× before falling back to NFW-comparable values at f_c ≈ 0.5, which is not physically plausible.
• kappa_s should remain positive across the full f_c ∈ (0, 1) range.

Actual Behaviour 2

f_c=0.01: kappa_s= 4.615e-03  sum(kappa_grid)= 0.000e+00
f_c= 0.1: kappa_s= 9.652e-03  sum(kappa_grid)= 0.000e+00
f_c= 0.2: kappa_s= 9.754e-02  sum(kappa_grid)= 0.000e+00
f_c= 0.5: kappa_s=-5.118e-03  sum(kappa_grid)= 0.000e+00

NFW   |alpha|(r=1")=3.762e-03
cNFW f_c=0.01: |alpha|(r=1")=4.346e-03   (1.16 × NFW)
cNFW f_c= 0.1: |alpha|(r=1")=8.156e-03   (2.17 × NFW)
cNFW f_c= 0.2: |alpha|(r=1")=7.490e-02   (19.9 × NFW)
cNFW f_c= 0.5: |alpha|(r=1")=3.156e-03   (0.84 × NFW)

The f_c → 0 limit reduces cleanly to NFW (within ~20 % of the same deflection at the same c), so the analytic formulae are presumably correct in the small-core limit. The f_c ≈ 0.2 regime — exactly the physically interesting one for SIDM core radii — produces both a silently zero convergence map and a deflection magnitude inflated by an order of magnitude, then a sign flip in kappa_s for larger cores.

Multiple RuntimeWarning: invalid value encountered in sqrt / arctanh / arctan messages are emitted from autogalaxy.profiles.mass.dark.cnfw at lines 182, 183, 192, 193, 213, 214, 230, 231 when these routines are evaluated, which suggests the integration-region masks (radius > theta vs radius < theta) or branch selection in the cored-NFW analytic formulae may be at fault.

Environment

• OS: macOS arm64
• Python version: Python 3.12
• Package version: autolens 2026.5.21.1 (with autogalaxy dependency)

[Jammy2211]

Investigation Findings — 2026-05-27

Three distinct bugs identified, two now fixed on feature/cnfw-boundary-fix:

Bug 1: convergence_2d_from returns zeros — Already fixed

This was resolved in #449 (Phase 3 of mass profiles refactoring, merged 2026-05-26). cNFW and cNFWSph now compute convergence via MGE decomposition.

Bug 2: F_func boundary condition at θ = r_s — Fixed in this branch

cNFWSph.F_func used strict < / > masks:

mask1 = (theta > 0) & (theta < radius)   # strict <
mask2 = theta > radius                     # strict >

When theta == radius (grid point exactly at the scale radius), neither branch applies and F stays at its initial value of zero, producing a wildly wrong deflection (e.g. −0.708 instead of +0.319 at θ = r_s = 3.0″).

Fix: Changed to theta <= radius so the boundary is handled by the first branch (the arctanh expression evaluates to zero at the boundary, giving the correct limiting value). After fix: analytic matches MGE to 0.01% at θ = r_s.

Bug 3: Peñarrubia MCR integral crosses zero for f_c > ~0.18 — Mitigated with warning + clamp

The Peñarrubia (2012) mass-concentration integral factor:

F(c, f_c) = [f_c² ln(1+c/f_c) + (1−2f_c) ln(1+c)] / (1+f_c)² − c/[(1+c)(1−f_c)]

crosses zero at f_c ≈ 0.18 for typical Ludlow16 concentrations (c ≈ 7–12). Beyond this, ρ_0 = M₂₀₀ / (4π r_s³ F) diverges through infinity and becomes negative, producing:

• Deflections that blow up by ~20× at f_c ≈ 0.2
• Negative kappa_s for f_c > ~0.22
• Sign-flipped deflection field (light deflected toward lens centre)

This is a limitation of the Peñarrubia parameterisation, not a code bug per se — the cored NFW model as parameterised here cannot support cores larger than ~18% of the scale radius.

Mitigation: Added a runtime warning when the MCR factor is non-positive, and clamped to a small positive value to prevent NaN/Inf propagation. The warning explicitly states the valid f_c range.

How we found it

The Phase 1 self-consistency test suite (#124) and the new MGE potential (#449) gave us two independent reference implementations to compare against the analytic cNFWSph deflections. The boundary bug was found by sweeping θ from 0.1 to 10.0 at θ = r_s intervals and comparing cNFWSph (analytic) vs cNFW(q=1) (MGE).