Miscentering: exact single-offset kernel

This note derives the miscentered lensing observables implemented in clenspy.lensing.miscentering, for a single (delta-function) offset. The key result is a closed-form reduction of the miscentered aperture mean to a single smooth integral over the known centered mean profile — no nested cumulative integral, no integrable cusp, no endpoint singularity.

1. Definitions

Symbols used throughout:

  • \(\Sigma(R)\) — centered (azimuthally symmetric) surface density of the halo, at projected halo-centric radius \(R\).

  • \(\bar\Sigma({<}R)\) — mean of \(\Sigma\) inside the disk of radius \(R\) about the halo center,

\[ \bar\Sigma({<}R) \;=\; \frac{2}{R^{2}} \int_0^{R} \Sigma(R')\, R'\, dR' . \]

  • \(\Delta\Sigma(R) = \bar\Sigma({<}R) - \Sigma(R)\) — centered excess surface density. For NFW, \(\Sigma\) and \(\Delta\Sigma\) (hence \(\bar\Sigma = \Delta\Sigma + \Sigma\)) have the Wright & Brainerd (2000) closed forms.

  • \(R_{\rm mis}\) — distance between the true halo center and the assumed (observed) center.

  • \(R\) — projected radius measured from the assumed center.

  • \(u(t)\) — halo-centric distance of a point at radius \(R\) from the assumed center and azimuth \(t\) (law of cosines),

\[ u(t) \;=\; \sqrt{R^{2} + R_{\rm mis}^{2} - 2\,R\,R_{\rm mis}\cos t }. \]

Subscript “mis” denotes the corresponding quantity measured about the assumed center for a halo offset by \(R_{\rm mis}\).

2. Miscentered surface density

\(\Sigma_{\rm mis}\) is the azimuthal average of the centered profile over the circle of radius \(R\) about the assumed center (Yang et al. 2006; Johnston et al. 2007):

\[ \Sigma_{\rm mis}(R \mid R_{\rm mis}) \;=\; \frac{1}{\pi} \int_0^{\pi} \Sigma\!\big(u(t)\big)\, dt . \]

The integrand is smooth except when \(R = R_{\rm mis}\), where \(u(t) \to 0\) as \(t \to 0\) and NFW’s \(\Sigma \sim -\ln u\) gives an integrable logarithmic endpoint singularity (handled in §5).

3. Miscentered aperture mean: aperture-mass identity

The naive route to \(\bar\Sigma_{\rm mis}({<}R)\) is the cumulative integral of §2’s result,

\[ \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;=\; \frac{2}{R^{2}} \int_0^{R} \Sigma_{\rm mis}(R' \mid R_{\rm mis})\, R'\, dR' , \]

a nested double integral whose inner integrand has a cusp at \(R' = R_{\rm mis}\). Both problems disappear by computing the aperture mass directly in halo-centric coordinates. The mass inside the aperture (disk of radius \(R\) about the assumed center) is a weighted integral of the known \(\Sigma\) over halo-centric annuli:

\[ \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;=\; \frac{1}{\pi R^{2}} \int_0^{\infty} \Sigma(u)\, u\, \Lambda(u)\, du , \]

where \(\Lambda(u)\) is the overlap angle — the angular extent of the halo-centric circle of radius \(u\) that lies inside the aperture. With \(A(u)\) the cosine of the half-opening angle,

\[\begin{split} A(u) \;=\; \frac{u^{2} + R_{\rm mis}^{2} - R^{2}}{2\,u\,R_{\rm mis}}, \qquad \Lambda(u) \;=\; \begin{cases} 2\pi, & u \le R - R_{\rm mis} \quad (\text{only if } R_{\rm mis} < R),\\[1mm] 2\arccos A(u), & |R - R_{\rm mis}| \le u \le R + R_{\rm mis},\\[1mm] 0, & \text{otherwise.} \end{cases} \end{split}\]

This is a single integral, but \(\Sigma(u)\) still carries its central \(\ln u\) divergence when \(R_{\rm mis} < R\) (the region \(u \to 0\) contributes with weight \(2\pi\)).

4. Integration by parts onto the known mean

Since

\[ \frac{d}{du}\!\left[\frac{u^{2}}{2}\,\bar\Sigma({<}u)\right] \;=\; \Sigma(u)\, u , \]

split §3 at \(u_- = |R - R_{\rm mis}|\) and \(u_+ = R + R_{\rm mis}\) and integrate the \([u_-, u_+]\) piece by parts. The endpoint values of the overlap angle are \(\Lambda(u_+) = 0\) always, and \(\Lambda(u_-) = 2\pi\) if \(R_{\rm mis} < R\) (the circle \(u = u_-\) lies entirely inside the aperture) or \(\Lambda(u_-) = 0\) if \(R_{\rm mis} > R\). When \(R_{\rm mis} < R\), the by-parts boundary term at \(u_-\),

\[ -\,\frac{u_-^{2}}{2}\,\bar\Sigma({<}u_-)\,\Lambda(u_-) \;=\; -\,\pi\, u_-^{2}\,\bar\Sigma({<}u_-), \]

cancels the inner-disk contribution \(\int_0^{u_-} \Sigma(u)\, u \cdot 2\pi\, du = \pi\, u_-^{2}\, \bar\Sigma({<}u_-)\) exactly; when \(R_{\rm mis} > R\) both pieces are individually zero. In either case no boundary term survives:

\[ \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;=\; -\frac{1}{\pi R^{2}} \int_{u_-}^{u_+} \frac{u^{2}}{2}\, \bar\Sigma({<}u)\, \Lambda'(u)\, du , \]

with

\[ \Lambda'(u) \;=\; -\,\frac{2\,A'(u)}{\sqrt{1 - A^{2}(u)}} \;=\; -\,\frac{u^{2} + R^{2} - R_{\rm mis}^{2}} {u^{2}\, R_{\rm mis}\, \sqrt{1 - A^{2}(u)}} . \]

The central divergence of \(\Sigma\) is gone — only the smooth \(\bar\Sigma\) appears (near \(u=0\) the combination \(u^2 \bar\Sigma({<}u) \to 0\) for any profile with finite central mass) — but \(\Lambda'(u)\) still has inverse-square-root singularities at both endpoints, since

\[ 1 - A^{2}(u) \;=\; \frac{\big(u_+^{2} - u^{2}\big)\big(u^{2} - u_-^{2}\big)}{4\,u^{2}\,R_{\rm mis}^{2}} . \]

5. The substitution that removes everything

Factor the endpoint singularity using the identity above:

\[ \frac{u^{2}}{2}\,\Lambda'(u)\, du \;=\; -\,\frac{u\,\big(u^{2} + R^{2} - R_{\rm mis}^{2}\big)} {\sqrt{\big(u_+^{2} - u^{2}\big)\big(u^{2} - u_-^{2}\big)}}\; du , \]

and substitute the law-of-cosines parametrization

\[ u(t)^{2} \;=\; R^{2} + R_{\rm mis}^{2} - 2\,R\,R_{\rm mis}\cos t \;=\; \tfrac{1}{2}\big(u_+^{2} + u_-^{2}\big) - \tfrac{1}{2}\big(u_+^{2} - u_-^{2}\big)\cos t , \qquad t \in [0, \pi]. \]

Then

\[ u_+^{2} - u^{2} = \tfrac{1}{2}\big(u_+^{2} - u_-^{2}\big)(1 + \cos t), \qquad u^{2} - u_-^{2} = \tfrac{1}{2}\big(u_+^{2} - u_-^{2}\big)(1 - \cos t), \]

so the square root becomes \(\tfrac{1}{2}(u_+^{2} - u_-^{2})\sin t\), while \(2\,u\,du = \tfrac{1}{2}(u_+^{2} - u_-^{2})\sin t\, dt\) — the Jacobian cancels the singular factor exactly, leaving

\[ \boxed{\; \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;=\; \frac{1}{2\pi R^{2}} \int_0^{\pi} \Big[\, u(t)^{2} + R^{2} - R_{\rm mis}^{2} \,\Big]\, \bar\Sigma\!\big({<}u(t)\big)\, dt \;} \]

with \(u(t)\) exactly the nodes already used for \(\Sigma_{\rm mis}\) in §2. The integrand is smooth for all \(R \ne R_{\rm mis}\); at \(R = R_{\rm mis}\) the factor \(u^{2}\,\bar\Sigma({<}u) \to 0\) suppresses the (integrable) \(\ln u\) tail of \(\bar\Sigma\).

The miscentered excess surface density follows as

\[ \Delta\Sigma_{\rm mis}(R \mid R_{\rm mis}) \;=\; \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;-\; \Sigma_{\rm mis}(R \mid R_{\rm mis}) . \]

6. Analytic checks

Centered limit. \(R_{\rm mis} \to 0\): \(u(t) \to R\), the bracket is \(2R^{2}\), and the boxed formula returns \(\bar\Sigma({<}R)\).

Uniform sheet. \(\bar\Sigma \equiv \Sigma_0\): using \(\int_0^\pi \cos t\, dt = 0\),

\[ \frac{\Sigma_0}{2\pi R^{2}} \int_0^{\pi} \big(2R^{2} - 2RR_{\rm mis}\cos t\big)\, dt \;=\; \Sigma_0 \qquad \text{for every } R_{\rm mis}, \]

so a constant surface density is invariant under miscentering, as it must be — on both sides of \(R_{\rm mis} = R\), confirming the boundary cancellation of §4.

Point mass. \(\bar\Sigma({<}u) = M/(\pi u^{2})\) for a point mass \(M\) at the true center. With the standard integral

\[ \int_0^{\pi} \frac{dt}{a - b\cos t} \;=\; \frac{\pi}{\sqrt{a^{2}-b^{2}}}, \qquad a = R^{2} + R_{\rm mis}^{2},\; b = 2RR_{\rm mis},\; \sqrt{a^{2}-b^{2}} = \big|R^{2} - R_{\rm mis}^{2}\big| , \]

the boxed formula gives

\[\begin{split} \bar\Sigma_{\rm mis}({<}R \mid R_{\rm mis}) \;=\; \frac{M}{2\pi^{2} R^{2}} \left[\pi + \big(R^{2} - R_{\rm mis}^{2}\big) \frac{\pi}{\big|R^{2} - R_{\rm mis}^{2}\big|}\right] = \begin{cases} \dfrac{M}{\pi R^{2}}, & R_{\rm mis} < R,\\[2mm] 0, & R_{\rm mis} > R, \end{cases} \end{split}\]

i.e. exactly the enclosed point mass — zero when the true center lies outside the aperture.

7. The sign of \(\Delta\Sigma_{\rm mis}\) — do not clamp

\(\Delta\Sigma_{\rm mis}\) is genuinely signed: positive for \(R_{\rm mis} \ll R\) and negative for \(R_{\rm mis} \gtrsim R\) (dimensionless NFW, \(r_s = 1\), \(\Sigma = f(x)\): at \(R = 1\) it crosses from \(+0.0760\) at \(R_{\rm mis}=0.8\) to \(-0.0628\) at \(R_{\rm mis}=1.0\)). Two facts show the negative lobe is physical:

  1. The point-mass result of §6 gives \(\Delta\Sigma_{\rm mis} = 0\) exactly for \(R_{\rm mis} > R\); the negative values for an extended profile are a finite-profile gradient effect — the halo’s mass leaks into the aperture and rises toward its edge, so \(\Sigma_{\rm mis}(R) > \bar\Sigma_{\rm mis}({<}R)\).

  2. A population of halos with offsets uniform over the plane is a uniform mass sheet, which has zero lensing contrast. Hence the mean-field cancellation

\[ \int_0^{\infty} \Delta\Sigma_{\rm mis}(R \mid R_{\rm mis})\; 2\pi R_{\rm mis}\, dR_{\rm mis} \;=\; 0 , \]

and it is the negative lobe that delivers it (numerically the truncated integral falls to \(0.3\%\) of its \(L_1\) norm by \(R_{\rm mis}^{\max} = 60\,r_s\)). Clamping \(\Delta\Sigma_{\rm mis} \ge 0\) pins the integral at roughly \(+50\%\) of the \(L_1\) norm — it breaks the cancellation rather than protecting it.

8. Numerical scheme and accuracy

Both integrals of §2 and §5 are evaluated by fixed-order Gauss–Legendre quadrature after the node-clustering map

\[ t \;=\; \pi s^{2}, \qquad s \in [0, 1], \qquad dt = 2\pi s\, ds , \]

which concentrates nodes at \(t = 0\), where \(u \to |R - R_{\rm mis}|\) and the integrand of §2 develops its logarithmic behavior when \(R \approx R_{\rm mis}\). The radii are evaluated in the cancellation-free form

\[ u^{2} \;=\; \big(R - R_{\rm mis}\big)^{2} + 4\,R\,R_{\rm mis}\sin^{2}(t/2) . \]

One further ingredient is needed: the Wright & Brainerd closed forms for the NFW kernels \(f(x) = \Sigma/(2 r_s \rho_s)\) and \(g(x) = \Delta\Sigma/(r_s \rho_s)\) are \(0/0\) at \(x = 1\) and lose \(\sim 10^{-16}/|x-1|\) to cancellation nearby — and the clustered nodes land precisely there whenever \(u\) crosses the scale radius. NfwProfile therefore switches to the exact Taylor expansions for \(|x - 1| \le 10^{-2}\) (both one-sided forms continue to the same analytic series):

\[ f(1+d) = \tfrac{1}{3} - \tfrac{2}{5}d + \tfrac{13}{35}d^{2} - \tfrac{20}{63}d^{3} + \cdots, \qquad g(1+d) = \sum_{n\ge 0} \Big[\, p_n + (-1)^{n+1}\, 4(n{+}1)\ln 2 \,\Big] d^{n}, \]

with \(p_0 = \tfrac{10}{3},\; p_1 = -\tfrac{88}{15},\; p_2 = \tfrac{296}{35},\; p_3 = -\tfrac{3508}{315}, \ldots\) (nine terms kept; truncation \(\sim 10^{-19}\), direct forms outside the window accurate to \(\lesssim 2\times 10^{-14}\)).

Validation against mpmath (30-digit) adaptive quadrature of the same integrals, dimensionless NFW, \(n_{\rm nodes} = 128\)–\(256\):

quantity

typical relative error

worst case tested

\(\bar\Sigma_{\rm mis}\)

\(\lesssim 10^{-13}\)

\(10^{-10}\) at \(R = R_{\rm mis}\)

\(\Sigma_{\rm mis}\)

\(\lesssim 10^{-13}\)

\(2\times 10^{-10}\) at \(R = R_{\rm mis}\)

The scheme is a fixed matrix–vector product per radius vector (one profile evaluation on an \((n_R \times n_{\rm nodes})\) grid), with no adaptive quadrature anywhere.

9. Tabulation: making it a lookup

The DES Y3 pipeline (y3_cluster_cpp) does not evaluate §5’s quadrature per call. It reads a precomputed table, and the design is worth adopting.

9.1 The scaling that makes one table enough

\(\Delta\Sigma_{\rm mis}\) is universal in units of the scale radius. Writing \(x = R/r_s\), \(x_{\rm mis} = R_{\rm mis}/r_s\), the whole profile is

\[ \Delta\Sigma_{\rm mis}(R \mid R_{\rm mis}, M, c) = \underbrace{2\, r_s\, \delta_c\, \rho_{\rm ref}}_{\text{amplitude}} \;\times\; \underbrace{u(x, x_{\rm mis})}_{\text{dimensionless, tabulated once}} \]

with

\[ r_{200} = \left[\frac{3M}{800\pi\rho_{\rm ref}}\right]^{1/3}, \qquad r_s = \frac{r_{200}}{c}, \qquad \delta_c = \frac{200\,c^3/3}{\ln(1+c) - c/(1+c)} . \]

Mass, concentration and cosmology enter only through the analytic prefactor. One dimensionless grid therefore serves every halo, and the runtime cost collapses to an interpolation.

Note \(\rho_{\rm ref}\) sets both the boundary \(r_{200}\) and the amplitude \(\rho_s = \delta_c\rho_{\rm ref}\). Splitting those two conventions is a silent normalisation error; keep one \(\rho_{\rm ref}\) (see NfwProfile.rho_ref).

9.2 Store the value, not its log

The natural instinct is to tabulate \(\log\Delta\Sigma_{\rm mis}\) and exponentiate. This is wrong, and for the reason in §7: a log table is structurally \(\ge 0\), so it deletes the entire negative lobe at \(R_{\rm mis} > R\). The lobe is what makes the mean-field (random-point) projection term cancel — a uniformly distributed halo population is a uniform mass sheet and must give zero contrast. With the lobe zeroed it never cancels.

Store the signed linear value. The shipped table_..._deltasigma_signed_single.txt is 47.6% negative entries; that is the physics, not a bug.

9.3 Put the cusp on a grid line

\(\Delta\Sigma_{\rm mis}\) has a cusp along \(x = x_{\rm mis}\) — it crosses zero and turns over where the halo centre sits on the aperture circle. Tabulating on axes \((\ln x_{\rm mis}, \ln x)\) puts that ridge at a \(45^\circ\) angle to the grid, so bilinear interpolation cuts across it and smooths the peak off.

Measured against this module’s quadrature (which converges to 11 digits against adaptive quad), the shipped \(250 \times 1000\) table is accurate to \(\sim\!10^{-4}\) off the diagonal but degrades badly on it:

\(x = x_{\rm mis}\)

rel. error, \((\ln x_{\rm mis}, \ln x)\) axes

rel. error, ratio axes

0.01

6.5 (sign wrong)

\(8\times10^{-10}\)

0.1

0.33

\(1.3\times10^{-5}\)

0.37

0.10

\(3.9\times10^{-6}\)

1.0

0.064

\(9.0\times10^{-6}\)

10

0.12

\(6.3\times10^{-6}\)

100

0.42

\(2.9\times10^{-5}\)

The stored values are correct — at grid nodes they match this module to \(10^{-8}\)–\(10^{-5}\). The error is purely interpolation across the cusp: at \(x_{\rm mis} = 0.136\), adjacent \(x\) nodes hold \(-1.372\times10^{-2}\) and \(+1.028\times10^{-3}\), so the function crosses zero between two neighbouring samples and bilinear has nothing to work with. The two axes are also mismatched — 1000 nodes in \(x\) against 250 in \(x_{\rm mis}\) — so the ridge almost never lands on a grid line.

Fix: tabulate against the ratio. Use axes \(\big(\ln x_{\rm mis},\; \ln(x/x_{\rm mis})\big)\) with \(\ln(x/x_{\rm mis}) = 0\) an exact node. The cusp then lies along a grid line and interpolation never crosses it. The right-hand column above is a \(250 \times 241\) ratio grid — four times fewer points than the shipped table and four to nine orders of magnitude more accurate where it matters.

\(x = x_{\rm mis}\) is not a corner case: with \(R_\lambda = (\lambda/100)^{0.2}\), \(\tau_{\rm mis} = 0.17\) and \(c = 4\), a \(\lambda \sim 25\), \(M \sim 10^{14}\) cluster sits at \(x_{\rm mis} \approx 0.4\) — inside the fitted radial range, and in the projection term the “offset” sweeps all \(\theta\), so the ridge is crossed constantly.

9.4 What CLensPy does

CLensPy is table-only at runtime. It never solves these integrals during evaluation:

Module

Role

clenspy/data/nfw_miscentering.npz

the packaged grid, 1.3 MB

selection/miscentering.py

the runtime path — interpolation only

selection/miscentering_kernel.py

the §5 quadrature, offline generator

tools/make_miscentering_table.py

rebuilds the table

lensing/miscentering.py

MiscenteringProfile, reads the table

Domain \(x \in [10^{-3}, 5\times10^{3}]\), \(x_{\rm mis} \in [10^{-3}, 5\times10^{2}]\) — the DES Y3 extents, with \(x_{\rm mis}\) carried one decade lower. 300 × 841 nodes.

Which generator runs where

cluster_toolkit (the marcpaterno fork, in the y3cl_je_macos env) produces everything it can do accurately: 195 of the 300 rows, \(x_{\rm mis} \ge 0.1\). The by-parts quadrature produces the remaining 105.

The split is forced, not a preference. \(\Delta\Sigma_{\rm mis}\) is the difference \(\bar\Sigma_{\rm mis} - \Sigma_{\rm mis}\), and that difference collapses as \(x_{\rm mis}\) shrinks — at \(x_{\rm mis} = 10^{-3}\) it is \(6.6\times10^{-7}\) of \(\Sigma\) itself. ct forms it by subtracting two nearly-equal numbers, so it is cancellation-limited there:

\(x_{\rm mis}\)

ct

by-parts

ct error

1e-3

+1.01e-04

−4.33e-06

24×, wrong sign

5e-3

+2.79e-05

−7.80e-05

wrong sign

1e-2

−2.07e-04

−2.60e-04

2.0e-01

1e-1

−9.473e-03

−9.489e-03

1.6e-03

1e+0

−6.2834e-02

−6.2841e-02

1.0e-04

Refining ct’s Rsigma grid does not help — 30k → 1M points moved \(x_{\rm mis} = 10^{-2}\) only from 0.35 to 0.14, and moved \(5\times10^{-3}\) not at all. The sign error is the serious part, for §7’s reason. The seam at \(x_{\rm mis} = 0.1\) is where the two agree to \(1.65\times10^{-3}\), asserted at build time.

Reaching \(x_{\rm mis} = 10^{-3}\) with the by-parts kernel required fixing the small-\(x\) NFW kernels first (halo/nfw.py): the azimuthal ring passes through the halo centre whenever \(R = R_{\rm mis}\), and both \(f\) and \(\bar{g}\) were unusable there. See that module’s _gbarNfw.

Accuracy of the shipped table

Against the by-parts quadrature, which is self-converged to \(10^{-11}\):

region

\(\hat\Sigma\) med / max

\(\widehat{\Delta\Sigma}\) med / max

sign flips

\(x_{\rm mis} < 0.1\), global

5.7e-04 / 2.9e-03

4.3e-04 / 2.2e-03

0/120

\(x_{\rm mis} < 0.1\), on cusp

5.3e-08 / 9.2e-06

5.2e-04 / 1.1e-03

0/120

\(x_{\rm mis}\) 0.1–10, global

4.3e-04 / 2.8e-03

7.1e-04 / 1.2e-02

0/120

\(x_{\rm mis}\) 0.1–10, on cusp

5.2e-05 / 2.2e-04

3.1e-04 / 1.4e-03

0/120

\(x_{\rm mis}\) 10–500, on cusp

1.6e-04 / 2.3e-04

1.7e-03 / 1.7e-02

0/120

against 25/60 sign flips for the natural axes at equal budget. \(r_{\rm mis} = 0\) short-circuits to the analytic centred profile, so that limit is exact rather than interpolated.

One known corner: at \(\ln q \lesssim -9\) (i.e. \(x\) four or more decades below \(x_{\rm mis}\)) the values fall to \(10^{-10}\)–\(10^{-15}\) while the row scale is \(10^{-2}\), below what float32 storage resolves, so the deep wing carries sign noise. \(|\Delta\Sigma_{\rm mis}| / \Sigma < 10^{-6}\) there and it contributes nothing to any observable.

Untabulated profiles raise

require_tabulated_profile rejects anything but NfwProfile with MiscenteringTableError, at construction rather than on first evaluation. There is deliberately no quadrature fallback — falling back would reintroduce the per-call integral the table exists to remove. Einasto is the case in practice: its miscentered profile is not universal in one shape parameter the way NFW’s is, since the index \(n\) enters as a third axis, so the NFW grid cannot be reused and no Einasto grid exists.

10. References

  • Wright, C. O. & Brainerd, T. G. 2000, ApJ, 534, 34 — closed-form NFW \(\Sigma\), \(\Delta\Sigma\).

  • Yang, X. et al. 2006, MNRAS, 373, 1159 — azimuthal-average form of \(\Sigma_{\rm mis}\).

  • Johnston, D. E. et al. 2007, arXiv:0709.1159 — miscentering formalism for stacked cluster lensing.

  • Simet, M. et al. 2017, MNRAS, 466, 3103 — offset distributions for redMaPPer clusters.