diff --git a/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md b/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md
new file mode 100644
index 0000000..1360e82
--- /dev/null
+++ b/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md
@@ -0,0 +1,53 @@
+# Analytic Local Residual Moments
+
+## Objective
+
+Test the simplest analytic/asymptotic local-residual approximation for the 2D log
+Ewald split: full-space Taylor moments of the compact residual. Use it for both
+interior windows and trim-intersecting windows first, then measure where the
+boundary error appears.
+
+## Scope
+
+- Add a full-space analytic local residual estimator to the near-field split
+  sweep.
+- Compare analytic local residual values against high-order folded local
+  residual references.
+- Report interior, boundary-near, exterior-near, and vertex-near target errors.
+- Run the comparison on `ipa` and update the report.
+
+## Non-Goals
+
+- Implement boundary-aware moments.
+- Implement curved CAD fixtures.
+- Replace direct local residual references.
+
+## Acceptance Criteria
+
+- [x] Report identifies when full-space analytic moments are accurate enough and
+  when trim boundaries dominate the error.
+- [x] `make dev` passes before handoff.
+
+## Checklist
+
+- [x] Add analytic local estimator and report fields.
+- [x] Run the comparison on `ipa`.
+- [x] Update docs with conclusions.
+- [x] Validate with `make dev`.
+
+## Outcomes
+
+- The leading full-space moment `sigma^2 rho(x)/4` works for interior windows
+  when the target is well away from the trim. At `sigma=0.08`, interior analytic
+  local relative errors ranged from `2.6e-6` to `6.4e-4`.
+- The same full-space moment fails for trim-intersecting windows. Median
+  trim-near analytic local relative errors were `0.31`, `0.79`, and `1.41` for
+  `sigma=0.08`, `0.16`, and `0.32`, respectively.
+- Exterior-near targets are especially bad because the full-space moment includes
+  material outside the source region.
+- The next step is a window/boundary classifier: interior windows use analytic
+  moments; trim-intersecting windows need boundary-aware corrections.
+
+## Validation
+
+Ran `make dev` successfully after implementation and report updates.
diff --git a/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md
new file mode 100644
index 0000000..4b07f66
--- /dev/null
+++ b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md
@@ -0,0 +1,71 @@
+# Box-Code Near-Field Template Experiment
+
+## Objective
+
+Evaluate whether folded-decomposition charts can support reusable near-field
+template corrections for box-code and Volumential-style workflows.
+
+## Scope
+
+- Derive the mapped-kernel form for a 2D fan folded piece.
+- Focus the near-field form on point targets in physical space.
+- Identify singular factors that can be separated from smooth geometry payloads.
+- Frame the reuse hypothesis as a functional expansion in smooth metric fields
+  `M(z)`, not just as raw runtime geometry data.
+- Add a small analytic experiment driver with deterministic tests.
+- Document feasibility limits and the next CUTKIT/Volumential boundary.
+
+## Non-Goals
+
+- No production Volumential integration in this branch.
+- No CAD dependency for the first prototype.
+- No claim that arbitrary cut shapes have finite exact correction tables.
+
+## Acceptance Criteria
+
+- The derivation states the mapped self-interaction form and its split.
+- The derivation records the point-target form, target/source geometry taxonomy,
+  and measurement criteria from issue 61.
+- The prototype records a point-target template experiment with a direct
+  high-order reference value.
+- Tests check the experiment against an analytic or independently computed
+  invariant.
+- Documentation explains whether the idea is feasible and what data must remain
+  runtime geometry payload.
+
+## Completed Checklist
+
+- [x] Read issue 61 and repository orientation docs.
+- [x] Work through the 2D fan-map singularity derivation.
+- [x] Add an experiment module or script for the near-field template prototype.
+- [x] Add tests for the prototype invariants.
+- [x] Update docs and doc index if new documentation is added.
+- [x] Run targeted tests and `make dev`.
+- [x] Move this plan to completed with outcomes.
+
+## Outcomes
+
+- The near-field template idea is feasible as a CUTKIT experiment: singular and
+  nearly singular folded-piece interactions can be expressed on fixed template
+  domains.
+- The reusable singular part is a template-space log singular model. The metric
+  field `M(z) = DT(z)^T DT(z)` and higher-order terms are smooth geometry data.
+- The practical reuse hypothesis is that a functional expansion in `M(z)` and
+  smooth-remainder data covers the folded-decomposition cases with few modes.
+- The first prototype covers analytic straight fan maps, point-target
+  integrals, a log-kernel scale-law check, and a diagonal metric-remainder
+  sample.
+
+## Validation
+
+- `uv run --extra dev python -m pytest tests/evals/test_nearfield_templates.py tests/test_script_wrappers.py`
+- `uv run --extra dev python scripts/run_nearfield_template_experiment.py --order 6`
+- `make dev`
+
+## Remaining Risks
+
+- Template corrections may be reusable only after parameterized compression, not
+  as shape-independent lookup tables.
+- Self-interaction singular quadrature may need specialized rules beyond the
+  first direct high-order reference prototype.
+- Curved analytic pieces and adjacent-piece cases still need follow-up coverage.
diff --git a/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md b/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md
new file mode 100644
index 0000000..cd7ce7e
--- /dev/null
+++ b/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md
@@ -0,0 +1,65 @@
+# DMK Split Near-Field Experiments
+
+## Objective
+
+Run a first numerical check of the preferred near-field direction: split the 2D
+log kernel into a smooth Ewald/heat-kernel part plus a localized singular
+residual, evaluate the smooth part with folded fan quadrature, and measure when
+the residual may require boundary-aware handling.
+
+## Scope
+
+- Extend the existing straight-fan near-field prototype with a 2D log Ewald split.
+- Compare ordinary folded quadrature convergence for the full kernel versus the
+  smoothed kernel.
+- Compare current robust interior-style seeds against a simple node-barycenter
+  seed on CAD-like fan fixtures.
+- Run the numerical sweep on `ipa` and summarize results in
+  `docs/nearfield-template-experiments.md` and
+  `docs/nearfield-dmk-split-report.md`.
+
+## Non-Goals
+
+- Implement the final boundary-aware correction tables.
+- Add a production near-field API.
+- Replace existing folded decomposition seed policy.
+
+## Acceptance Criteria
+
+- [x] Experiment script emits machine-readable and Markdown summaries.
+- [x] Results include interior, boundary-near, and exterior-near target cases.
+- [x] Report states whether the smooth split improves quadrature convergence and
+  whether seed choice materially changes fan quality in the fixtures.
+- [x] `make dev` passes before handoff.
+
+## Checklist
+
+- [x] Add Ewald-split utilities and experiment dataclasses.
+- [x] Add CLI options for the DMK-style sweep and report output.
+- [x] Run the sweep on `ipa`.
+- [x] Add a concise report to the near-field notes.
+- [x] Validate with `make dev`.
+
+## Outcomes
+
+- The 2D log Ewald split made the smoothed folded-quadrature path much easier:
+  at `order=24`, smooth-part errors were at least about `70x` smaller than direct
+  full-log errors across the tested cases.
+- Median improvement was `1.8e3`, `6.6e2`, and `1.4e3` for `sigma` values
+  `0.08`, `0.16`, and `0.32`, respectively.
+- Boundary-near targets also benefited strongly, so the smooth remainder is not
+  the bottleneck in these fixtures.
+- The node-barycenter seed is viable but not uniformly better than the current
+  robust interior anchor; seed choice should remain a measured quality parameter.
+- The next unresolved issue is the localized singular residual when its window
+  intersects a rigid trim boundary.
+
+## Validation
+
+Ran `make dev` successfully after the implementation and report updates.
+
+## Risks
+
+- The straight-fan fixtures may understate CAD curvature effects.
+- Boundary-aware residual behavior requires a later local-residual experiment with
+  curved or trim-intersecting fixtures to be conclusive.
diff --git a/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md b/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md
new file mode 100644
index 0000000..3cd5c5b
--- /dev/null
+++ b/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md
@@ -0,0 +1,52 @@
+# Local Residual Near-Field Experiment
+
+## Objective
+
+Measure the compact singular residual from the DMK-style 2D log split directly,
+to decide whether ordinary folded quadrature is adequate or whether analytic or
+boundary-aware precomputed handling is needed.
+
+## Scope
+
+- Extend the near-field split sweep with local-residual folded quadrature values.
+- Report full, smooth, local, and reconstructed errors against high-order
+  references.
+- Run the residual-inclusive sweep on `ipa`.
+- Update the near-field DMK report with conclusions.
+
+## Non-Goals
+
+- Implement analytic residual moments.
+- Implement boundary-aware precomputed correction tables.
+- Add curved CAD fan fixtures.
+
+## Acceptance Criteria
+
+- [x] Residual-inclusive report identifies whether direct local folded quadrature
+  is the new bottleneck.
+- [x] Results include the same target and seed cases as the smooth split sweep.
+- [x] `make dev` passes before handoff.
+
+## Checklist
+
+- [x] Add local residual measurements to experiment dataclasses and CLI tables.
+- [x] Run the sweep on `ipa`.
+- [x] Update reports and notes.
+- [x] Validate with `make dev`.
+
+## Outcomes
+
+- Direct ordinary folded quadrature of the compact residual is not competitive as
+  the primary method.
+- At `order=24`, smooth-part median relative errors are near `1e-7`, while local
+  residual median relative errors range from `4.2e-4` to `8.9e-3` depending on
+  `sigma`.
+- Reconstructing `smooth + local` with direct local quadrature reproduces the
+  direct full-log error, confirming that the residual is the remaining hard part.
+- The next implementation choice should be a special local residual method:
+  analytic/asymptotic moments for interior windows or boundary-aware precomputed
+  folded-fan corrections for trim-intersecting windows.
+
+## Validation
+
+Ran `make dev` successfully after the implementation and report updates.
diff --git a/docs/exec-plans/completed/20260613-local-model-catalogue.md b/docs/exec-plans/completed/20260613-local-model-catalogue.md
new file mode 100644
index 0000000..c035e39
--- /dev/null
+++ b/docs/exec-plans/completed/20260613-local-model-catalogue.md
@@ -0,0 +1,37 @@
+# Local Model Catalogue
+
+## Objective
+
+Record the local analytic/asymptotic model catalogue for compact near-field
+residual corrections in 2D and 3D.
+
+## Scope
+
+- Catalog one-feature residual windows: interiors, smooth boundaries, corners,
+  faces, edges, and vertices.
+- Derive the common Taylor/modal moment expansion used by these models.
+- Record the 2D log residual radial moments used by the current experiments.
+- Link the catalogue from the near-field experiment documentation.
+
+## Acceptance Criteria
+
+- The catalogue states the routing rule that each residual support should see at
+  most one geometric singular feature.
+- 2D and 3D local models include moment formulas and leading terms where known.
+- The main near-field docs point to the catalogue as the next design reference.
+- Documentation validation passes.
+
+## Checklist
+
+- [x] Add local model catalogue.
+- [x] Add common Taylor/modal expansion.
+- [x] Add 2D log residual moment derivation.
+- [x] Link from docs index and near-field reports.
+- [x] Run validation.
+
+## Outcome
+
+Added `docs/nearfield-local-model-catalogue.md`, covering the one-feature local
+models and their expansion structure. The next implementation step is a
+window/feature classifier followed by the 2D straight-boundary half-plane moment
+experiment.
diff --git a/docs/exec-plans/completed/20260613-local-model-moment-tests.md b/docs/exec-plans/completed/20260613-local-model-moment-tests.md
new file mode 100644
index 0000000..77e32a3
--- /dev/null
+++ b/docs/exec-plans/completed/20260613-local-model-moment-tests.md
@@ -0,0 +1,64 @@
+# Local Model Moment Tests
+
+## Objective
+
+Turn the 2D straight-feature local model catalogue into executable checks for the
+log Ewald residual moments.
+
+## Scope
+
+- Add a radial moment helper for the 2D local log residual.
+- Add a sector monomial moment helper for full-plane, half-plane, and wedge
+  models at the local feature point.
+- Test the helpers against known full-plane values and independent polar
+  quadrature references.
+- Add a practical curved cut-cell check using an exact rational quadratic trim
+  and folded curve quadrature.
+- Add an exact target-centered curved-boundary model with ray clipping and finite
+  radial Ewald moments.
+- Add a high-order Taylor graph boundary model and compare it with the exact
+  circular boundary model.
+- Add a precomputation-style interpolation check for high-order smooth-boundary
+  moments.
+
+## Acceptance Criteria
+
+- Full-plane leading and second moments match closed-form values.
+- Half-plane and wedge sector moments match direct numerical quadrature within
+  the reference quadrature tolerance.
+- On a curved cut-cell sample, the half-plane model beats the full-plane model
+  for a target on the smooth trim.
+- The exact curved-boundary model agrees under refinement to better than
+  `1e-10` relative error.
+- The Taylor boundary model reaches below `1e-10` relative error at sufficiently
+  high geometry order for tested local windows.
+- The precomputed table interpolation reaches below `1e-10` relative error on
+  intermediate local-window samples.
+- The near-field targeted tests and full validation pass.
+
+## Checklist
+
+- [x] Implement radial moment helper.
+- [x] Implement sector monomial moment helper.
+- [x] Add full-plane moment tests.
+- [x] Add half-plane and wedge reference tests.
+- [x] Add curved cut-cell practical test.
+- [x] Add machine-precision curved-boundary model refinement test.
+- [x] Add high-order Taylor boundary model convergence test.
+- [x] Add precomputed boundary-table interpolation test.
+- [x] Run full validation.
+
+## Outcome
+
+Added 2D log-residual moment helpers for straight local models and tests covering
+full-plane, half-plane, and wedge sectors. The numerical reference uses a
+log-radius polar quadrature to avoid direct endpoint integration of the
+logarithmic singularity. Added a curved rational cut-cell experiment showing that
+the half-plane model is only a leading diagnostic for a target on a smooth curved
+trim; the full-plane moment incorrectly includes outside-domain mass, while
+machine precision requires exact or high-order curved boundary ray clipping with
+finite radial moments. The high-order Taylor graph model reaches the `1e-10`
+relative target on the rational quarter-circle trim at order 8 for `sigma=0.02`
+and at order 6 for `sigma=0.005, 0.01`. A first interpolation-table check over
+the normalized boundary scale also reaches the `1e-10` target, while direct
+online graph-strip quadrature is not accurate enough near the target endpoint.
diff --git a/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md b/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md
new file mode 100644
index 0000000..738e3ac
--- /dev/null
+++ b/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md
@@ -0,0 +1,49 @@
+# Complement-Subtraction Singular Quadrature
+
+## Objective
+
+Record and validate the strategy
+
+```text
+int_{Omega cap B} K(x,y) p(y) dy = int_B K(x,y) p(y) dy - int_{B \ Omega} K(x,y) p(y) dy
+```
+
+for cut-box near-field moments, using existing full-box singular tables for the
+first term and folded decomposition for the smooth complement term.
+
+## Scope
+
+- Document the mathematical route and validity conditions.
+- Identify when the complement integral is smooth and when it still needs a
+  singular/near-singular fallback.
+- Connect the route to existing full-box moment tables and folded fan pieces.
+- Keep this as a design note first; production implementation comes after the
+  classifier/table API is settled.
+
+## Acceptance Criteria
+
+- `docs/nearfield-template-experiments.md` explains complement subtraction as a
+  candidate cut-box singular quadrature route.
+- `docs/nearfield-local-model-catalogue.md` records the smoothness conditions,
+  endpoint-node caveats, and fallback cases.
+- `docs/index.md` remains current if a new standalone document is added.
+- `make dev` passes before handoff.
+
+## Checklist
+
+- [x] Add complement-subtraction derivation and algorithm sketch.
+- [x] Record smooth-complement conditions and failure modes.
+- [x] Validate documentation with `make dev`.
+- [x] Move this plan to `docs/exec-plans/completed/` with outcomes.
+
+## Outcomes
+
+- Added the full-box-minus-complement identity to
+  `docs/nearfield-template-experiments.md` as a candidate cut-box singular
+  quadrature route.
+- Recorded the key condition: the complement folded integral is smooth only when
+  the target is separated from `B \ Omega`.
+- Added routing rules and the open-quadrature caveat to
+  `docs/nearfield-local-model-catalogue.md`, including when to fall back to
+  boundary moments, wedge/cone models, refinement, or target-centered Duffy.
+- Validation passed with `make dev`.
diff --git a/docs/index.md b/docs/index.md
index ff80f0e..c3371d8 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -23,6 +23,9 @@ Repository-local documentation is the system of record for CUTKIT.
 - `docs/cad-box-batch-interface.md`: consumer-facing CAD object-or-arrays clip/integrate interface
 - `docs/volumential-handoff.md`: CUTKIT far/near primitive handoff contract for volumential workflows
 - `docs/meshmode-cut-overlay.md`: meshmode cut-overlay contract and validation semantics
+- `docs/nearfield-template-experiments.md`: folded fan near-field template derivation and prototype
+- `docs/nearfield-dmk-split-report.md`: DMK-style smooth/local split experiment results
+- `docs/nearfield-local-model-catalogue.md`: local residual model catalogue and moment expansions
 - `docs/poisson-benchmarks.md`: Poisson-oriented benchmark usage and profiles
 - `docs/poisson-galerkin-solver.md`: immersed Poisson Galerkin solve + validation workflow
 - `docs/formdsl-parity-benchmarks.md`: shared IGA + DG-SEM formdsl parity benchmark usage
diff --git a/docs/nearfield-dmk-split-report.md b/docs/nearfield-dmk-split-report.md
new file mode 100644
index 0000000..645fd42
--- /dev/null
+++ b/docs/nearfield-dmk-split-report.md
@@ -0,0 +1,193 @@
+# DMK-Style Near-Field Split Report
+
+Date: 2026-06-12
+
+Remote runner: `ipa`
+
+Command:
+
+```bash
+PYTHONPATH=src python3 scripts/run_nearfield_template_experiment.py \
+  --dmk-split-report \
+  --orders 4,6,8,10,14,18,24 \
+  --sigmas 0.08,0.16,0.32 \
+  --reference-order 128
+```
+
+Raw JSON and Markdown artifacts were generated on `ipa` in the remote run
+directory and copied to the local temporary workspace for summarization. A second
+residual-inclusive run was generated with the same parameters after adding direct
+quadrature of `G_{\mathrm{local},\sigma}`. The raw artifacts are not checked in
+because the concise tables below capture the relevant results.
+
+## Setup
+
+The experiment uses two straight-edge polygon fixtures, each integrated as a
+signed sum of folded fans from a shared seed to each boundary edge. It compares:
+
+- full 2D log-kernel folded quadrature;
+- smooth Ewald/heat-kernel complement folded quadrature;
+- the implied localized residual, measured as `full_reference - smooth_reference`.
+
+The split is
+
+$$
+G(r)=G_{\mathrm{smooth},\sigma}(r)+G_{\mathrm{local},\sigma}(r),
+\qquad
+G_{\mathrm{local},\sigma}(r)=\frac{1}{4\pi}E_1(r^2/\sigma^2),
+$$
+
+where `E1` is the exponential integral. The smoothed complement has a removable
+limit at `r=0`, so ordinary folded quadrature no longer sees the log singularity.
+
+Targets cover interior, boundary-near-inside, boundary-near-outside, and
+vertex-near-inside cases. Seed modes are the current robust interior anchor and a
+simple node barycenter.
+
+## Summary
+
+At the highest tested order, `order=24`, the smoothed folded quadrature was
+consistently more accurate than direct full-log quadrature.
+
+sigma | min improvement | median improvement | max improvement | median local/full
+--- | --- | --- | --- | ---
+0.08 | 9.252e+01 | 1.805e+03 | 3.986e+04 | 1.347e-02
+0.16 | 7.045e+01 | 6.599e+02 | 5.466e+04 | 4.753e-02
+0.32 | 1.179e+02 | 1.359e+03 | 4.176e+04 | 1.530e-01
+
+Target class | min improvement | median improvement | max improvement
+--- | --- | --- | ---
+boundary_near_inside | 9.327e+01 | 1.158e+04 | 5.466e+04
+boundary_near_outside | 7.045e+01 | 8.439e+02 | 1.479e+04
+interior | 1.479e+02 | 1.415e+03 | 4.176e+04
+vertex_near_inside | 9.366e+01 | 5.099e+02 | 1.890e+03
+
+Here `improvement = full_abs_error / smooth_abs_error` against the same
+high-order reference. The smallest improvement is still about `70x`; median
+improvements are hundreds to tens of thousands depending on target class and
+window width.
+
+Relative errors for the same `order=24` samples show the same pattern:
+
+sigma | full rel err median | smooth rel err median | local rel err median | reconstructed rel err median
+--- | ---: | ---: | ---: | ---:
+0.08 | 5.691e-05 | 9.370e-08 | 8.892e-03 | 5.691e-05
+0.16 | 5.691e-05 | 7.368e-08 | 1.629e-03 | 5.691e-05
+0.32 | 5.691e-05 | 6.844e-08 | 4.190e-04 | 5.691e-05
+
+Target class | full rel err median | smooth rel err median | local rel err median | reconstructed rel err median
+--- | ---: | ---: | ---: | ---:
+boundary_near_inside | 1.147e-03 | 7.614e-08 | 1.494e-02 | 1.147e-03
+boundary_near_outside | 4.122e-05 | 5.124e-08 | 1.629e-03 | 4.122e-05
+interior | 4.225e-05 | 4.332e-08 | 8.944e-04 | 4.225e-05
+vertex_near_inside | 1.053e-04 | 1.680e-07 | 1.235e-03 | 1.053e-04
+
+The residual-inclusive run evaluates
+
+$$
+\int G_{\mathrm{local},\sigma}(|x-y|)\rho(y)\,dy
+$$
+
+by the same ordinary folded quadrature. Its error is essentially the full-log
+error, while the reconstructed error from `smooth + local` matches the direct
+full-log error. This means the current direct local-residual quadrature has not
+solved the near-field problem; it simply moves the difficult part into an
+isolated compact term.
+
+## Seed Comparison
+
+fixture | seed mode | det ratio | edge-distance ratio
+--- | --- | --- | ---
+convex_quad | interior_anchor | 1.718e+00 | 1.510e+00
+convex_quad | node_barycenter | 2.220e+00 | 1.550e+00
+skew_quad | interior_anchor | 1.425e+00 | 1.474e+00
+skew_quad | node_barycenter | 1.697e+00 | 1.264e+00
+
+The node barycenter is viable on these convex fixtures but is not uniformly
+better. It improves the edge-distance ratio on `skew_quad` but worsens the
+boundary-determinant ratio on both fixtures. This supports treating the seed as
+a measured quality parameter rather than replacing the robust interior anchor
+yet.
+
+## Interpretation
+
+The split validates the preferred direction for the first experiment:
+
+- Ordinary folded decomposition is effective for the smoothed kernel.
+- Direct folded quadrature of the compact local residual is the new bottleneck;
+  at `order=24`, local residual median relative error ranges from `4.2e-4` to
+  `8.9e-3` depending on `sigma`, versus smooth-part median relative errors near
+  `1e-7`.
+- The singularity removal, not seed tuning, is the dominant improvement for the
+  smooth path in these fixtures.
+- Wider windows move more of the full interaction into the localized residual:
+  median `|local/full|` grows from `1.3e-2` at `sigma=0.08` to `1.5e-1` at
+  `sigma=0.32`.
+- Boundary-near targets still benefit strongly in the smooth path, so the next
+  question is not whether folded quadrature can handle the smooth remainder; it
+  can. The remaining question is how to evaluate the localized residual with an
+  analytic, asymptotic, or precomputed boundary-aware method.
+
+## Analytic Full-Space Local Moment
+
+A follow-up run tested the simplest analytic local residual estimator for
+`G_{\mathrm{local},\sigma}`: replace the clipped local domain by the full plane
+and use the leading Taylor moment
+
+$$
+u_{\mathrm{local},\sigma}(x)
+\approx \frac{\sigma^2}{4}\rho(x).
+$$
+
+For the current polynomial test density, the isotropic Laplacian correction
+vanishes, so this is the natural first full-space analytic estimate. It was used
+unchanged for interior, boundary-near, exterior-near, and vertex-near targets.
+
+At `order=24`, the analytic local relative errors were:
+
+sigma | case set | analytic local rel err min | median | max | analytic reconstructed rel err median
+--- | --- | ---: | ---: | ---: | ---:
+0.08 | all | 2.594e-06 | 2.610e-01 | 4.969e+00 | 4.023e-03
+0.08 | interior only | 2.594e-06 | 1.869e-04 | 6.371e-04 | 2.342e-06
+0.08 | trim-near only | 1.892e-01 | 3.145e-01 | 4.969e+00 | 4.891e-03
+0.16 | all | 2.446e-04 | 6.148e-01 | 2.378e+00 | 3.038e-02
+0.16 | interior only | 2.446e-04 | 8.227e-04 | 1.466e-03 | 4.383e-05
+0.16 | trim-near only | 4.213e-01 | 7.880e-01 | 2.378e+00 | 3.997e-02
+0.32 | all | 5.745e-02 | 1.060e+00 | 1.623e+00 | 1.532e-01
+0.32 | interior only | 5.745e-02 | 7.254e-02 | 8.766e-02 | 1.414e-02
+0.32 | trim-near only | 6.560e-01 | 1.409e+00 | 1.623e+00 | 1.866e-01
+
+By target class, over all tested `sigma` values and seeds:
+
+Target class | analytic local rel err min | median | max | analytic reconstructed rel err median
+--- | ---: | ---: | ---: | ---:
+boundary_near_inside | 1.892e-01 | 4.799e-01 | 7.530e-01 | 2.186e-02
+boundary_near_outside | 1.452e+00 | 2.141e+00 | 4.969e+00 | 5.235e-02
+interior | 2.594e-06 | 9.924e-04 | 8.766e-02 | 4.383e-05
+vertex_near_inside | 2.171e-01 | 7.880e-01 | 1.609e+00 | 3.997e-02
+
+This confirms the expected split:
+
+- Full-space analytic moments work for interior windows when `sigma` is small
+  compared with distance to the physical boundary.
+- The same full-space moments are not acceptable for trim-intersecting windows.
+  They include outside-domain mass and can be worse than direct local quadrature.
+- Larger `sigma` improves smooth-part quadrature but worsens the full-space local
+  approximation unless the target is far from every trim boundary.
+
+## Next Experiment
+
+Add a local-residual method that separates two cases:
+
+- interior support: use analytic Taylor/radial moments when the window is well
+  separated from the physical boundary;
+- boundary-intersecting support: build boundary-aware moments or a small
+  precomputed folded-fan correction table.
+
+The residual-inclusive and analytic-moment runs together indicate that direct
+local folded quadrature is not competitive as the primary method, while
+full-space analytic moments are only valid for interior windows. The next
+implementation should therefore add a window/boundary classifier and route
+interior windows to analytic moments, with trim-intersecting windows routed to a
+boundary-aware correction experiment. The candidate local models and moment
+expansions are catalogued in `docs/nearfield-local-model-catalogue.md`.
diff --git a/docs/nearfield-local-model-catalogue.md b/docs/nearfield-local-model-catalogue.md
new file mode 100644
index 0000000..c70504d
--- /dev/null
+++ b/docs/nearfield-local-model-catalogue.md
@@ -0,0 +1,775 @@
+# Near-Field Local Model Catalogue
+
+This note catalogs the local analytic/asymptotic models for compact residual
+kernels after a DMK/Ewald-style split. The guiding rule is:
+
+> choose the residual support so each local correction sees at most one geometric
+> singular feature.
+
+If a residual window sees multiple unrelated trim arcs, corners, faces, or
+vertices, it should be shrunk, the source geometry should be subdivided, or the
+case should fall back to a precomputed/direct local correction.
+
+## Common Expansion Form
+
+Let the compact residual kernel be radial:
+
+$$
+K_\sigma(z) = K_{\mathrm{local},\sigma}(|z|),
+\qquad z = y-x.
+$$
+
+The local residual contribution is
+
+$$
+u_{\mathrm{local},\sigma}(x)
+= \int_{\Omega_x} K_\sigma(y-x)a(y)\,dy,
+$$
+
+where `a` denotes density times any smooth geometric/Jacobian factor and
+`Omega_x` is the source region inside the effective residual support around the
+target.
+
+Choose a local model feature `F` and local coordinates `z`. Expand the smooth
+amplitude at the target or closest feature point:
+
+$$
+a(x+z)
+= \sum_{|\alpha|\le K}
+\frac{1}{\alpha!}\partial^\alpha a(x)z^\alpha
++ R_{K+1}(z).
+$$
+
+Then
+
+$$
+u_{\mathrm{local},\sigma}(x)
+= \sum_{|\alpha|\le K}
+\frac{1}{\alpha!}\partial^\alpha a(x)
+M_\alpha(F,\sigma)
++ \mathcal R_{K+1},
+$$
+
+with feature moments
+
+$$
+M_\alpha(F,\sigma)
+= \int_{\Omega_F} K_\sigma(z)z^\alpha\,dz.
+$$
+
+For curved features, `Omega_F` is itself expanded around a leading flat, wedge,
+or cone model. The resulting moments are products or perturbations of radial
+moments and angular moments.
+
+Equivalently, introduce the scaled coordinate
+
+$$
+z = \sigma q,
+\qquad K_\sigma(z)=\sigma^{-s}k(|q|),
+$$
+
+where `s` is the kernel scaling power. Then
+
+$$
+M_\alpha(F,\sigma)
+= \sigma^{d-s+|\alpha|}
+\int_{\widehat\Omega_F} k(|q|)q^\alpha\,dq,
+$$
+
+up to corrections from finite support and curved geometry. For the 2D log
+residual used here, `s = 0`, so the leading mass is `O(sigma^2)`. For 3D kernels
+whose residuals scale like `sigma^{-1} k(r/sigma)`, the leading mass is also
+`O(sigma^2)`; the exact power should be attached to the selected kernel split.
+
+In a local modal implementation, write the smooth amplitude as
+
+$$
+a(x+z) \approx \sum_{\beta\in\mathcal B_K} c_\beta \phi_\beta(z),
+$$
+
+where the modal basis may be monomials, orthogonal polynomials on the local
+feature model, or chart-jet modes. The correction is then
+
+$$
+u_{\mathrm{local},\sigma}(x)
+\approx \sum_{\beta\in\mathcal B_K} c_\beta
+\mu_\beta(F,\sigma),
+\qquad
+\mu_\beta(F,\sigma)=\int_{\Omega_F}K_\sigma(z)\phi_\beta(z)\,dz.
+$$
+
+The Taylor form below is the monomial version of this modal expansion.
+
+## Feature Catalogue
+
+| Dimension | Feature count in support | Leading model | Parameters |
+| --- | --- | --- | --- |
+| 2D | none | plane | target offset, amplitude jet |
+| 2D | one smooth trim | half-plane | signed distance, tangent, curvature jet |
+| 2D | one corner/junction | wedge | offset, opening angle, edge curvature jets |
+| 3D | none | space | target offset, amplitude jet |
+| 3D | one smooth face | half-space | signed distance, frame, principal curvature jet |
+| 3D | one edge | dihedral wedge | offset, dihedral angle, face/edge jets |
+| 3D | one vertex | polyhedral cone | offset, incident-face cone, edge/face jets |
+
+The catalogue is complete for one geometric singularity because local CAD/cut
+features stratify by codimension: interiors, smooth boundary strata, edges or
+corners, and vertices. Multi-feature windows are deliberately excluded from one
+model and must be split or handled by fallback.
+
+## 2D Log Residual Moments
+
+For the 2D log split used in the first experiments,
+
+$$
+K_\sigma(r)
+= \frac{1}{4\pi}E_1(r^2/\sigma^2).
+$$
+
+Define radial moments
+
+$$
+R_m(\sigma)
+= \int_0^\infty K_\sigma(r)r^{m+1}\,dr.
+$$
+
+Using
+
+$$
+\int_0^\infty u^q E_1(u)\,du = \frac{\Gamma(q+1)}{q+1},
+$$
+
+and `u = r^2/sigma^2`,
+
+$$
+R_m(\sigma)
+= \frac{\sigma^{m+2}}{8\pi}
+\frac{\Gamma\left(\frac{m}{2}+1\right)}{\frac{m}{2}+1}.
+$$
+
+In particular,
+
+$$
+2\pi R_0(\sigma) = \frac{\sigma^2}{4},
+$$
+
+which gives the leading full-plane estimate used in the experiments:
+
+$$
+u_{\mathrm{local},\sigma}(x)
+\approx \frac{\sigma^2}{4}a(x).
+$$
+
+For finite support truncation at `r <= c sigma`, replace `R_m` by
+
+$$
+R_m(c,\sigma)
+= \int_0^{c\sigma} K_\sigma(r)r^{m+1}\,dr.
+$$
+
+These truncated moments can be evaluated with recurrence/incomplete-gamma forms
+or precomputed as one-dimensional functions of `c`.
+
+## 2D Local Models
+
+### Interior Model
+
+Condition: the residual support lies inside the source region and away from all
+trim features.
+
+Model domain:
+
+$$
+\Omega_F = \mathbb R^2.
+$$
+
+The moments are full angular moments:
+
+$$
+M_{ij}^{\mathrm{int}}
+= \int_0^{2\pi}\int_0^\infty
+K_\sigma(r)(r\cos\theta)^i(r\sin\theta)^j r\,dr\,d\theta.
+$$
+
+Odd moments vanish by symmetry. The first terms are
+
+$$
+\begin{aligned}
+u_{\mathrm{int}}(x)
+&= 2\pi R_0 a(x)
++ \frac{1}{2}\left(\partial_{xx}a+\partial_{yy}a\right)
+\pi R_2 + O(\sigma^6) \\
+&= \frac{\sigma^2}{4}a(x)
++ \frac{\sigma^4}{32}\Delta a(x)
++ O(\sigma^6)
+\end{aligned}
+$$
+
+for the 2D log residual.
+
+### Smooth Boundary Model
+
+Condition: the residual support intersects one smooth trim segment and no corner.
+
+Use local tangent-normal coordinates at the closest trim point `p`:
+
+$$
+z = \tau t + \eta n,
+\qquad x-p = d_t\tau + d_n n.
+$$
+
+For a flat leading model with the source on `eta >= -d_n`, use target-centered
+coordinates
+
+$$
+z = s(\cos\theta,\sin\theta),
+$$
+
+with angular/radial domain
+
+$$
+s\sin\theta \ge -d_n.
+$$
+
+Equivalently, along a ray `theta`, the lower radial limit is
+
+$$
+s_0(\theta;d_n)=
+\begin{cases}
+0, & \sin\theta \ge 0, \\
+\frac{-d_n}{\sin\theta}, & \sin\theta < 0,
+\end{cases}
+$$
+
+for an interior target with `d_n >= 0`. Exterior targets reverse which angular
+sector contributes. The moments are
+
+$$
+M_{ij}^{\mathrm{half}}(d_n,\sigma)
+= \int_{\Theta(d_n)}
+\int_{s_0(\theta)}^\infty
+K_\sigma(s)(s\cos\theta)^i(s\sin\theta)^j s\,ds\,d\theta.
+$$
+
+For `d_n = 0`, this reduces to half-plane moments. The leading term is one half
+of the full-plane mass:
+
+$$
+u_{\mathrm{half}}(x)
+\approx \frac{\sigma^2}{8}a(x)
+$$
+
+when the target is exactly on a straight boundary.
+
+For curved trims, express the boundary as a graph in tangent-normal coordinates:
+
+$$
+\eta = h(\tau)
+= \frac{\kappa}{2}\tau^2
++ \frac{h_3}{6}\tau^3
++ \cdots.
+$$
+
+The source condition becomes
+
+$$
+d_n + s\sin\theta
+\ge h(d_t+s\cos\theta).
+$$
+
+Expand the boundary perturbation around the flat model. This produces curvature
+corrections of the schematic form
+
+$$
+M_\alpha^{\mathrm{bdry}}
+= M_\alpha^{\mathrm{half}}
++ \kappa C_{\alpha,1}(d/\sigma)\sigma^{|\alpha|+3}
++ h_3 C_{\alpha,2}(d/\sigma)\sigma^{|\alpha|+4}
++ \cdots.
+$$
+
+The functions `C` are universal coefficient maps in normalized target offset and
+can be tabulated.
+
+### Corner Model
+
+Condition: the residual support intersects one polygon corner or one trim
+junction, and no other geometric singular feature.
+
+The leading model is a wedge with opening angle `omega`:
+
+$$
+\Omega_F = \{(r,\theta): r\ge 0,
+\theta_0\le \theta\le \theta_1\},
+\qquad \omega = \theta_1-\theta_0.
+$$
+
+For a target at the corner, moments are separable:
+
+$$
+M_{ij}^{\mathrm{wedge}}
+= R_{i+j}(\sigma)
+\int_{\theta_0}^{\theta_1}
+(\cos\theta)^i(\sin\theta)^j\,d\theta.
+$$
+
+The leading mass term is
+
+$$
+u_{\mathrm{wedge}}(x)
+\approx \frac{\omega}{2\pi}\frac{\sigma^2}{4}a(x).
+$$
+
+For a target offset from the corner, rays have angle-dependent entry and exit
+limits from the two boundary rays. The same formula applies with radial limits
+`s_0(theta), s_1(theta)`. Curved edges meeting at the corner introduce two graph
+perturbations and a corner-angle perturbation series.
+
+## 3D Local Models
+
+The same taxonomy applies in 3D. Let
+
+$$
+z = r\omega,
+\qquad \omega\in S^2.
+$$
+
+Define radial moments
+
+$$
+R_m^{3D}(\sigma)
+= \int_0^\infty K_\sigma(r)r^{m+2}\,dr.
+$$
+
+Then feature moments are radial moments times angular moments over subsets of
+the sphere, plus curvature corrections.
+
+### Interior Model
+
+Condition: support lies inside the source volume and away from faces, edges, and
+vertices.
+
+Model domain:
+
+$$
+\Omega_F = \mathbb R^3.
+$$
+
+Moments are
+
+$$
+M_\alpha^{\mathrm{int},3D}
+= R_{|\alpha|}^{3D}(\sigma)
+\int_{S^2}\omega^\alpha\,d\omega.
+$$
+
+Odd moments vanish. The first correction is proportional to the Laplacian of the
+amplitude.
+
+### Smooth Face Model
+
+Condition: support intersects one smooth boundary face and no edge.
+
+Use local coordinates `(u,v,n)` at the closest face point. The flat leading model
+is a half-space:
+
+$$
+n \ge -d_n.
+$$
+
+Moments are spherical cap moments:
+
+$$
+M_\alpha^{\mathrm{face}}
+= \int_{S^2}\int_{s_0(\omega)}^\infty
+K_\sigma(s)(s\omega)^\alpha s^2\,ds\,d\omega.
+$$
+
+For a target exactly on a flat face, the leading mass term is half of the
+full-space mass. Curvature corrections use the face graph
+
+$$
+n = h(u,v)
+= \frac{1}{2}(\kappa_1 u^2 + \kappa_2 v^2)
++ \cdots.
+$$
+
+The coefficient maps depend on normalized offset `d/sigma` and the principal
+curvatures scaled by `sigma`.
+
+### Edge Model
+
+Condition: support intersects one edge where two faces meet and no vertex.
+
+The leading model is a dihedral wedge extruded along the edge tangent. In local
+cylindrical coordinates around the edge,
+
+$$
+z = z_e e + r(\cos\theta n_1 + \sin\theta n_2),
+\qquad \theta_0\le\theta\le\theta_1.
+$$
+
+Equivalently, the angular domain on `S^2` is a lune determined by the two face
+planes. Moments are
+
+$$
+M_\alpha^{\mathrm{edge}}
+= \int_{\Omega_{\mathrm{dihedral}}}
+K_\sigma(|z|)z^\alpha\,dz.
+$$
+
+For a target on a straight edge, the leading mass is the solid-angle fraction of
+the full-space mass:
+
+$$
+u_{\mathrm{edge}}(x)
+\approx \frac{\Omega_{\mathrm{edge}}}{4\pi}M_0^{\mathrm{full}}a(x),
+$$
+
+where `Omega_edge` is the solid angle of the dihedral wedge. Higher terms depend
+on angular monomial moments of the spherical lune. Curved faces and curved edge
+tangents add perturbation terms from the two surface graphs and the edge curve.
+
+### Vertex Model
+
+Condition: support intersects one vertex and no other independent feature.
+
+The leading model is a polyhedral cone:
+
+$$
+\Omega_F = \{r\omega: r\ge 0,
+\omega\in A\subset S^2\}.
+$$
+
+Moments are
+
+$$
+M_\alpha^{\mathrm{vertex}}
+= R_{|\alpha|}^{3D}(\sigma)
+\int_A \omega^\alpha\,d\omega.
+$$
+
+The leading mass term is again the solid-angle fraction of the full-space mass:
+
+$$
+u_{\mathrm{vertex}}(x)
+\approx \frac{|A|}{4\pi}M_0^{\mathrm{full}}a(x).
+$$
+
+For curved CAD vertices, the cone receives perturbations from all incident face
+charts and edge curves. These are likely table-driven rather than hand-derived in
+early experiments.
+
+## Routing Rules
+
+Each residual correction should be routed by the local feature set inside the
+effective support radius `R_sigma = c sigma`:
+
+1. No physical feature in support: use interior moments.
+2. Exactly one smooth boundary face or trim segment: use smooth-boundary moments.
+3. Exactly one corner in 2D or edge in 3D: use wedge or dihedral moments.
+4. Exactly one vertex in 3D: use cone moments.
+5. More than one feature or ambiguous classification: shrink `sigma`, subdivide
+   the source region, or use a direct/precomputed local fallback.
+
+The normalized parameters for reusable tables are:
+
+- normalized target offset `d/sigma`;
+- opening angle or solid-angle data for wedges/cones;
+- curvature coefficients scaled by powers of `sigma`;
+- amplitude and Jacobian jets;
+- chart orientation and sign.
+
+## Complement-Subtraction Route
+
+For cut boxes with an ambient tensor-product density representation, a useful
+alternative to direct boundary-aware singular moments is to subtract the outside
+of-domain complement from an existing full-box singular table. For a box `B`,
+physical region `Omega_B = Omega cap B`, source mode `p_alpha`, and target `x`,
+
+$$
+\int_{\Omega_B} K(x,y)p_\alpha(y)\,dy
+= \int_B K(x,y)p_\alpha(y)\,dy
+- \int_{B\setminus\Omega}K(x,y)p_\alpha(y)\,dy.
+$$
+
+The first term is the standard interior-box moment. The complement term can be
+evaluated by folded decomposition as an ordinary smooth integral if
+
+$$
+\operatorname{dist}(x, B\setminus\Omega) > 0.
+$$
+
+This condition holds for targets strictly inside the cut cell with a positive
+distance to the cut boundary. It fails, or becomes numerically fragile, when the
+target lies on the physical boundary, when the complement closure contains the
+target, or when the target-complement distance is small relative to the requested
+accuracy and quadrature order.
+
+For the small-distance case, the preferred fix is a QBX-style expansion of the
+complement potential rather than direct high-resolution quadrature at the target.
+Although the integrand `K(x,y)` is nearly singular when `x` approaches the cut
+boundary, the complement potential is smooth on the physical side away from the
+complement sources. Choose an expansion center `c` in `Omega_B` and a radius `r`
+such that
+
+$$
+|x-c| < r < \operatorname{dist}(c, B\setminus\Omega).
+$$
+
+Then compute local expansion coefficients for
+
+$$
+u_{\mathrm{comp}}(x)=\int_{B\setminus\Omega}K(x,y)p_\alpha(y)\,dy
+$$
+
+by folded quadrature over the complement with target/center `c`, where the
+coefficient integrands are smooth. The expansion is then evaluated at the actual
+near-boundary target `x` and subtracted from the full-box singular moment.
+
+This is a volume-potential analogue of QBX: the expansion is not used to represent
+the singular self term on `Omega_B`; that part is already in the full-box table.
+It is used only to avoid nearly singular evaluation of the smooth complement
+field near the boundary.
+
+Routing for this route is:
+
+1. If the target is target-separated from the complement, use full-box table minus
+   folded complement.
+2. If the target is close to the complement boundary but covered by a safe
+   expansion ball, use full-box table minus complement-QBX evaluation.
+3. If no safe complement expansion center exists and the target lies on one
+   smooth boundary feature, use smooth-boundary moments or a boundary table.
+4. If the target lies on an edge, corner, or vertex, route to the corresponding
+   wedge/cone model or refine.
+5. If multiple complement pieces approach the target, shrink the box/window,
+   subdivide the geometry, or use a target-centered fallback.
+
+The complement folded quadrature should use open quadrature rules on fan or cone
+coordinates. Lobatto endpoint nodes are useful for interpolation and element
+interfaces, but they are undesirable for this smooth-complement integral because a
+fan endpoint, clipped-boundary endpoint, or coning apex may lie on the target or
+on a nearly singular geometric feature. Open rules keep all quadrature nodes in
+the smooth interior of each complement piece.
+
+This route is attractive for List 1 matrix generation because the singular part
+is geometry independent and can reuse the same full-box tables as uncut boxes.
+The geometry-dependent object is the complement moment matrix, which is smooth
+under the separation condition and can be built by folded quadrature, by
+QBX-style local expansion coefficients, or compressed over cut-geometry
+parameters.
+
+## Precomputed Smooth-Boundary Moment Scheme
+
+For a target on or near one smooth boundary feature, introduce local
+tangent-normal coordinates at the closest boundary point:
+
+$$
+y = x + s\tau + n\nu,
+$$
+
+where `tau` is tangent and `nu` is the inward normal. The flat model is the
+half-plane `n >= 0`. The curved boundary is represented as a local graph
+
+$$
+n = h(s)
+= c_2 s^2 + c_3 s^3 + c_4 s^4 + \cdots.
+$$
+
+The local source domain is
+
+$$
+n \ge h(s).
+$$
+
+For a compact residual kernel `K_sigma`, the local correction is
+
+$$
+I_\sigma
+= \int_{n\ge h(s)} K_\sigma(s,n)a(s,n)\,ds\,dn,
+$$
+
+where `a` is density times Jacobian and any smooth chart factor. Rewrite this as
+a flat half-plane moment minus the curved strip excluded by the boundary:
+
+$$
+I_\sigma
+= \int_{n\ge 0} K_\sigma(s,n)a(s,n)\,ds\,dn
+- \int_{0\le n\le h(s)}K_\sigma(s,n)a(s,n)\,ds\,dn.
+$$
+
+This identity is the conceptual split. The final implementation should not do a
+naive online tensor quadrature over the strip; the logarithmic endpoint at the
+target makes that too slow for machine precision. Instead, normalize by the
+residual scale:
+
+$$
+s = \sigma q,
+\qquad n = \sigma p.
+$$
+
+Then the scaled boundary is
+
+$$
+p = H(q;\lambda)
+= \lambda_2 q^2 + \lambda_3 q^3 + \lambda_4 q^4 + \cdots,
+$$
+
+with normalized geometry coefficients
+
+$$
+\lambda_m = c_m\sigma^{m-1}.
+$$
+
+For the 2D log residual, `K_sigma(s,n) = k(q,p)` after scaling, and
+
+$$
+I_\sigma
+= \sigma^2
+\int_{p\ge H(q;\lambda)} k(q,p)
+a(x+\sigma(q\tau+p\nu))\,dq\,dp.
+$$
+
+Expand the smooth amplitude in the same local coordinates:
+
+$$
+a(x+\sigma(q\tau+p\nu))
+= \sum_{i+j\le A} a_{ij}\sigma^{i+j}q^i p^j
++ O(\sigma^{A+1}).
+$$
+
+The correction becomes a finite contraction:
+
+$$
+I_\sigma
+\approx
+\sum_{i+j\le A}a_{ij}\sigma^{2+i+j}
+M_{ij}(\lambda),
+$$
+
+where the normalized smooth-boundary moments are
+
+$$
+M_{ij}(\lambda)
+= \int_{p\ge H(q;\lambda)} k(q,p)q^i p^j\,dq\,dp.
+$$
+
+These `M_ij(lambda)` are the reusable objects. They can be represented in either
+of two equivalent ways:
+
+1. tabulate `M_ij(lambda)` over normalized boundary-jet parameters and
+   interpolate;
+2. expand the moment function itself in geometry parameters:
+
+$$
+M_{ij}(\lambda)
+= M_{ij}^{\mathrm{half}}
++ \lambda_2 M_{ij}^{(2)}
++ \lambda_3 M_{ij}^{(3)}
++ \lambda_2^2 M_{ij}^{(22)}
++ \cdots.
+$$
+
+The current prototype tests the first representation on a rational quarter-circle
+trim. Seven precomputed samples in the normalized scale interpolate intermediate
+smooth-boundary moment values below `1e-10` relative error.
+
+### Offline Table Generation
+
+For each kernel split, expansion order, and model family:
+
+1. choose a normalized parameter box for `lambda`;
+2. choose a monomial or modal basis `q^i p^j` up to amplitude order `A`;
+3. for each table node in `lambda`, evaluate `M_ij(lambda)` using a robust
+   high-order method such as exact ray clipping and finite radial moments;
+4. store interpolation coefficients or a compressed modal representation;
+5. validate the table against direct high-order local references on curved CAD
+   fixtures.
+
+The expensive ray-exit solves belong in this offline generation step, not in the
+main runtime path.
+
+### Runtime Evaluation
+
+For each target/source local residual interaction:
+
+1. classify the support as interior, one smooth boundary, corner/wedge, or
+   fallback;
+2. if smooth-boundary, find the closest local boundary point and frame
+   `(tau, nu)`;
+3. compute the scaled boundary coefficients `lambda_m = c_m sigma^(m-1)`;
+4. compute the amplitude jet coefficients `a_ij` for density, Jacobian, and chart
+   factors;
+5. interpolate or evaluate the precomputed moments `M_ij(lambda)`;
+6. return the contraction
+
+$$
+\sum_{i+j\le A}a_{ij}\sigma^{2+i+j}M_{ij}(\lambda).
+$$
+
+For higher dimensions, the same idea applies with local coordinates split into
+tangential variables and one normal variable. Smooth faces use graph moments;
+edges use wedge/dihedral graph moments; vertices use cone-sector moments.
+
+## Experiment Order
+
+The recommended implementation sequence is:
+
+1. 2D interior full-space moments. Done for the leading term.
+2. 2D straight-boundary half-plane moments.
+3. 2D polygon-corner wedge moments.
+4. 2D curved-boundary perturbations.
+5. 3D face half-space moments.
+6. 3D edge dihedral moments.
+7. 3D vertex cone moments.
+
+At each stage, compare the analytic model against the high-order local residual
+reference and only add precomputed tables where closed-form or low-dimensional
+moment evaluation is not sufficient.
+
+## Implemented Checks
+
+The first executable checks cover straight-feature 2D moments for the log Ewald
+residual:
+
+- full-plane mass and second moments;
+- boundary-point half-plane sector moments;
+- corner-point wedge sector moments.
+
+The implementation exposes the radial moment and angular-sector monomial moment
+helpers in `cutkit.evals.nearfield_templates`. Tests compare closed-form sector
+moments against direct polar quadrature of the local residual kernel.
+
+A second practical check uses an exact rational quadratic quarter-circle cut cell
+with one curved trim and two straight cell faces. For targets on the smooth curved
+trim and `sigma = 0.02, 0.04`, the leading flat half-plane local model is only a
+zeroth-order diagnostic: it is percent-level accurate, while the full-plane model
+is worse than `100%` relative error because it includes outside-domain mass.
+
+The production-quality boundary-aware model must include the local curved
+geometry. For the same rational curved cut cell, the implemented target-centered
+model clips each ray against the exact curved cell and evaluates the radial Ewald
+moments analytically along each ray. The remaining angular integral is adaptively
+resolved, and the test requires refinement agreement below `1e-10` relative
+error. This is the practical conclusion: flat half-plane moments are useful for
+routing and leading-order estimates, but machine precision requires exact or
+high-order curved boundary geometry plus finite radial moments.
+
+The first high-order boundary-model experiment confirms this route on the same
+smooth rational trim. The local circular boundary is replaced by an order-`K`
+Taylor graph in tangent-normal coordinates. Along each target-centered ray, the
+model solves for the positive graph exit distance and uses finite radial Ewald
+moments. For `sigma = 0.02`, the relative error against the exact circular
+boundary model decreases from `O(1e-6)` at quadratic order to `O(1e-11)` at order
+6 and below `1e-10` at order 8. For smaller windows `sigma = 0.005, 0.01`, order
+6 is already below `1e-10` relative error.
+
+The first precomputation-style check tabulates the high-order smooth-boundary
+moment value as a function of the normalized local boundary scale. Seven table
+samples over `sigma = 0.004..0.028` interpolate intermediate samples at
+`sigma = 0.010, 0.018, 0.026` with maximum relative error below `1e-10`. This is
+not yet the final modal table format, but it shows that the online ray solve can
+move into an offline/table-generation phase. A direct tensor graph-strip
+quadrature was also tested and is not sufficient for machine precision near the
+target logarithmic endpoint; the precomputed object should be a normalized moment
+or coefficient table, not a naive online strip quadrature.
diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md
new file mode 100644
index 0000000..bdd22b3
--- /dev/null
+++ b/docs/nearfield-template-experiments.md
@@ -0,0 +1,698 @@
+# Near-Field Template Experiments
+
+This note records the first mathematical check for issue 61: whether folded
+decomposition charts can support a reusable near-field template story for later
+box-code or Volumential work.
+
+## Bezier Fan Chart
+
+For the mathematical model, start with the chart type expected from CAD-backed
+folded decomposition: a fan from a seed point `V` to a smooth Bezier trim curve
+`C(t)`. A degree-`p` Bezier edge is
+
+$$
+\begin{aligned}
+C(t) &= \sum_{a=0}^p B_a^p(t)P_a, \\
+B_a^p(t) &= \binom{p}{a}(1-t)^{p-a}t^a,
+\qquad 0\le t\le 1.
+\end{aligned}
+$$
+
+The folded chart is
+
+$$
+\begin{aligned}
+T(r,t) &= (1-r)V + rC(t),
+\qquad 0\le r,t\le 1.
+\end{aligned}
+$$
+
+Its Jacobian is
+
+$$
+J(r,t) = r\det(C(t)-V,C'(t)).
+$$
+
+The factor `r` is the Duffy-style apex degeneracy already used by
+`cutkit.quadrature.folded2d`. It is not a new singularity for physical
+integration; it cancels area at the apex. The straight-edge case is only the
+degree-1 specialization and should be treated as a smoke-test fixture, not as
+the primary derivation.
+
+## Mapped Kernel
+
+For the 2D Laplace kernel
+
+$$
+G(x,y) = -\frac{1}{2\pi}\log |x-y|,
+$$
+
+the point-target near-field integral is
+
+$$
+u(x) = \int_{[0,1]^2} G(x,T_e(r,t))\rho(T_e(r,t))J_e(r,t)\,dr\,dt.
+$$
+
+This is the primary target model for the box-code experiment. The source is a
+folded piece, but targets are physical-space discretization nodes in the source
+box and neighboring near-field boxes. There is no target folded piece in the
+first near-field experiment.
+
+For reusable tables, express the target relative to the same physical box or
+chart payload. If `x0` is a chart/box reference point and `H` is a box-scale map,
+write
+
+$$
+x = x_0 + H\tau,
+$$
+
+where `tau` is the template-space location of a target node in the containing or
+neighboring box. The source-to-point integral is then a fixed-domain integral in
+the source chart coordinates and the target offset parameter `tau`:
+
+$$
+u(\tau) = \int_{[0,1]^2}G(x_0+H\tau,T_e(\xi))\rho(T_e(\xi))J_e(\xi)\,d\xi.
+$$
+
+The experiment should classify target/source geometry as:
+
+- target node on or very near the source folded piece;
+- target node in the source piece's containing box;
+- target node in an edge-neighboring box;
+- target node in a vertex-neighboring box;
+- well-separated target node, where ordinary tensor-product quadrature should already work.
+
+## Point-Target Singular Split
+
+For a target point near the source chart, choose a template coordinate `z_x`
+whose mapped source point `T(z_x)` is closest to `x`, or otherwise represents the
+local source point responsible for the near-singular behavior. Write nearby
+source coordinates as
+
+$$
+\begin{aligned}
+\xi &= z_x + \delta, \\
+d &= x - T(z_x).
+\end{aligned}
+$$
+
+Here `delta` is the local source displacement away from the nearest source point,
+and `d` is the physical target offset from that point. The kernel depends on the
+physical displacement from the source point to the target:
+
+$$
+x - T(\xi) = x - T(z_x+\delta).
+$$
+
+The purpose of the singular split is to approximate this displacement accurately
+near `delta=0`, while keeping the approximation expressed on the fixed source
+template. Instead of keeping only the first metric term, replace the chart by its
+order-`K` Taylor jet around `z_x`:
+
+$$
+\begin{aligned}
+T(z_x+\delta)
+&= T(z_x) + \sum_{1\le |\alpha|\le K} \\
+&\quad \frac{1}{\alpha!}\partial^\alpha T(z_x)\delta^\alpha \\
+&\quad + R_{K+1}(\delta).
+\end{aligned}
+$$
+
+Substituting this jet into `x - T(z_x+delta)` gives a polynomial approximation to
+the target-to-source displacement. This is the critical object to tabulate
+against:
+
+$$
+\begin{aligned}
+P_K(\delta;z_x,d)
+&= d - \sum_{1\le |\alpha|\le K} \\
+&\quad \frac{1}{\alpha!}\partial^\alpha T(z_x)\delta^\alpha.
+\end{aligned}
+$$
+
+Thus `P_K` is not an extra geometric map; it is the high-order local model of
+the physical vector `x - T(xi)`. Its coefficients are exactly the chart jet at
+`z_x` plus the target offset `d`.
+
+The singular model for the point-target kernel is then
+
+$$
+\begin{aligned}
+G_K^{\mathrm{sing}}(x,z_x,\delta)
+&= -\frac{1}{2\pi}\log |P_K(\delta;z_x,d)|.
+\end{aligned}
+$$
+
+For an on-surface target, `d=0`. For an off-surface target, `d` carries the
+normal and tangential target offset. Increasing `K` moves curvature and mixed
+terms from the remainder into the singular model. If `C(t)` is a degree-`p`
+Bezier curve, the fan map `T(r,t)=(1-r)V+rC(t)` is a polynomial of total degree
+`p+1` in `(r,t)`, so the Taylor jet terminates once `K >= p+1`. For
+rational/NURBS-derived charts, `K` is a truncation order chosen by the requested
+accuracy.
+
+The reusable part is a family of template-space singular model integrals or
+correction operators parameterized by the finite jet data
+$\{\partial^\alpha T(z_x)\}$ and the target offset `d`. The runtime payload remains
+map coefficients, Jacobian data, source coefficients, target-node offsets, and
+high-order correction or interpolation data.
+
+## Metric-Field Expansion Tables
+
+The table-building objective is to separate fixed singular template integrals
+from per-chart smooth geometry and per-target offset coefficients. The table is
+not built from the metric term alone. It is built from the high-order
+target-to-source displacement polynomial `P_K`.
+
+For a target node represented by offset coordinates `tau`, write source points
+near the closest chart point as `xi = z_x + delta`. The singular distance model is
+
+$$
+|x-T(z_x+\delta)|^2 \approx |P_K(\delta;z_x,d)|^2.
+$$
+
+The coefficients of `P_K` are the runtime geometry payload:
+
+$$
+\mathcal J_K(z_x,d)
+= \left(d,\{\partial^\alpha T(z_x):1\le |\alpha|\le K\}\right).
+$$
+
+This includes both on-surface and off-surface targets. On-surface targets have
+`d=0`; off-surface targets have nonzero `d`, which may include both normal and
+tangential offset components. The metric field is only the first quadratic part
+of this larger jet. In the special case `K=1` and `d=0`, the model reduces to
+
+$$
+\begin{aligned}
+|P_1(\delta;z_x,0)|^2 &= \delta^T M(z_x)\delta, \\
+M(z_x) &= DT(z_x)^TDT(z_x).
+\end{aligned}
+$$
+
+That special case is useful for sanity checks, but it is not the intended table
+construction. The practical table construction chooses a reference jet
+$\mathcal J_0$ for a geometry/target-offset bin and expands the actual jet around
+it:
+
+$$
+\mathcal J_K(z_x,d) = \mathcal J_0 + \Delta\mathcal J(z_x,d).
+$$
+
+Equivalently, write the displacement polynomial as a reference model plus a
+smooth perturbation:
+
+$$
+\begin{aligned}
+P_K(\delta;z_x,d)
+&= P_K^0(\delta) + \Delta P_K(\delta;z_x,d).
+\end{aligned}
+$$
+
+Then the log kernel can be expanded around the reference displacement:
+
+$$
+\begin{aligned}
+\log |P_K(\delta;z_x,d)|
+&= \log |P_K^0(\delta)| \\
+&\quad + \frac{1}{2}\log\left(1+R_K(\delta;z_x,d)\right),
+\end{aligned}
+$$
+
+where
+
+$$
+\begin{aligned}
+R_K(\delta;z_x,d)
+&= \frac{N_K(\delta;z_x,d)}{|P_K^0(\delta)|^2}, \\
+N_K(\delta;z_x,d)
+&= 2P_K^0(\delta)\cdot\Delta P_K(\delta;z_x,d) \\
+&\quad + |\Delta P_K(\delta;z_x,d)|^2.
+\end{aligned}
+$$
+
+When each bin is chosen so the perturbation ratio is controlled, this expression
+can be expanded in powers, interpolation modes, or a low-rank basis in the jet
+perturbation coefficients. The resulting singular model has the schematic form
+
+$$
+\begin{aligned}
+G_K^{\mathrm{sing}}(x,z_x,\delta)
+&\approx \sum_\alpha c_\alpha(z_x,\tau)S_\alpha(\delta;\mathcal J_0).
+\end{aligned}
+$$
+
+Here `S_alpha` are fixed template functions for the selected reference jet, and
+`c_alpha(z_x,tau)` are smooth functions of the target offset `d`, the metric
+entries, and all higher-order chart derivatives included in `P_K`. For Bezier
+fan charts, these jet coefficients are smooth functions of the seed point and
+Bezier control points.
+
+For one target node and one source density expansion mode `rho_j`, the singular
+contribution becomes
+
+$$
+\begin{aligned}
+u_j^{\mathrm{sing}}(\tau)
+&\approx \sum_\alpha \int_{[0,1]^2} \\
+&\quad \rho_j(\xi)c_\alpha(z_x,\tau) \\
+&\quad \times S_\alpha(\xi-z_x;\mathcal J_0)J(\xi)\,d\xi.
+\end{aligned}
+$$
+
+If the smooth per-chart/per-target factor is also expanded in a template basis
+`p_beta`,
+
+$$
+\begin{aligned}
+c_\alpha(z_x,\tau)J(\xi)
+&\approx \sum_\beta q_{\alpha\beta}(\tau)p_\beta(\xi),
+\end{aligned}
+$$
+
+then the runtime evaluation uses precomputed source-template tables:
+
+$$
+\begin{aligned}
+u_j^{\mathrm{sing}}(\tau)
+&\approx \sum_{\alpha,\beta}q_{\alpha\beta}(\tau)T_{j\alpha\beta}.
+\end{aligned}
+$$
+
+$$
+\begin{aligned}
+T_{j\alpha\beta}
+&= \int_{[0,1]^2} \\
+&\quad \rho_j(\xi)p_\beta(\xi)S_\alpha(\xi-z_x;\mathcal J_0)\,d\xi.
+\end{aligned}
+$$
+
+The tensors `T_{j alpha beta}` are precomputed on fixed source template domains,
+or tabulated over a small target-offset grid. At runtime, each folded chart and
+target node only supply the coefficients `q_{alpha beta}(tau)` from smooth
+jet/Jacobian fields, target offset data, and any higher-order correction or
+remainder representation.
+
+## Gaussian Window And Asymptotic Local Part
+
+A more controlled near-field experiment is to split the kernel before deciding
+what must be tabled. Choose a length scale `sigma`, usually tied to the local box
+size or fan diameter, and a smooth radial window `w_sigma(r)` that is near one at
+`r=0` and rapidly decays to numerical zero for `r` larger than a few `sigma`.
+For the 2D Laplace kernel, write
+
+$$
+\begin{aligned}
+G(x,y) &= G_{\mathrm{sing},\sigma}(x,y) + G_{\mathrm{smooth},\sigma}(x,y), \\
+G_{\mathrm{sing},\sigma}(x,y) &= w_\sigma(|x-y|)G(x,y), \\
+G_{\mathrm{smooth},\sigma}(x,y) &= \left(1-w_\sigma(|x-y|)\right)G(x,y).
+\end{aligned}
+$$
+
+The smooth part can be handled by ordinary folded-decomposition quadrature on
+the source fan:
+
+$$
+\begin{aligned}
+u_{\mathrm{smooth},\sigma}(x)
+&= \int_{[0,1]^2}G_{\mathrm{smooth},\sigma}(x,T(\xi)) \\
+&\quad \times \rho(T(\xi))J(\xi)\,d\xi.
+\end{aligned}
+$$
+
+The singular-windowed part is numerically local:
+
+$$
+\begin{aligned}
+u_{\mathrm{sing},\sigma}(x)
+&= \int_{[0,1]^2}G_{\mathrm{sing},\sigma}(x,T(\xi)) \\
+&\quad \times \rho(T(\xi))J(\xi)\,d\xi.
+\end{aligned}
+$$
+
+Because `G_{sing,sigma}` is compactly supported to numerical tolerance, this
+term matters only when the physical target point lies inside the fan piece or
+within the chosen window radius of it. For target nodes in neighboring boxes that
+are outside this support, the singular local contribution is skipped and the
+ordinary folded quadrature of the smooth part is sufficient.
+
+This changes the hard problem substantially. The windowed singular part only
+needs a special local treatment on a restricted target set:
+
+- target points on or inside the source fan;
+- target points within a few `sigma` of the fan boundary;
+- target offsets whose support intersects the fan under the high-order local
+  displacement model `P_K`.
+
+Following the DMK-style strategy, this local treatment should not be ordinary
+folded quadrature of the singular kernel. Jiang and Greengard's continuous DMK
+construction splits the finest-level local interaction into Gaussian-sum
+convolutions, accelerated by separation of variables, plus a highly localized
+correction handled by Taylor/asymptotic expansion of the source density and
+low-dimensional radial or spherical moments. For CUTKIT, the analogous local
+coordinates should build those moments from `P_K`, the Gaussian-window scale, and
+smooth expansions of the density and Jacobian. Precomputation, if used, should
+target reusable asymptotic moments/coefficient maps for this local expansion,
+not a generic table for the whole near-field integral.
+
+One important difference from the continuous DMK setting is the role of physical
+boundaries. DMK's continuous-source derivation treats smooth densities on a
+rectangular box, so artificial leaf-box boundaries do not introduce local jumps
+and physical boundary effects are compatible with smooth tapering or extension of
+the density. A CUTKIT cut boundary is rigid: the effective density is multiplied
+by a domain indicator and can jump across the trim. The first experiment should
+therefore try the DMK-inspired split directly in folded-decomposition
+coordinates: evaluate the smoothed kernel with ordinary folded quadrature and
+evaluate the localized correction only for targets whose window support reaches a
+source fan. Measurements should then decide whether a special boundary-aware
+local expansion is actually required for windows intersecting the cut boundary,
+or whether direct local folded quadrature is adequate at the requested
+accuracies.
+
+### Full-Box Minus Complement Route
+
+A complementary route is available for cut boxes when the source density is
+represented on the ambient box basis. Let `B` be a source leaf box and let
+`Omega_B = Omega cap B`. For a target `x in Omega_B` and a local source mode
+`p_j(y)`, the desired cut-box moment is
+
+$$
+I_j^{\Omega_B}(x)
+= \int_{\Omega_B}K(x,y)p_j(y)\,dy.
+$$
+
+If a Volumential-style or boxcode-style singular table already provides the
+full-box moment
+
+$$
+I_j^B(x)=\int_B K(x,y)p_j(y)\,dy,
+$$
+
+then
+
+$$
+I_j^{\Omega_B}(x)
+= I_j^B(x)-I_j^{B\setminus\Omega}(x),
+\qquad
+I_j^{B\setminus\Omega}(x)
+=\int_{B\setminus\Omega}K(x,y)p_j(y)\,dy.
+$$
+
+The first term contains the point singularity and is handled by the existing
+interior-box table. The second term is smooth when the target is separated from
+the complement `B \setminus Omega`. In that case, folded decomposition can
+integrate the complement with ordinary high-order quadrature, using open rules on
+fan coordinates so no quadrature node lands on a boundary endpoint or fan apex
+that could coincide with the target.
+
+This identity turns a singular cut-cell integral into:
+
+1. one reusable singular full-box table lookup;
+2. one smooth folded-complement integral over the outside-of-domain portion of the
+   same box;
+3. a subtraction in the same source basis and target normalization.
+
+It is especially attractive for local List 1 matrices on cut boxes because the
+full-box table is geometry independent and the geometry dependence moves to a
+smooth complement payload. It is not a replacement for boundary-singular local
+models when the target lies on the cut boundary, the complement closure contains
+the target, or the target-complement distance is too small for ordinary folded
+quadrature to be efficient. Those cases still need boundary-aware moment tables,
+target-centered Duffy fallback, smaller boxes, or a residual window that sees only
+one feature.
+
+There is a useful intermediate case before falling back to singular boundary
+moments. The complement contribution
+
+$$
+u_{\mathrm{comp}}(x)=\int_{B\setminus\Omega}K(x,y)p_j(y)\,dy
+$$
+
+is smooth as a function of `x` throughout the physical cut cell, even when `x` is
+very close to the cut boundary and direct quadrature of the nearly singular
+integrand is expensive. This is the same opportunity used by QBX: choose an
+expansion center `c` inside `Omega_B`, separated from the complement by a safe
+radius, compute a local expansion of `u_comp` about `c`, and evaluate that
+expansion at near-boundary target nodes covered by the expansion ball.
+
+In this version, the complement term becomes:
+
+1. folded quadrature of expansion coefficients at a center with safe separation
+   from `B \setminus Omega`;
+2. local Taylor, harmonic, or kernel-specific QBX evaluation at the physical
+   target nodes;
+3. subtraction from the full-box singular table.
+
+This preserves the main advantage of the complement route: the singularity is
+still handled by the full-box table, while the geometry-dependent complement is
+handled through smooth coefficient integrals. The expansion center and order must
+be chosen so the expansion ball covers the target but does not cross the
+complement source region.
+
+Target-centered Duffy refolding is useful as a diagnostic or single-target
+fallback, but it should not be treated as the main reusable strategy. Folded
+decomposition does not require a fixed seed; for one point target `x` inside the
+cut region, one can choose the decomposition anchor to be `x` itself and refold
+the local source panel around that target. Each resulting fan has the form
+
+$$
+\begin{aligned}
+T_x(r,t) &= (1-r)x + rC(t).
+\end{aligned}
+$$
+
+The point singularity is then exactly at the Duffy apex `r=0`, and the Jacobian
+contributes the usual factor `r`. For the logarithmic kernel, the local integrand
+has the form
+
+$$
+\begin{aligned}
+G(x,T_x(r,t))J_x(r,t)
+&\sim -\frac{1}{2\pi}\log(r|C(t)-x|)\,r\,\det(C(t)-x,C'(t)),
+\end{aligned}
+$$
+
+which is integrable in the Duffy coordinate, but not smooth. High-order accuracy
+can still require many radial nodes, and the refolding must be rebuilt for each
+target point. This makes target-centered Duffy unsuitable as the primary
+near-field template mechanism. The main experiment should instead test whether a
+shared DMK-inspired local correction, parameterized by target offsets and chart
+jets, can amortize over many physical target nodes. Target-centered Duffy remains
+a useful reference calculation for validating the local correction on individual
+interior targets.
+
+The smooth remainder no longer needs singular quadrature or local asymptotics; it
+is evaluated directly with the same signed folded quadrature machinery used for
+far-field source clouds. The open design choices are the window family, the scale
+`sigma`, the asymptotic expansion order, and the criterion used to skip the local
+singular treatment for targets outside the window support. This DMK-style split
+is the preferred experiment direction: ordinary folded decomposition handles the
+smoothed kernel, while only the highly localized residual needs special local
+treatment.
+
+For the Bezier fan map
+
+$$
+T(r,t) = (1-r)V + rC(t),
+$$
+
+the metric entries are explicit:
+
+$$
+\begin{aligned}
+T_r &= C(t)-V, \\
+T_t &= rC'(t).
+\end{aligned}
+$$
+
+$$
+M(r,t) =
+\left[
+\begin{matrix}
+|C(t)-V|^2 & r(C(t)-V)\cdot C'(t) \\
+r(C(t)-V)\cdot C'(t) & r^2|C'(t)|^2
+\end{matrix}
+\right].
+$$
+
+For Bezier or piecewise-Bezier CAD trims, these entries are smooth
+low-parameter functions of `(r,t)`, the seed `V`, and the Bezier control points
+`P_a`. They are the first terms of the larger displacement-jet payload. The
+mathematical reason precomputed template tables may apply is not that `M` alone
+is universal, but that the full finite jet $\mathcal J_K$, the Jacobian, and the
+target offset vary smoothly across the folded-decomposition chart family.
+
+The seed point should also be treated as an experimental degree of freedom.
+CUTKIT's current 2D folded quadrature chooses a deterministic strictly interior
+anchor for robustness, but near-field quality may improve with other stable
+choices. The first alternative should be simple: use the barycenter of the trim
+nodes or sampled curve nodes as `V`, then compare it with the current interior
+anchor. Later seed-quality criteria can minimize fan aspect ratio, Jacobian
+variation, signed-weight cancellation, or the spread of the local jet payloads.
+
+## Feasibility Result
+
+The answer is more favorable with a Gaussian-window split. Singular or nearly
+singular source-folded-piece to physical-point interactions can be moved to fixed
+source template domains, and the local singular treatment can be restricted to
+targets inside or very close to the fan piece. The preferred decomposition is
+therefore not a broad near-field table. It is a DMK-inspired split: direct folded
+quadrature for the smoothed kernel and a compact local correction for the
+singular residual. When the local window lies in a smooth interior neighborhood,
+the correction should be analytically derived from Taylor/asymptotic moments. If
+the local window intersects a rigid cut boundary, the correction may need
+boundary-aware precomputation over folded fan templates.
+
+The practical hypothesis is now about the compactness of the full point-target
+payload, not only the metric field. For each near target, the relevant runtime
+data is
+
+$$
+\mathcal J_K(z_x,d),\quad J(\xi),\quad \rho(\xi),\quad \tau,
+$$
+
+where $\mathcal J_K$ contains the target offset and all chart derivatives used by
+`P_K`. Because Bezier fan maps have smooth, low-parameter structure away from
+collapsed seed faces, these jet coefficients should vary smoothly across the
+folded-decomposition cases produced by CAD-like trims. A functional expansion in
+the jet data, target-offset data, and Jacobian/density factors may therefore
+cover enough practical cases to make precomputed near-field tables worthwhile.
+The Gaussian-window split improves the odds because the local asymptotic
+treatment does not need to represent weakly near or well-separated target
+interactions; those move to the smooth folded-quadrature path. Target-centered
+Duffy refolding should be used only as a validation baseline or fallback because
+it is target-specific and leaves a nonsmooth integrable kernel.
+
+The open numerical question is whether the observed set of jets and target
+offsets is compact or low-rank enough after binning by reference jet
+$\mathcal J_0$.
+If many asymptotic terms, bins, boundary-aware tables, or local coefficient modes
+are required even after windowing, the technique may not be worthwhile even
+though the template formulation is mathematically valid.
+
+## Measurements
+
+The first DMK-style smooth/local split sweep is summarized in
+`docs/nearfield-dmk-split-report.md`. It confirms that ordinary folded
+decomposition is effective for the smoothed kernel: at `order=24`, the smoothed
+folded-quadrature error was at least about `70x` smaller than direct full-log
+quadrature error across the tested interior, boundary-near, and vertex-near
+targets. The residual-inclusive follow-up shows that direct ordinary folded
+quadrature of the compact local residual is the remaining bottleneck: its error
+tracks the direct full-log error while the smooth-part error is near `1e-7`
+relative. The next unresolved question is therefore which special residual path
+to implement first: analytic/asymptotic interior moments or boundary-aware
+precomputed folded-fan corrections for trim-intersecting windows.
+
+The first analytic-moment check used the full-space leading residual moment
+`sigma^2 rho(x)/4`. It works for interior windows when the target is well away
+from the trim: at `order=24`, interior analytic local relative errors were
+`2.6e-6` to `6.4e-4` for `sigma=0.08`. The same formula fails for
+trim-intersecting windows, where median analytic local relative errors were
+`0.31` for `sigma=0.08`, `0.79` for `sigma=0.16`, and `1.41` for `sigma=0.32`.
+The next implementation should therefore classify whether the residual window is
+interior to the source region. Interior windows can use analytic moments;
+trim-intersecting windows need boundary-aware corrections. The detailed local
+model catalogue and moment expansions are recorded in
+`docs/nearfield-local-model-catalogue.md`.
+
+The full-box-minus-complement route should be tested alongside these residual
+models. Its expected sweet spot is a cut-box target that is interior to
+`Omega cap B` but close enough to the cut boundary that direct folded quadrature
+of the singular kernel is poor. The full-box singular table supplies the singular
+moment exactly for the ambient basis, and the complement integral should converge
+as a smooth folded integral as long as the complement is target-separated. For
+near-boundary targets where direct complement quadrature is nearly singular, test
+the complement-QBX variant: compute expansion coefficients from folded complement
+quadrature at a separated center and evaluate the expansion at the target.
+
+For each geometry and interaction case, record:
+
+- reference value;
+- ordinary pulled-forward tensor-product quadrature error;
+- template singular-correction error;
+- Gaussian-window smooth-part quadrature error;
+- Gaussian-window singular asymptotic expansion error;
+- higher-order correction/remainder approximation error after the asymptotic
+  local part is removed;
+- sensitivity to the window scale `sigma` and support cutoff;
+- dependence on quadrature order;
+- dependence on target offset, including on-surface, near-surface, containing-box,
+  and neighbor-box target nodes;
+- dependence on seed location and Bezier curve coefficients, starting with the
+  current robust interior anchor versus the barycenter of trim or sampled curve
+  nodes;
+- size, rank, and smoothness of the local asymptotic coefficient data;
+- how many window/asymptotic orders, bins, and jet modes are needed for the
+  observed folded-decomposition cases;
+- whether the localized singular residual can be handled by analytic interior
+  moments, or whether boundary-aware precomputed fan-template corrections are
+  needed when the window intersects a trim boundary;
+- whether a full-box singular moment minus a smooth folded complement integral is
+  more accurate or cheaper than direct boundary-aware residual tables for
+  target-separated cut-box interactions;
+- whether complement-QBX expansions remove the close-to-boundary resolution
+  burden in the complement integral, and how expansion center placement/order
+  should be chosen relative to the cut boundary;
+- whether metric-only organization is sufficient for any subfamily, or whether
+  higher-order jet coefficients dominate the correction size.
+- how often source-box and neighbor-box target nodes actually require the
+  singular local treatment after the window-support test.
+- how much accuracy target-centered Duffy refolding provides as a single-target
+  baseline for the same node count, and when it is too expensive to use as a
+  fallback.
+
+The first geometry family should be CAD-independent but CAD-realistic: quadratic
+and cubic Bezier fan charts with seed locations chosen to produce positive,
+negative, and folded orientations. Straight fan pieces are useful only as
+low-level smoke tests. OpenCascade should be used later as a source of
+stress-test Bezier/NURBS-derived fixtures, not as a dependency of the numerical
+question.
+
+## Prototype
+
+`cutkit.evals.nearfield_templates` implements:
+
+- `FanTemplateMap2D` for the current straight-edge smoke-test fan map.
+- `point_target_laplace_potential(...)` for point-target source integrals with
+  polynomial template densities.
+- `diagonal_remainder_sample(...)` for the metric singular split.
+- `run_nearfield_template_experiment(...)` for the baseline feasibility check.
+
+The script wrapper is:
+
+```bash
+uv run python scripts/run_nearfield_template_experiment.py --order 12
+```
+
+The current smoke-test experiment checks the exact scale law for the 2D log
+kernel. If both the source fan and physical target point are scaled by `lambda`,
+then
+
+$$
+\begin{aligned}
+u_\lambda(\lambda x)
+&= \lambda^2\left(u(x) - \frac{\log(\lambda)m_\rho}{2\pi}\right).
+\end{aligned}
+$$
+
+where `m_rho` is the density-weighted signed source mass on the unscaled fan
+piece. It also records diagonal remainder samples that shrink as the source and
+target template coordinates coalesce and a point-target low-order-vs-reference
+error.
+
+## Next Decision
+
+The idea is feasible as a CUTKIT experiment, with these limits:
+
+- Far-field source clouds can remain ordinary signed quadrature sources.
+- Near-field correction reuse should target a Gaussian-windowed point-target
+  local asymptotic expansion for targets inside or very close to each fan piece.
+  The smooth remainder should use ordinary folded-decomposition quadrature.
+  The singular residual should use analytic moments where the local window is
+  interior to the source region; if it intersects a trim boundary, test whether
+  boundary-aware precomputed folded-fan corrections are necessary.
+  Target-centered Duffy refolding should be kept as a single-target validation
+  baseline or fallback, not as the main reusable correction strategy.
+- Volumential should still own tree/list composition; CUTKIT should export the
+  local geometry/operator payloads needed by those lists.
diff --git a/scripts/run_nearfield_template_experiment.py b/scripts/run_nearfield_template_experiment.py
new file mode 100755
index 0000000..6fbf473
--- /dev/null
+++ b/scripts/run_nearfield_template_experiment.py
@@ -0,0 +1,211 @@
+#!/usr/bin/env python3
+"""Run the folded fan near-field template experiment."""
+
+from __future__ import annotations
+
+import argparse
+import json
+from dataclasses import asdict
+from pathlib import Path
+
+from cutkit.evals import (
+    DmkSplitExperimentReport,
+    run_dmk_split_nearfield_experiment,
+    run_nearfield_template_experiment,
+)
+
+
+def _parse_int_tuple(value: str) -> tuple[int, ...]:
+    items = tuple(int(item.strip()) for item in value.split(",") if item.strip())
+    if not items:
+        raise argparse.ArgumentTypeError("expected at least one integer")
+    return items
+
+
+def _parse_float_tuple(value: str) -> tuple[float, ...]:
+    items = tuple(float(item.strip()) for item in value.split(",") if item.strip())
+    if not items:
+        raise argparse.ArgumentTypeError("expected at least one float")
+    return items
+
+
+def _parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "--order",
+        type=int,
+        default=12,
+        help="Gauss-Legendre order for the target template coordinates.",
+    )
+    parser.add_argument(
+        "--scale-factor",
+        type=float,
+        default=1.75,
+        help="Positive geometry scale factor for the log-kernel scale-law check.",
+    )
+    parser.add_argument(
+        "--dmk-split-report",
+        action="store_true",
+        help="Run the DMK-style smooth/local split experiment instead of the baseline.",
+    )
+    parser.add_argument(
+        "--orders",
+        type=_parse_int_tuple,
+        default=(4, 6, 8, 10, 14, 18),
+        help="Comma-separated quadrature orders for --dmk-split-report.",
+    )
+    parser.add_argument(
+        "--sigmas",
+        type=_parse_float_tuple,
+        default=(0.08, 0.16, 0.32),
+        help="Comma-separated Ewald split widths for --dmk-split-report.",
+    )
+    parser.add_argument(
+        "--reference-order",
+        type=int,
+        default=96,
+        help="Reference quadrature order for --dmk-split-report.",
+    )
+    parser.add_argument(
+        "--json-output",
+        type=Path,
+        help="Optional JSON output path for --dmk-split-report.",
+    )
+    parser.add_argument(
+        "--markdown-output",
+        type=Path,
+        help="Optional Markdown output path for --dmk-split-report.",
+    )
+    return parser.parse_args()
+
+
+def _format_dmk_split_report(report: DmkSplitExperimentReport) -> str:
+    lines: list[str] = []
+    lines.append("# DMK-Style Near-Field Split Report")
+    lines.append("")
+    lines.append(f"reference_order: `{report.reference_order}`")
+    lines.append(f"orders: `{', '.join(str(order) for order in report.orders)}`")
+    lines.append(f"sigmas: `{', '.join(f'{sigma:g}' for sigma in report.sigmas)}`")
+    lines.append("")
+    lines.append("## Seed Quality")
+    lines.append("")
+    lines.append(
+        "fixture | seed_mode | seed | min_det | det_ratio | min_edge_distance | edge_distance_ratio"
+    )
+    lines.append("--- | --- | --- | --- | --- | --- | ---")
+    for sample in report.seed_quality:
+        lines.append(
+            f"{sample.fixture} | {sample.seed_mode} | {sample.seed} | "
+            f"{sample.min_abs_boundary_det:.6e} | "
+            f"{sample.max_to_min_abs_boundary_det:.6e} | "
+            f"{sample.min_edge_distance:.6e} | "
+            f"{sample.max_to_min_edge_distance:.6e}"
+        )
+
+    lines.append("")
+    lines.append("## Highest-Order Error Summary")
+    lines.append("")
+    lines.append(
+        "fixture | seed_mode | target | sigma | full_error | smooth_error | local_error | analytic_local_error | reconstructed_error | analytic_reconstructed_error | smooth_improvement | local_to_full"
+    )
+    lines.append(
+        "--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---"
+    )
+    max_order = max(report.orders)
+    for sample in report.samples:
+        if sample.order != max_order:
+            continue
+        improvement = sample.full_abs_error / max(sample.smooth_abs_error, 1.0e-300)
+        lines.append(
+            f"{sample.fixture} | {sample.seed_mode} | {sample.target_label} | "
+            f"{sample.sigma:.3g} | {sample.full_abs_error:.6e} | "
+            f"{sample.smooth_abs_error:.6e} | {sample.local_abs_error:.6e} | "
+            f"{sample.analytic_local_abs_error:.6e} | "
+            f"{sample.reconstructed_abs_error:.6e} | "
+            f"{sample.analytic_reconstructed_abs_error:.6e} | "
+            f"{improvement:.6e} | "
+            f"{sample.local_to_full_ratio:.6e}"
+        )
+
+    lines.append("")
+    lines.append("## Error By Order")
+    lines.append("")
+    lines.append(
+        "fixture | seed_mode | target | sigma | order | full_error | smooth_error | local_error | analytic_local_error | reconstructed_error | analytic_reconstructed_error | smooth_improvement"
+    )
+    lines.append(
+        "--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---"
+    )
+    for sample in report.samples:
+        improvement = sample.full_abs_error / max(sample.smooth_abs_error, 1.0e-300)
+        lines.append(
+            f"{sample.fixture} | {sample.seed_mode} | {sample.target_label} | "
+            f"{sample.sigma:.3g} | {sample.order} | {sample.full_abs_error:.6e} | "
+            f"{sample.smooth_abs_error:.6e} | {sample.local_abs_error:.6e} | "
+            f"{sample.analytic_local_abs_error:.6e} | "
+            f"{sample.reconstructed_abs_error:.6e} | "
+            f"{sample.analytic_reconstructed_abs_error:.6e} | {improvement:.6e}"
+        )
+    lines.append("")
+    return "\n".join(lines)
+
+
+def _run_dmk_split_report(args: argparse.Namespace) -> int:
+    report = run_dmk_split_nearfield_experiment(
+        orders=args.orders,
+        sigmas=args.sigmas,
+        reference_order=args.reference_order,
+    )
+    markdown = _format_dmk_split_report(report)
+    if args.json_output is not None:
+        args.json_output.parent.mkdir(parents=True, exist_ok=True)
+        args.json_output.write_text(json.dumps(asdict(report), indent=2) + "\n")
+    if args.markdown_output is not None:
+        args.markdown_output.parent.mkdir(parents=True, exist_ok=True)
+        args.markdown_output.write_text(markdown)
+    print(markdown)
+    return 0
+
+
+def main() -> int:
+    args = _parse_args()
+    if args.dmk_split_report:
+        return _run_dmk_split_report(args)
+
+    result = run_nearfield_template_experiment(
+        order=args.order,
+        scale_factor=args.scale_factor,
+    )
+
+    print("field | value")
+    print("--- | ---")
+    print(f"order | {result.order}")
+    print(f"near_point_target | {result.near_point_target}")
+    print(f"point_target_reference | {result.point_target_reference:.16e}")
+    print(f"point_target_low_order | {result.point_target_low_order:.16e}")
+    print(f"point_target_abs_error | {result.point_target_abs_error:.16e}")
+    print(f"signed_area | {result.signed_area:.16e}")
+    print(f"density_mass | {result.density_mass:.16e}")
+    print(f"scale_factor | {result.scale_factor:.16e}")
+    print(f"scaled_near_point_target | {result.scaled_near_point_target}")
+    print(
+        f"scaled_point_target_potential | {result.scaled_point_target_potential:.16e}"
+    )
+    print(
+        "expected_scaled_point_target_potential | "
+        f"{result.expected_scaled_point_target_potential:.16e}"
+    )
+    print(f"scaled_abs_error | {result.scaled_abs_error:.16e}")
+    print("")
+    print("delta | physical_distance | model_distance | remainder")
+    print("--- | --- | --- | ---")
+    for sample in result.diagonal_remainders:
+        print(
+            f"{sample.delta:.1e} | {sample.physical_distance:.16e} | "
+            f"{sample.model_distance:.16e} | {sample.remainder:.16e}"
+        )
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/src/cutkit/evals/__init__.py b/src/cutkit/evals/__init__.py
index 110c262..1d3bd92 100644
--- a/src/cutkit/evals/__init__.py
+++ b/src/cutkit/evals/__init__.py
@@ -74,6 +74,40 @@
     run_poisson_galerkin_benchmark,
     solve_trimmed_poisson_galerkin,
 )
+from cutkit.evals.nearfield_templates import (
+    DiagonalRemainderSample,
+    CurvedBoundaryLocalModelReport,
+    CurvedBoundaryLocalModelSample,
+    DmkSplitExperimentReport,
+    DmkSplitSample,
+    FanTemplateMap2D,
+    GraphBoundaryModelReport,
+    GraphBoundaryModelSample,
+    NearfieldTemplateExperiment,
+    PrecomputedBoundaryTableReport,
+    PrecomputedBoundaryTableSample,
+    SeedQualitySample,
+    TaylorBoundaryModelReport,
+    TaylorBoundaryModelSample,
+    diagonal_remainder_sample,
+    ewald_log_local_kernel,
+    ewald_log_local_radial_moment_2d,
+    ewald_log_local_sector_monomial_moment_2d,
+    ewald_log_smooth_kernel,
+    expected_scaled_laplace_point_potential,
+    exponential_integral_e1,
+    laplace_log_kernel,
+    metric_model_distance,
+    point_target_laplace_potential,
+    point_target_smooth_ewald_potential,
+    run_curved_boundary_local_model_experiment,
+    run_dmk_split_nearfield_experiment,
+    run_graph_boundary_model_experiment,
+    run_nearfield_template_experiment,
+    run_precomputed_boundary_table_experiment,
+    run_taylor_boundary_model_experiment,
+    template_density_mass,
+)
 
 __all__ = [
     "case_to_fixture_payload",
@@ -138,6 +172,38 @@
     "build_poisson_galerkin_geometry_snapshot",
     "run_poisson_galerkin_benchmark",
     "solve_trimmed_poisson_galerkin",
+    "DiagonalRemainderSample",
+    "CurvedBoundaryLocalModelReport",
+    "CurvedBoundaryLocalModelSample",
+    "DmkSplitExperimentReport",
+    "DmkSplitSample",
+    "FanTemplateMap2D",
+    "GraphBoundaryModelReport",
+    "GraphBoundaryModelSample",
+    "NearfieldTemplateExperiment",
+    "PrecomputedBoundaryTableReport",
+    "PrecomputedBoundaryTableSample",
+    "SeedQualitySample",
+    "TaylorBoundaryModelReport",
+    "TaylorBoundaryModelSample",
+    "diagonal_remainder_sample",
+    "ewald_log_local_kernel",
+    "ewald_log_local_radial_moment_2d",
+    "ewald_log_local_sector_monomial_moment_2d",
+    "ewald_log_smooth_kernel",
+    "expected_scaled_laplace_point_potential",
+    "exponential_integral_e1",
+    "laplace_log_kernel",
+    "metric_model_distance",
+    "point_target_laplace_potential",
+    "point_target_smooth_ewald_potential",
+    "run_curved_boundary_local_model_experiment",
+    "run_dmk_split_nearfield_experiment",
+    "run_graph_boundary_model_experiment",
+    "run_nearfield_template_experiment",
+    "run_precomputed_boundary_table_experiment",
+    "run_taylor_boundary_model_experiment",
+    "template_density_mass",
     "FormDslParityBenchmark",
     "FormDslParityRow",
     "FormDslMultipatchStressBenchmark",
diff --git a/src/cutkit/evals/nearfield_templates.py b/src/cutkit/evals/nearfield_templates.py
new file mode 100644
index 0000000..18e905f
--- /dev/null
+++ b/src/cutkit/evals/nearfield_templates.py
@@ -0,0 +1,1763 @@
+"""Near-field template experiments for folded 2D fan pieces."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from math import atan2, cos, erf, exp, gamma, hypot, isfinite, log, pi, sin, sqrt
+from typing import Callable
+
+from cutkit.geometry import (
+    CurveEdge2D,
+    CurveLoop2D,
+    CurveTrimmedPanel2D,
+    PanelLoop2D,
+    TrimmedPanel2D,
+)
+from cutkit.quadrature import folded_curve_quadrature_rule, gauss_legendre_01
+from cutkit.topology import select_interior_anchor
+
+Point2D = tuple[float, float]
+TemplateDensity2D = Callable[[float, float], float]
+_COINCIDENT_TOL = 1.0e-14
+_EULER_GAMMA = 0.5772156649015329
+
+
+def unit_template_density(_r: float, _t: float) -> float:
+    """Return the constant template density used by default experiments."""
+
+    return 1.0
+
+
+def _sub(lhs: Point2D, rhs: Point2D) -> Point2D:
+    return (lhs[0] - rhs[0], lhs[1] - rhs[1])
+
+
+def _dot(lhs: Point2D, rhs: Point2D) -> float:
+    return lhs[0] * rhs[0] + lhs[1] * rhs[1]
+
+
+def _cross(lhs: Point2D, rhs: Point2D) -> float:
+    return lhs[0] * rhs[1] - lhs[1] * rhs[0]
+
+
+def _scale_point(point: Point2D, scale: float) -> Point2D:
+    return (scale * point[0], scale * point[1])
+
+
+@dataclass(frozen=True)
+class FanTemplateMap2D:
+    """Straight-edge folded fan chart ``T(r, t) = (1-r)V + r C(t)``."""
+
+    vertex: Point2D
+    edge_start: Point2D
+    edge_end: Point2D
+
+    def point(self, r: float, t: float) -> Point2D:
+        """Map ``(r, t)`` in ``[0, 1]^2`` to physical space."""
+
+        one_minus_t = 1.0 - t
+        curve_x = one_minus_t * self.edge_start[0] + t * self.edge_end[0]
+        curve_y = one_minus_t * self.edge_start[1] + t * self.edge_end[1]
+        return (
+            (1.0 - r) * self.vertex[0] + r * curve_x,
+            (1.0 - r) * self.vertex[1] + r * curve_y,
+        )
+
+    @property
+    def boundary_det(self) -> float:
+        """Return ``det(C(t)-V, C'(t))``, constant for a straight edge."""
+
+        return _cross(
+            _sub(self.edge_start, self.vertex),
+            _sub(self.edge_end, self.edge_start),
+        )
+
+    @property
+    def signed_area(self) -> float:
+        """Return the oriented physical area of the fan piece."""
+
+        return 0.5 * self.boundary_det
+
+    def signed_jacobian(self, r: float) -> float:
+        """Return the signed Jacobian of the ``(r, t)`` fan chart."""
+
+        return r * self.boundary_det
+
+    def local_metric(
+        self,
+        r: float,
+        t: float,
+    ) -> tuple[tuple[float, float], tuple[float, float]]:
+        """Return ``DT(r,t)^T DT(r,t)`` for the fan chart."""
+
+        radial = _sub(self.point(1.0, t), self.vertex)
+        tangent_base = _sub(self.edge_end, self.edge_start)
+        tangent = (r * tangent_base[0], r * tangent_base[1])
+        return (
+            (_dot(radial, radial), _dot(radial, tangent)),
+            (_dot(tangent, radial), _dot(tangent, tangent)),
+        )
+
+    def scaled(self, factor: float) -> FanTemplateMap2D:
+        """Return the same template scaled about the origin."""
+
+        return FanTemplateMap2D(
+            vertex=_scale_point(self.vertex, factor),
+            edge_start=_scale_point(self.edge_start, factor),
+            edge_end=_scale_point(self.edge_end, factor),
+        )
+
+
+@dataclass(frozen=True)
+class DiagonalRemainderSample:
+    """One mapped-kernel singular split sample near the diagonal."""
+
+    delta: float
+    physical_distance: float
+    model_distance: float
+    laplace_kernel: float
+    metric_model_kernel: float
+    remainder: float
+
+
+@dataclass(frozen=True)
+class NearfieldTemplateExperiment:
+    """Summary of one fan-chart point-target experiment."""
+
+    fan: FanTemplateMap2D
+    order: int
+    near_point_target: Point2D
+    point_target_reference: float
+    point_target_low_order: float
+    point_target_abs_error: float
+    scaled_near_point_target: Point2D
+    scaled_point_target_potential: float
+    expected_scaled_point_target_potential: float
+    scaled_abs_error: float
+    signed_area: float
+    density_mass: float
+    scale_factor: float
+    diagonal_remainders: tuple[DiagonalRemainderSample, ...]
+
+
+@dataclass(frozen=True)
+class DmkSplitSample:
+    """One quadrature-convergence sample for the DMK-style log-kernel split."""
+
+    fixture: str
+    seed_mode: str
+    seed: Point2D
+    target_label: str
+    target: Point2D
+    sigma: float
+    order: int
+    full_value: float
+    smooth_value: float
+    local_value: float
+    analytic_local_value: float
+    reconstructed_value: float
+    analytic_reconstructed_value: float
+    full_reference: float
+    smooth_reference: float
+    local_reference: float
+    full_abs_error: float
+    smooth_abs_error: float
+    local_abs_error: float
+    analytic_local_abs_error: float
+    reconstructed_abs_error: float
+    analytic_reconstructed_abs_error: float
+    local_to_full_ratio: float
+
+
+@dataclass(frozen=True)
+class SeedQualitySample:
+    """Geometry quality summary for one folded-panel seed."""
+
+    fixture: str
+    seed_mode: str
+    seed: Point2D
+    min_abs_boundary_det: float
+    max_abs_boundary_det: float
+    max_to_min_abs_boundary_det: float
+    min_edge_distance: float
+    max_edge_distance: float
+    max_to_min_edge_distance: float
+
+
+@dataclass(frozen=True)
+class DmkSplitExperimentReport:
+    """Report payload for the DMK-style folded near-field split experiment."""
+
+    reference_order: int
+    orders: tuple[int, ...]
+    sigmas: tuple[float, ...]
+    seed_quality: tuple[SeedQualitySample, ...]
+    samples: tuple[DmkSplitSample, ...]
+
+
+@dataclass(frozen=True)
+class CurvedBoundaryLocalModelSample:
+    """Practical local-model check on one curved CAD-style cut cell."""
+
+    sigma: float
+    order: int
+    reference_order: int
+    target: Point2D
+    folded_value: float
+    reference_value: float
+    curved_boundary_value: float
+    half_plane_value: float
+    full_plane_value: float
+    folded_abs_error: float
+    curved_boundary_abs_error: float
+    half_plane_abs_error: float
+    full_plane_abs_error: float
+    curved_boundary_rel_error: float
+    half_plane_rel_error: float
+    full_plane_rel_error: float
+
+
+@dataclass(frozen=True)
+class CurvedBoundaryLocalModelReport:
+    """Report for local residual models on an exact curved cut cell."""
+
+    sigmas: tuple[float, ...]
+    orders: tuple[int, ...]
+    reference_order: int
+    samples: tuple[CurvedBoundaryLocalModelSample, ...]
+
+
+@dataclass(frozen=True)
+class TaylorBoundaryModelSample:
+    """One high-order Taylor boundary model sample for a smooth trim."""
+
+    sigma: float
+    boundary_order: int
+    exact_value: float
+    taylor_value: float
+    abs_error: float
+    rel_error: float
+
+
+@dataclass(frozen=True)
+class TaylorBoundaryModelReport:
+    """Convergence report for high-order smooth-boundary Taylor models."""
+
+    sigmas: tuple[float, ...]
+    boundary_orders: tuple[int, ...]
+    angular_order: int
+    samples: tuple[TaylorBoundaryModelSample, ...]
+
+
+@dataclass(frozen=True)
+class GraphBoundaryModelSample:
+    """One graph-strip smooth-boundary model sample."""
+
+    sigma: float
+    boundary_order: int
+    strip_order: int
+    exact_value: float
+    graph_value: float
+    abs_error: float
+    rel_error: float
+
+
+@dataclass(frozen=True)
+class GraphBoundaryModelReport:
+    """Convergence report for graph-strip boundary corrections."""
+
+    sigmas: tuple[float, ...]
+    boundary_orders: tuple[int, ...]
+    strip_orders: tuple[int, ...]
+    angular_order: int
+    samples: tuple[GraphBoundaryModelSample, ...]
+
+
+@dataclass(frozen=True)
+class PrecomputedBoundaryTableSample:
+    """One interpolated precomputed boundary-table sample."""
+
+    sigma: float
+    exact_value: float
+    table_value: float
+    abs_error: float
+    rel_error: float
+
+
+@dataclass(frozen=True)
+class PrecomputedBoundaryTableReport:
+    """Report for interpolation over precomputed smooth-boundary moments."""
+
+    table_sigmas: tuple[float, ...]
+    eval_sigmas: tuple[float, ...]
+    boundary_order: int
+    angular_order: int
+    samples: tuple[PrecomputedBoundaryTableSample, ...]
+
+
+def laplace_log_kernel(target: Point2D, source: Point2D) -> float:
+    """Return the 2D Laplace fundamental solution ``-log(|x-y|)/(2*pi)``."""
+
+    distance = hypot(target[0] - source[0], target[1] - source[1])
+    if distance <= 0.0:
+        raise ValueError("Laplace log kernel is singular at coincident points")
+    return -log(distance) / (2.0 * pi)
+
+
+def exponential_integral_e1(x: float, *, tol: float = 1.0e-15) -> float:
+    """Return ``E1(x) = integral_x^inf exp(-t)/t dt`` for ``x > 0``.
+
+    This small stdlib-only evaluator is intended for experiment-scale kernel
+    splitting, not as a general special-functions package.
+    """
+
+    if x <= 0.0 or not isfinite(x):
+        raise ValueError("x must be positive and finite")
+
+    if x <= 4.0:
+        term_power = 1.0
+        factorial = 1.0
+        series = 0.0
+        for k in range(1, 256):
+            term_power *= -x
+            factorial *= k
+            term = term_power / (k * factorial)
+            series += term
+            if abs(term) <= tol * max(1.0, abs(series)):
+                break
+        return -_EULER_GAMMA - log(x) - series
+
+    term = 1.0
+    series = 1.0
+    previous_abs = abs(term)
+    for k in range(1, 256):
+        term *= -k / x
+        current_abs = abs(term)
+        if current_abs > previous_abs:
+            break
+        series += term
+        if current_abs <= tol * max(1.0, abs(series)):
+            break
+        previous_abs = current_abs
+    return exp(-x) * series / x
+
+
+def ewald_log_local_kernel(target: Point2D, source: Point2D, *, sigma: float) -> float:
+    """Return the localized 2D log Ewald residual ``E1(r^2/sigma^2)/(4*pi)``."""
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    distance = hypot(target[0] - source[0], target[1] - source[1])
+    if distance <= 0.0:
+        raise ValueError("local log kernel is singular at coincident points")
+    return exponential_integral_e1((distance / sigma) ** 2) / (4.0 * pi)
+
+
+def ewald_log_smooth_kernel(target: Point2D, source: Point2D, *, sigma: float) -> float:
+    """Return the smooth complement of the localized 2D log Ewald residual."""
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    distance = hypot(target[0] - source[0], target[1] - source[1])
+    if distance <= _COINCIDENT_TOL:
+        return (_EULER_GAMMA - 2.0 * log(sigma)) / (4.0 * pi)
+    return laplace_log_kernel(target, source) - ewald_log_local_kernel(
+        target,
+        source,
+        sigma=sigma,
+    )
+
+
+def ewald_log_local_radial_moment_2d(*, sigma: float, monomial_degree: int) -> float:
+    """Return ``int_0^inf K_sigma(r) r^(monomial_degree + 1) dr``.
+
+    Here ``K_sigma(r) = E1(r^2/sigma^2)/(4*pi)`` is the localized 2D log
+    residual. Multiplying this radial moment by the corresponding angular
+    monomial moment gives sector moments for full-plane, half-plane, and wedge
+    local models.
+    """
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    if monomial_degree < 0:
+        raise ValueError("monomial_degree must be nonnegative")
+
+    half_degree = 0.5 * monomial_degree
+    return (
+        sigma ** (monomial_degree + 2)
+        * gamma(half_degree + 1.0)
+        / (8.0 * pi * (half_degree + 1.0))
+    )
+
+
+def ewald_log_local_sector_monomial_moment_2d(
+    *,
+    sigma: float,
+    x_power: int,
+    y_power: int,
+    theta_start: float,
+    theta_end: float,
+    angular_order: int = 96,
+) -> float:
+    """Return a 2D local residual monomial moment over an angular sector.
+
+    The sector has vertex at the target and angle range ``[theta_start,
+    theta_end]``. This is the leading straight-feature model for full-plane,
+    half-plane, and wedge residual corrections.
+    """
+
+    if x_power < 0 or y_power < 0:
+        raise ValueError("monomial powers must be nonnegative")
+    if angular_order < 1:
+        raise ValueError("angular_order must be positive")
+    if theta_end < theta_start:
+        raise ValueError("theta_end must be at least theta_start")
+
+    radial = ewald_log_local_radial_moment_2d(
+        sigma=sigma,
+        monomial_degree=x_power + y_power,
+    )
+    if theta_end == theta_start:
+        return 0.0
+
+    nodes, weights = gauss_legendre_01(angular_order)
+    width = theta_end - theta_start
+    angular = 0.0
+    for node, weight in zip(nodes, weights, strict=True):
+        theta = theta_start + width * node
+        angular += weight * width * (cos(theta) ** x_power) * (sin(theta) ** y_power)
+    return radial * angular
+
+
+def metric_model_distance(
+    fan: FanTemplateMap2D,
+    *,
+    r: float,
+    t: float,
+    delta_r: float,
+    delta_t: float,
+) -> float:
+    """Return the local metric distance for a template-space displacement."""
+
+    metric = fan.local_metric(r, t)
+    distance_squared = (
+        metric[0][0] * delta_r * delta_r
+        + 2.0 * metric[0][1] * delta_r * delta_t
+        + metric[1][1] * delta_t * delta_t
+    )
+    if distance_squared <= 0.0:
+        raise ValueError("metric model distance must be positive")
+    return sqrt(distance_squared)
+
+
+def template_density_mass(
+    fan: FanTemplateMap2D,
+    *,
+    order: int,
+    density: TemplateDensity2D | None = None,
+) -> float:
+    """Return ``integral density(r,t) * J(r,t) dr dt`` on one fan chart."""
+
+    if order < 1:
+        raise ValueError("order must be positive")
+    if density is None:
+        density = unit_template_density
+
+    nodes, weights = gauss_legendre_01(order)
+    total = 0.0
+    for r, wr in zip(nodes, weights, strict=True):
+        jacobian = fan.signed_jacobian(r)
+        for t, wt in zip(nodes, weights, strict=True):
+            total += density(r, t) * jacobian * wr * wt
+    return total
+
+
+def point_target_laplace_potential(
+    fan: FanTemplateMap2D,
+    target: Point2D,
+    *,
+    order: int,
+    source_density: TemplateDensity2D | None = None,
+) -> float:
+    """Approximate a folded-piece source integral at one point target."""
+
+    if order < 1:
+        raise ValueError("order must be positive")
+    if source_density is None:
+        source_density = unit_template_density
+
+    nodes, weights = gauss_legendre_01(order)
+    total = 0.0
+    for r, wr in zip(nodes, weights, strict=True):
+        jacobian = fan.signed_jacobian(r)
+        for t, wt in zip(nodes, weights, strict=True):
+            source = fan.point(r, t)
+            if hypot(target[0] - source[0], target[1] - source[1]) <= _COINCIDENT_TOL:
+                raise ValueError(
+                    "point target coincides with a source quadrature node; "
+                    "use a singular point-target reference rule"
+                )
+            total += (
+                laplace_log_kernel(target, source)
+                * source_density(r, t)
+                * jacobian
+                * wr
+                * wt
+            )
+    return total
+
+
+def point_target_smooth_ewald_potential(
+    fan: FanTemplateMap2D,
+    target: Point2D,
+    *,
+    order: int,
+    sigma: float,
+    source_density: TemplateDensity2D | None = None,
+) -> float:
+    """Approximate the smoothed 2D log-kernel source integral at one target."""
+
+    if order < 1:
+        raise ValueError("order must be positive")
+    if source_density is None:
+        source_density = unit_template_density
+
+    nodes, weights = gauss_legendre_01(order)
+    total = 0.0
+    for r, wr in zip(nodes, weights, strict=True):
+        jacobian = fan.signed_jacobian(r)
+        for t, wt in zip(nodes, weights, strict=True):
+            source = fan.point(r, t)
+            total += (
+                ewald_log_smooth_kernel(target, source, sigma=sigma)
+                * source_density(r, t)
+                * jacobian
+                * wr
+                * wt
+            )
+    return total
+
+
+def _edge_distance(point: Point2D, start: Point2D, end: Point2D) -> float:
+    edge = _sub(end, start)
+    length_sq = _dot(edge, edge)
+    if length_sq <= 0.0:
+        return hypot(point[0] - start[0], point[1] - start[1])
+    rel = _sub(point, start)
+    tau = max(0.0, min(1.0, _dot(rel, edge) / length_sq))
+    closest = (start[0] + tau * edge[0], start[1] + tau * edge[1])
+    return hypot(point[0] - closest[0], point[1] - closest[1])
+
+
+def _node_barycenter(nodes: tuple[Point2D, ...]) -> Point2D:
+    scale = 1.0 / len(nodes)
+    return (
+        scale * sum(point[0] for point in nodes),
+        scale * sum(point[1] for point in nodes),
+    )
+
+
+def _panel_fans(
+    nodes: tuple[Point2D, ...], *, seed: Point2D
+) -> tuple[FanTemplateMap2D, ...]:
+    return tuple(
+        FanTemplateMap2D(
+            vertex=seed,
+            edge_start=nodes[index],
+            edge_end=nodes[(index + 1) % len(nodes)],
+        )
+        for index in range(len(nodes))
+    )
+
+
+def _panel_seed_quality(
+    *,
+    fixture: str,
+    seed_mode: str,
+    nodes: tuple[Point2D, ...],
+    seed: Point2D,
+) -> SeedQualitySample:
+    fans = _panel_fans(nodes, seed=seed)
+    dets = tuple(abs(fan.boundary_det) for fan in fans)
+    distances = tuple(
+        _edge_distance(seed, nodes[index], nodes[(index + 1) % len(nodes)])
+        for index in range(len(nodes))
+    )
+    min_det = min(dets)
+    min_distance = min(distances)
+    max_det = max(dets)
+    max_distance = max(distances)
+    return SeedQualitySample(
+        fixture=fixture,
+        seed_mode=seed_mode,
+        seed=seed,
+        min_abs_boundary_det=min_det,
+        max_abs_boundary_det=max_det,
+        max_to_min_abs_boundary_det=max_det / min_det
+        if min_det > 0.0
+        else float("inf"),
+        min_edge_distance=min_distance,
+        max_edge_distance=max_distance,
+        max_to_min_edge_distance=(
+            max_distance / min_distance if min_distance > 0.0 else float("inf")
+        ),
+    )
+
+
+def _physical_density(point: Point2D) -> float:
+    return 1.0 + 0.25 * point[0] - 0.15 * point[1] + 0.1 * point[0] * point[1]
+
+
+def _rational_quarter_circle_edge() -> CurveEdge2D:
+    """Return the exact rational quadratic arc used as a curved cut trim."""
+
+    p0 = (0.0, 0.5)
+    p1 = (0.0, 0.0)
+    p2 = (0.5, 0.0)
+    w0 = 1.0
+    w1 = 1.0 / sqrt(2.0)
+    w2 = 1.0
+
+    def basis(t: float) -> tuple[float, float, float]:
+        one_minus_t = 1.0 - t
+        return (one_minus_t * one_minus_t, 2.0 * one_minus_t * t, t * t)
+
+    def basis_derivative(t: float) -> tuple[float, float, float]:
+        return (-2.0 * (1.0 - t), 2.0 - 4.0 * t, 2.0 * t)
+
+    def weighted_sum(t: float) -> tuple[float, float, float]:
+        b0, b1, b2 = basis(t)
+        denominator = w0 * b0 + w1 * b1 + w2 * b2
+        numerator_x = w0 * b0 * p0[0] + w1 * b1 * p1[0] + w2 * b2 * p2[0]
+        numerator_y = w0 * b0 * p0[1] + w1 * b1 * p1[1] + w2 * b2 * p2[1]
+        return numerator_x, numerator_y, denominator
+
+    def weighted_sum_derivative(t: float) -> tuple[float, float, float]:
+        db0, db1, db2 = basis_derivative(t)
+        denominator_derivative = w0 * db0 + w1 * db1 + w2 * db2
+        numerator_x_derivative = w0 * db0 * p0[0] + w1 * db1 * p1[0] + w2 * db2 * p2[0]
+        numerator_y_derivative = w0 * db0 * p0[1] + w1 * db1 * p1[1] + w2 * db2 * p2[1]
+        return numerator_x_derivative, numerator_y_derivative, denominator_derivative
+
+    def point(t: float) -> Point2D:
+        numerator_x, numerator_y, denominator = weighted_sum(t)
+        return (numerator_x / denominator, numerator_y / denominator)
+
+    def tangent(t: float) -> Point2D:
+        numerator_x, numerator_y, denominator = weighted_sum(t)
+        dx_num, dy_num, denominator_derivative = weighted_sum_derivative(t)
+        denominator_sq = denominator * denominator
+        return (
+            (dx_num * denominator - numerator_x * denominator_derivative)
+            / denominator_sq,
+            (dy_num * denominator - numerator_y * denominator_derivative)
+            / denominator_sq,
+        )
+
+    return CurveEdge2D(evaluator=point, derivative=tangent)
+
+
+def _rational_quarter_circle_cut_cell() -> CurveTrimmedPanel2D:
+    """Return a curved CAD-style cut cell with one rational trim edge."""
+
+    curved = _rational_quarter_circle_edge()
+    p0 = curved.point(0.0)
+    p1 = curved.point(1.0)
+    cell_corner = (0.5, 0.5)
+    loop = CurveLoop2D(
+        (
+            curved,
+            CurveEdge2D.line(p1, cell_corner),
+            CurveEdge2D.line(cell_corner, p0),
+        )
+    )
+    return CurveTrimmedPanel2D(outer=loop)
+
+
+def _curve_panel_local_residual(
+    panel: CurveTrimmedPanel2D,
+    target: Point2D,
+    *,
+    sigma: float,
+    order: int,
+) -> float:
+    if order < 1:
+        raise ValueError("order must be positive")
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+
+    folded = folded_curve_quadrature_rule(panel, order=order)
+    total = 0.0
+    for source, weight in zip(folded.rule.points, folded.rule.weights, strict=True):
+        total += (
+            ewald_log_local_kernel(target, source, sigma=sigma)
+            * _physical_density(source)
+            * weight
+        )
+    return total
+
+
+def _quarter_circle_ray_exit_distance(target: Point2D, direction: Point2D) -> float:
+    """Return ray exit length for the exact rational quarter-circle cut cell."""
+
+    center = (0.5, 0.5)
+    rel = _sub(target, center)
+    circle_dot = _dot(rel, direction)
+    if circle_dot >= 0.0:
+        return 0.0
+
+    exit_distance = -2.0 * circle_dot
+    dx, dy = direction
+    if dx > 0.0:
+        exit_distance = min(exit_distance, (0.5 - target[0]) / dx)
+    if dy > 0.0:
+        exit_distance = min(exit_distance, (0.5 - target[1]) / dy)
+    return max(0.0, exit_distance)
+
+
+def _quarter_circle_exact_boundary_local_residual(
+    target: Point2D,
+    *,
+    sigma: float,
+    angular_order: int,
+) -> float:
+    """Integrate the local residual using exact target-centered ray clipping."""
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    if angular_order < 1:
+        raise ValueError("angular_order must be positive")
+
+    p0 = (0.0, 0.5)
+    p1 = (0.5, 0.0)
+    center = (0.5, 0.5)
+    tangent_low = -0.25 * pi
+    tangent_high = 0.75 * pi
+    intervals = (
+        (
+            tangent_low,
+            atan2(p1[1] - target[1], p1[0] - target[0]),
+            "circle_start",
+        ),
+        (
+            atan2(p1[1] - target[1], p1[0] - target[0]),
+            atan2(center[1] - target[1], center[0] - target[0]),
+            "right",
+        ),
+        (
+            atan2(center[1] - target[1], center[0] - target[0]),
+            atan2(p0[1] - target[1], p0[0] - target[0]),
+            "top",
+        ),
+        (
+            atan2(p0[1] - target[1], p0[0] - target[0]),
+            tangent_high,
+            "circle_end",
+        ),
+    )
+    total = 0.0
+    for theta_start, theta_end, exit_kind in intervals:
+
+        def integrand(theta: float, *, kind: str = exit_kind) -> float:
+            return _quarter_circle_ray_integrand(
+                target,
+                theta,
+                sigma=sigma,
+                exit_kind=kind,
+            )
+
+        total += _adaptive_simpson(
+            integrand,
+            theta_start,
+            theta_end,
+            abs_tol=1.0e-16,
+            max_depth=max(20, angular_order),
+        )
+    return total
+
+
+def _quarter_circle_ray_integrand(
+    target: Point2D,
+    theta: float,
+    *,
+    sigma: float,
+    exit_kind: str,
+) -> float:
+    direction = (cos(theta), sin(theta))
+    if exit_kind in ("circle_start", "circle_end"):
+        exit_distance = _quarter_circle_ray_exit_distance(target, direction)
+    elif exit_kind == "right":
+        exit_distance = (0.5 - target[0]) / direction[0]
+    else:
+        exit_distance = (0.5 - target[1]) / direction[1]
+    if exit_distance <= 0.0:
+        return 0.0
+    density_constant = _physical_density(target)
+    density_linear = (
+        0.25 * direction[0]
+        - 0.15 * direction[1]
+        + 0.1 * (target[0] * direction[1] + target[1] * direction[0])
+    )
+    density_quadratic = 0.1 * direction[0] * direction[1]
+    return (
+        density_constant
+        * _ewald_log_local_finite_radial_moment_2d(
+            sigma=sigma,
+            monomial_degree=0,
+            radius=exit_distance,
+        )
+        + density_linear
+        * _ewald_log_local_finite_radial_moment_2d(
+            sigma=sigma,
+            monomial_degree=1,
+            radius=exit_distance,
+        )
+        + density_quadratic
+        * _ewald_log_local_finite_radial_moment_2d(
+            sigma=sigma,
+            monomial_degree=2,
+            radius=exit_distance,
+        )
+    )
+
+
+def _adaptive_simpson(
+    function: Callable[[float], float],
+    start: float,
+    end: float,
+    *,
+    abs_tol: float,
+    max_depth: int,
+) -> float:
+    midpoint = 0.5 * (start + end)
+    f_start = function(start)
+    f_mid = function(midpoint)
+    f_end = function(end)
+    whole = _simpson_estimate(start, end, f_start, f_mid, f_end)
+    return _adaptive_simpson_recursive(
+        function,
+        start,
+        midpoint,
+        end,
+        f_start,
+        f_mid,
+        f_end,
+        whole,
+        abs_tol,
+        max_depth,
+    )
+
+
+def _simpson_estimate(
+    start: float,
+    end: float,
+    f_start: float,
+    f_mid: float,
+    f_end: float,
+) -> float:
+    return (end - start) * (f_start + 4.0 * f_mid + f_end) / 6.0
+
+
+def _adaptive_simpson_recursive(
+    function: Callable[[float], float],
+    start: float,
+    midpoint: float,
+    end: float,
+    f_start: float,
+    f_mid: float,
+    f_end: float,
+    whole: float,
+    abs_tol: float,
+    depth_remaining: int,
+) -> float:
+    left_mid = 0.5 * (start + midpoint)
+    right_mid = 0.5 * (midpoint + end)
+    f_left_mid = function(left_mid)
+    f_right_mid = function(right_mid)
+    left = _simpson_estimate(start, midpoint, f_start, f_left_mid, f_mid)
+    right = _simpson_estimate(midpoint, end, f_mid, f_right_mid, f_end)
+    delta = left + right - whole
+    if depth_remaining <= 0 or abs(delta) <= 15.0 * abs_tol:
+        return left + right + delta / 15.0
+    return _adaptive_simpson_recursive(
+        function,
+        start,
+        left_mid,
+        midpoint,
+        f_start,
+        f_left_mid,
+        f_mid,
+        left,
+        0.5 * abs_tol,
+        depth_remaining - 1,
+    ) + _adaptive_simpson_recursive(
+        function,
+        midpoint,
+        right_mid,
+        end,
+        f_mid,
+        f_right_mid,
+        f_end,
+        right,
+        0.5 * abs_tol,
+        depth_remaining - 1,
+    )
+
+
+def _ewald_log_local_finite_radial_moment_2d(
+    *,
+    sigma: float,
+    monomial_degree: int,
+    radius: float,
+) -> float:
+    """Return ``int_0^radius K_sigma(r) r^(monomial_degree + 1) dr``."""
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    if radius <= 0.0:
+        return 0.0
+    if monomial_degree not in (0, 1, 2):
+        raise ValueError("finite radial moment supports degrees 0, 1, and 2")
+
+    upper = (radius / sigma) ** 2
+    sqrt_upper = sqrt(upper)
+    exp_upper = exp(-upper) if upper < 745.0 else 0.0
+    e1_upper = exponential_integral_e1(upper)
+    if monomial_degree == 0:
+        lower_gamma = 1.0 - exp_upper
+        q_plus_one = 1.0
+    elif monomial_degree == 1:
+        lower_gamma = 0.5 * sqrt(pi) * erf(sqrt_upper) - sqrt_upper * exp_upper
+        q_plus_one = 1.5
+    else:
+        lower_gamma = 1.0 - exp_upper * (1.0 + upper)
+        q_plus_one = 2.0
+
+    integral_u = (upper**q_plus_one * e1_upper + lower_gamma) / q_plus_one
+    return sigma ** (monomial_degree + 2) * integral_u / (8.0 * pi)
+
+
+def _quarter_circle_boundary_frame() -> tuple[Point2D, Point2D, Point2D, float]:
+    panel = _rational_quarter_circle_cut_cell()
+    target = panel.outer.edges[0].point(0.5)
+    inv_sqrt2 = 1.0 / sqrt(2.0)
+    tangent = (inv_sqrt2, -inv_sqrt2)
+    inward_normal = (inv_sqrt2, inv_sqrt2)
+    return target, tangent, inward_normal, 0.5
+
+
+def _quarter_circle_taylor_graph_coefficients(
+    *, radius: float, boundary_order: int
+) -> tuple[float, ...]:
+    if boundary_order < 2 or boundary_order % 2 != 0:
+        raise ValueError("boundary_order must be a positive even order at least 2")
+
+    coefficients: list[float] = []
+    binomial = 1.0
+    for power in range(1, boundary_order // 2 + 1):
+        binomial *= (1.5 - float(power)) / float(power)
+        coefficient = -binomial * ((-1.0) ** power) / (radius ** (2 * power - 1))
+        coefficients.append(coefficient)
+    return tuple(coefficients)
+
+
+def _quarter_circle_taylor_graph_value(
+    s_coord: float,
+    *,
+    coefficients: tuple[float, ...],
+) -> float:
+    value = 0.0
+    s_squared = s_coord * s_coord
+    power = s_squared
+    for coefficient in coefficients:
+        value += coefficient * power
+        power *= s_squared
+    return value
+
+
+def _taylor_graph_exit_distance(
+    theta: float,
+    *,
+    coefficients: tuple[float, ...],
+    sigma: float,
+) -> float | None:
+    cos_theta = cos(theta)
+    sin_theta = sin(theta)
+    if sin_theta <= 0.0:
+        return 0.0
+    if abs(cos_theta) <= 1.0e-14:
+        return None
+
+    def boundary_residual(radius: float) -> float:
+        return radius * sin_theta - _quarter_circle_taylor_graph_value(
+            radius * cos_theta,
+            coefficients=coefficients,
+        )
+
+    low = 0.0
+    high = max(sigma, 1.0e-3)
+    while boundary_residual(high) > 0.0:
+        high *= 2.0
+        if high > 8.0:
+            return None
+
+    for _ in range(96):
+        mid = 0.5 * (low + high)
+        if boundary_residual(mid) > 0.0:
+            low = mid
+        else:
+            high = mid
+    return high
+
+
+def _local_boundary_ray_integral(
+    target: Point2D,
+    direction: Point2D,
+    *,
+    sigma: float,
+    exit_distance: float | None,
+) -> float:
+    density_constant = _physical_density(target)
+    density_linear = (
+        0.25 * direction[0]
+        - 0.15 * direction[1]
+        + 0.1 * (target[0] * direction[1] + target[1] * direction[0])
+    )
+    density_quadratic = 0.1 * direction[0] * direction[1]
+
+    def radial_moment(degree: int) -> float:
+        if exit_distance is None:
+            return ewald_log_local_radial_moment_2d(
+                sigma=sigma,
+                monomial_degree=degree,
+            )
+        return _ewald_log_local_finite_radial_moment_2d(
+            sigma=sigma,
+            monomial_degree=degree,
+            radius=exit_distance,
+        )
+
+    return (
+        density_constant * radial_moment(0)
+        + density_linear * radial_moment(1)
+        + density_quadratic * radial_moment(2)
+    )
+
+
+def _quarter_circle_exact_local_boundary_model(
+    *, sigma: float, angular_order: int
+) -> float:
+    target, tangent, inward_normal, radius = _quarter_circle_boundary_frame()
+    if angular_order < 1:
+        raise ValueError("angular_order must be positive")
+
+    def integrand(theta: float) -> float:
+        direction = (
+            cos(theta) * tangent[0] + sin(theta) * inward_normal[0],
+            cos(theta) * tangent[1] + sin(theta) * inward_normal[1],
+        )
+        exit_distance = 2.0 * radius * sin(theta)
+        return _local_boundary_ray_integral(
+            target,
+            direction,
+            sigma=sigma,
+            exit_distance=exit_distance,
+        )
+
+    return _adaptive_simpson(
+        integrand,
+        0.0,
+        pi,
+        abs_tol=1.0e-16,
+        max_depth=max(20, angular_order),
+    )
+
+
+def _quarter_circle_taylor_boundary_model(
+    *, sigma: float, boundary_order: int, angular_order: int
+) -> float:
+    target, tangent, inward_normal, radius = _quarter_circle_boundary_frame()
+    coefficients = _quarter_circle_taylor_graph_coefficients(
+        radius=radius,
+        boundary_order=boundary_order,
+    )
+
+    def integrand(theta: float) -> float:
+        direction = (
+            cos(theta) * tangent[0] + sin(theta) * inward_normal[0],
+            cos(theta) * tangent[1] + sin(theta) * inward_normal[1],
+        )
+        exit_distance = _taylor_graph_exit_distance(
+            theta,
+            coefficients=coefficients,
+            sigma=sigma,
+        )
+        return _local_boundary_ray_integral(
+            target,
+            direction,
+            sigma=sigma,
+            exit_distance=exit_distance,
+        )
+
+    return _adaptive_simpson(
+        integrand,
+        0.0,
+        pi,
+        abs_tol=1.0e-16,
+        max_depth=max(20, angular_order),
+    )
+
+
+def run_taylor_boundary_model_experiment(
+    *,
+    sigmas: tuple[float, ...] = (0.01, 0.02),
+    boundary_orders: tuple[int, ...] = (2, 4, 6, 8, 10, 12),
+    angular_order: int = 32,
+) -> TaylorBoundaryModelReport:
+    """Compare high-order local boundary Taylor models to exact circular geometry."""
+
+    samples: list[TaylorBoundaryModelSample] = []
+    for sigma in sigmas:
+        exact_value = _quarter_circle_exact_local_boundary_model(
+            sigma=sigma,
+            angular_order=angular_order,
+        )
+        denominator = max(abs(exact_value), 1.0e-300)
+        for boundary_order in boundary_orders:
+            taylor_value = _quarter_circle_taylor_boundary_model(
+                sigma=sigma,
+                boundary_order=boundary_order,
+                angular_order=angular_order,
+            )
+            abs_error = abs(taylor_value - exact_value)
+            samples.append(
+                TaylorBoundaryModelSample(
+                    sigma=sigma,
+                    boundary_order=boundary_order,
+                    exact_value=exact_value,
+                    taylor_value=taylor_value,
+                    abs_error=abs_error,
+                    rel_error=abs_error / denominator,
+                )
+            )
+    return TaylorBoundaryModelReport(
+        sigmas=sigmas,
+        boundary_orders=boundary_orders,
+        angular_order=angular_order,
+        samples=tuple(samples),
+    )
+
+
+def _local_half_plane_boundary_model(*, sigma: float, angular_order: int) -> float:
+    target, tangent, inward_normal, _radius = _quarter_circle_boundary_frame()
+
+    def integrand(theta: float) -> float:
+        direction = (
+            cos(theta) * tangent[0] + sin(theta) * inward_normal[0],
+            cos(theta) * tangent[1] + sin(theta) * inward_normal[1],
+        )
+        return _local_boundary_ray_integral(
+            target,
+            direction,
+            sigma=sigma,
+            exit_distance=None,
+        )
+
+    return _adaptive_simpson(
+        integrand,
+        0.0,
+        pi,
+        abs_tol=1.0e-16,
+        max_depth=max(20, angular_order),
+    )
+
+
+def _graph_strip_correction(
+    *,
+    sigma: float,
+    boundary_order: int,
+    strip_order: int,
+) -> float:
+    target, tangent, inward_normal, radius = _quarter_circle_boundary_frame()
+    coefficients = _quarter_circle_taylor_graph_coefficients(
+        radius=radius,
+        boundary_order=boundary_order,
+    )
+    nodes, weights = gauss_legendre_01(strip_order)
+    total = 0.0
+    for s_node, s_weight in zip(nodes, weights, strict=True):
+        s_coord = radius * (2.0 * s_node - 1.0)
+        ds_weight = 2.0 * radius * s_weight
+        graph_height = _quarter_circle_taylor_graph_value(
+            s_coord,
+            coefficients=coefficients,
+        )
+        if graph_height <= 0.0:
+            continue
+        for n_node, n_weight in zip(nodes, weights, strict=True):
+            n_coord = graph_height * n_node
+            source = (
+                target[0] + s_coord * tangent[0] + n_coord * inward_normal[0],
+                target[1] + s_coord * tangent[1] + n_coord * inward_normal[1],
+            )
+            if hypot(source[0] - target[0], source[1] - target[1]) <= _COINCIDENT_TOL:
+                continue
+            total += (
+                ewald_log_local_kernel(target, source, sigma=sigma)
+                * _physical_density(source)
+                * ds_weight
+                * graph_height
+                * n_weight
+            )
+    return total
+
+
+def _quarter_circle_graph_boundary_model(
+    *,
+    sigma: float,
+    boundary_order: int,
+    strip_order: int,
+    angular_order: int,
+) -> float:
+    return _local_half_plane_boundary_model(
+        sigma=sigma,
+        angular_order=angular_order,
+    ) - _graph_strip_correction(
+        sigma=sigma,
+        boundary_order=boundary_order,
+        strip_order=strip_order,
+    )
+
+
+def run_graph_boundary_model_experiment(
+    *,
+    sigmas: tuple[float, ...] = (0.005, 0.01, 0.02),
+    boundary_orders: tuple[int, ...] = (4, 6, 8),
+    strip_orders: tuple[int, ...] = (32, 64, 96),
+    angular_order: int = 32,
+) -> GraphBoundaryModelReport:
+    """Compare graph-strip boundary corrections against exact circular geometry."""
+
+    samples: list[GraphBoundaryModelSample] = []
+    for sigma in sigmas:
+        exact_value = _quarter_circle_exact_local_boundary_model(
+            sigma=sigma,
+            angular_order=angular_order,
+        )
+        denominator = max(abs(exact_value), 1.0e-300)
+        for boundary_order in boundary_orders:
+            for strip_order in strip_orders:
+                graph_value = _quarter_circle_graph_boundary_model(
+                    sigma=sigma,
+                    boundary_order=boundary_order,
+                    strip_order=strip_order,
+                    angular_order=angular_order,
+                )
+                abs_error = abs(graph_value - exact_value)
+                samples.append(
+                    GraphBoundaryModelSample(
+                        sigma=sigma,
+                        boundary_order=boundary_order,
+                        strip_order=strip_order,
+                        exact_value=exact_value,
+                        graph_value=graph_value,
+                        abs_error=abs_error,
+                        rel_error=abs_error / denominator,
+                    )
+                )
+    return GraphBoundaryModelReport(
+        sigmas=sigmas,
+        boundary_orders=boundary_orders,
+        strip_orders=strip_orders,
+        angular_order=angular_order,
+        samples=tuple(samples),
+    )
+
+
+def _lagrange_interpolate(
+    nodes: tuple[float, ...], values: tuple[float, ...], point: float
+) -> float:
+    if len(nodes) != len(values):
+        raise ValueError("nodes and values must have the same length")
+    total = 0.0
+    for i, node_i in enumerate(nodes):
+        basis = 1.0
+        for j, node_j in enumerate(nodes):
+            if i == j:
+                continue
+            basis *= (point - node_j) / (node_i - node_j)
+        total += values[i] * basis
+    return total
+
+
+def run_precomputed_boundary_table_experiment(
+    *,
+    table_sigmas: tuple[float, ...] = (
+        0.004,
+        0.008,
+        0.012,
+        0.016,
+        0.020,
+        0.024,
+        0.028,
+    ),
+    eval_sigmas: tuple[float, ...] = (0.010, 0.018, 0.026),
+    boundary_order: int = 8,
+    angular_order: int = 32,
+) -> PrecomputedBoundaryTableReport:
+    """Interpolate precomputed high-order smooth-boundary moment values."""
+
+    table_values = tuple(
+        _quarter_circle_taylor_boundary_model(
+            sigma=sigma,
+            boundary_order=boundary_order,
+            angular_order=angular_order,
+        )
+        for sigma in table_sigmas
+    )
+    samples: list[PrecomputedBoundaryTableSample] = []
+    for sigma in eval_sigmas:
+        exact_value = _quarter_circle_taylor_boundary_model(
+            sigma=sigma,
+            boundary_order=boundary_order,
+            angular_order=angular_order,
+        )
+        table_value = _lagrange_interpolate(table_sigmas, table_values, sigma)
+        abs_error = abs(table_value - exact_value)
+        samples.append(
+            PrecomputedBoundaryTableSample(
+                sigma=sigma,
+                exact_value=exact_value,
+                table_value=table_value,
+                abs_error=abs_error,
+                rel_error=abs_error / max(abs(exact_value), 1.0e-300),
+            )
+        )
+    return PrecomputedBoundaryTableReport(
+        table_sigmas=table_sigmas,
+        eval_sigmas=eval_sigmas,
+        boundary_order=boundary_order,
+        angular_order=angular_order,
+        samples=tuple(samples),
+    )
+
+
+def run_curved_boundary_local_model_experiment(
+    *,
+    sigmas: tuple[float, ...] = (0.02, 0.04, 0.08),
+    orders: tuple[int, ...] = (12, 20, 32),
+    reference_order: int = 96,
+) -> CurvedBoundaryLocalModelReport:
+    """Test local residual models on an exact curved rational cut cell."""
+
+    panel = _rational_quarter_circle_cut_cell()
+    trim = panel.outer.edges[0]
+    target = trim.point(0.5)
+
+    samples: list[CurvedBoundaryLocalModelSample] = []
+    for sigma in sigmas:
+        reference = _curve_panel_local_residual(
+            panel,
+            target,
+            sigma=sigma,
+            order=reference_order,
+        )
+        curved_boundary = _quarter_circle_exact_boundary_local_residual(
+            target,
+            sigma=sigma,
+            angular_order=reference_order,
+        )
+        density_at_target = _physical_density(target)
+        half_plane = (
+            ewald_log_local_sector_monomial_moment_2d(
+                sigma=sigma,
+                x_power=0,
+                y_power=0,
+                theta_start=0.0,
+                theta_end=pi,
+            )
+            * density_at_target
+        )
+        full_plane = (
+            ewald_log_local_sector_monomial_moment_2d(
+                sigma=sigma,
+                x_power=0,
+                y_power=0,
+                theta_start=0.0,
+                theta_end=2.0 * pi,
+            )
+            * density_at_target
+        )
+        denominator = max(abs(reference), 1.0e-300)
+        for order in orders:
+            folded_value = _curve_panel_local_residual(
+                panel,
+                target,
+                sigma=sigma,
+                order=order,
+            )
+            folded_error = abs(folded_value - reference)
+            curved_error = abs(curved_boundary - reference)
+            half_error = abs(half_plane - reference)
+            full_error = abs(full_plane - reference)
+            samples.append(
+                CurvedBoundaryLocalModelSample(
+                    sigma=sigma,
+                    order=order,
+                    reference_order=reference_order,
+                    target=target,
+                    folded_value=folded_value,
+                    reference_value=reference,
+                    curved_boundary_value=curved_boundary,
+                    half_plane_value=half_plane,
+                    full_plane_value=full_plane,
+                    folded_abs_error=folded_error,
+                    curved_boundary_abs_error=curved_error,
+                    half_plane_abs_error=half_error,
+                    full_plane_abs_error=full_error,
+                    curved_boundary_rel_error=curved_error / denominator,
+                    half_plane_rel_error=half_error / denominator,
+                    full_plane_rel_error=full_error / denominator,
+                )
+            )
+
+    return CurvedBoundaryLocalModelReport(
+        sigmas=sigmas,
+        orders=orders,
+        reference_order=reference_order,
+        samples=tuple(samples),
+    )
+
+
+def _analytic_full_space_local_residual(target: Point2D, *, sigma: float) -> float:
+    """Return the leading full-space Taylor moment for the 2D local residual.
+
+    For ``G_local = E1(r^2/sigma^2)/(4*pi)``, the zeroth full-space moment is
+    ``sigma^2/4``. The current experiment density has zero Laplacian, so the
+    second-order isotropic Taylor correction vanishes.
+    """
+
+    if sigma <= 0.0:
+        raise ValueError("sigma must be positive")
+    return (sigma * sigma / 4.0) * _physical_density(target)
+
+
+def _panel_potential(
+    fans: tuple[FanTemplateMap2D, ...],
+    target: Point2D,
+    *,
+    order: int,
+    sigma: float | None,
+    smooth_only: bool,
+    local_only: bool = False,
+) -> float:
+    if smooth_only and local_only:
+        raise ValueError("smooth_only and local_only are mutually exclusive")
+
+    total = 0.0
+    nodes, weights = gauss_legendre_01(order)
+    for fan in fans:
+        for r, wr in zip(nodes, weights, strict=True):
+            jacobian = fan.signed_jacobian(r)
+            for t, wt in zip(nodes, weights, strict=True):
+                source = fan.point(r, t)
+                if smooth_only:
+                    if sigma is None:
+                        raise ValueError("sigma is required for smooth split")
+                    kernel = ewald_log_smooth_kernel(target, source, sigma=sigma)
+                elif local_only:
+                    if sigma is None:
+                        raise ValueError("sigma is required for local split")
+                    kernel = ewald_log_local_kernel(target, source, sigma=sigma)
+                else:
+                    if (
+                        hypot(target[0] - source[0], target[1] - source[1])
+                        <= _COINCIDENT_TOL
+                    ):
+                        raise ValueError(
+                            "target coincides with a source quadrature node"
+                        )
+                    kernel = laplace_log_kernel(target, source)
+                total += kernel * _physical_density(source) * jacobian * wr * wt
+    return total
+
+
+def run_dmk_split_nearfield_experiment(
+    *,
+    orders: tuple[int, ...] = (4, 6, 8, 10, 14, 18),
+    sigmas: tuple[float, ...] = (0.08, 0.16, 0.32),
+    reference_order: int = 96,
+) -> DmkSplitExperimentReport:
+    """Run the DMK-style smooth/local split experiment on straight fan panels."""
+
+    fixtures: tuple[
+        tuple[str, tuple[Point2D, ...], tuple[tuple[str, Point2D], ...]], ...
+    ] = (
+        (
+            "convex_quad",
+            ((0.0, 0.0), (1.0, 0.0), (0.95, 0.42), (0.18, 0.86)),
+            (
+                ("interior", (0.55, 0.32)),
+                ("boundary_near_inside", (0.54, 0.025)),
+                ("boundary_near_outside", (0.54, -0.025)),
+                ("vertex_near_inside", (0.05, 0.045)),
+            ),
+        ),
+        (
+            "skew_quad",
+            ((0.0, 0.0), (1.15, 0.02), (0.72, 0.78), (0.08, 0.62)),
+            (
+                ("interior", (0.50, 0.34)),
+                ("boundary_near_inside", (0.56, 0.045)),
+                ("boundary_near_outside", (0.56, -0.025)),
+                ("vertex_near_inside", (0.06, 0.05)),
+            ),
+        ),
+    )
+
+    seed_quality: list[SeedQualitySample] = []
+    samples: list[DmkSplitSample] = []
+    for fixture_name, nodes, targets in fixtures:
+        panel = TrimmedPanel2D(outer=PanelLoop2D(nodes))
+        seed_options = (
+            ("interior_anchor", select_interior_anchor(panel)),
+            ("node_barycenter", _node_barycenter(nodes)),
+        )
+        for seed_mode, seed in seed_options:
+            seed_quality.append(
+                _panel_seed_quality(
+                    fixture=fixture_name,
+                    seed_mode=seed_mode,
+                    nodes=nodes,
+                    seed=seed,
+                )
+            )
+            fans = _panel_fans(nodes, seed=seed)
+            full_references = {
+                target_label: _panel_potential(
+                    fans,
+                    target,
+                    order=reference_order,
+                    sigma=None,
+                    smooth_only=False,
+                )
+                for target_label, target in targets
+            }
+            for sigma in sigmas:
+                smooth_references = {
+                    target_label: _panel_potential(
+                        fans,
+                        target,
+                        order=reference_order,
+                        sigma=sigma,
+                        smooth_only=True,
+                    )
+                    for target_label, target in targets
+                }
+                local_references = {
+                    target_label: _panel_potential(
+                        fans,
+                        target,
+                        order=reference_order,
+                        sigma=sigma,
+                        smooth_only=False,
+                        local_only=True,
+                    )
+                    for target_label, target in targets
+                }
+                for target_label, target in targets:
+                    full_reference = full_references[target_label]
+                    smooth_reference = smooth_references[target_label]
+                    local_reference = local_references[target_label]
+                    ratio_denominator = max(abs(full_reference), 1.0e-300)
+                    for order in orders:
+                        full_value = _panel_potential(
+                            fans,
+                            target,
+                            order=order,
+                            sigma=None,
+                            smooth_only=False,
+                        )
+                        smooth_value = _panel_potential(
+                            fans,
+                            target,
+                            order=order,
+                            sigma=sigma,
+                            smooth_only=True,
+                        )
+                        local_value = _panel_potential(
+                            fans,
+                            target,
+                            order=order,
+                            sigma=sigma,
+                            smooth_only=False,
+                            local_only=True,
+                        )
+                        analytic_local_value = _analytic_full_space_local_residual(
+                            target,
+                            sigma=sigma,
+                        )
+                        reconstructed_value = smooth_value + local_value
+                        analytic_reconstructed_value = (
+                            smooth_value + analytic_local_value
+                        )
+                        samples.append(
+                            DmkSplitSample(
+                                fixture=fixture_name,
+                                seed_mode=seed_mode,
+                                seed=seed,
+                                target_label=target_label,
+                                target=target,
+                                sigma=sigma,
+                                order=order,
+                                full_value=full_value,
+                                smooth_value=smooth_value,
+                                local_value=local_value,
+                                analytic_local_value=analytic_local_value,
+                                reconstructed_value=reconstructed_value,
+                                analytic_reconstructed_value=analytic_reconstructed_value,
+                                full_reference=full_reference,
+                                smooth_reference=smooth_reference,
+                                local_reference=local_reference,
+                                full_abs_error=abs(full_value - full_reference),
+                                smooth_abs_error=abs(smooth_value - smooth_reference),
+                                local_abs_error=abs(local_value - local_reference),
+                                analytic_local_abs_error=abs(
+                                    analytic_local_value - local_reference
+                                ),
+                                reconstructed_abs_error=abs(
+                                    reconstructed_value - full_reference
+                                ),
+                                analytic_reconstructed_abs_error=abs(
+                                    analytic_reconstructed_value - full_reference
+                                ),
+                                local_to_full_ratio=abs(local_reference)
+                                / ratio_denominator,
+                            )
+                        )
+
+    return DmkSplitExperimentReport(
+        reference_order=reference_order,
+        orders=orders,
+        sigmas=sigmas,
+        seed_quality=tuple(seed_quality),
+        samples=tuple(samples),
+    )
+
+
+def expected_scaled_laplace_point_potential(
+    base_potential: float,
+    source_mass: float,
+    scale_factor: float,
+) -> float:
+    """Return the exact 2D log-kernel scale law for a scaled source and target."""
+
+    if scale_factor <= 0.0:
+        raise ValueError("scale_factor must be positive")
+    return scale_factor**2 * (
+        base_potential - log(scale_factor) * source_mass / (2.0 * pi)
+    )
+
+
+def diagonal_remainder_sample(
+    fan: FanTemplateMap2D,
+    *,
+    r: float,
+    t: float,
+    delta: float,
+    direction: tuple[float, float] = (1.0, 0.5),
+) -> DiagonalRemainderSample:
+    """Compare the full mapped kernel with the local metric singular model."""
+
+    if delta <= 0.0:
+        raise ValueError("delta must be positive")
+    delta_r = delta * direction[0]
+    delta_t = delta * direction[1]
+    source = fan.point(r, t)
+    target = fan.point(r + delta_r, t + delta_t)
+    physical_distance = hypot(target[0] - source[0], target[1] - source[1])
+    model_distance = metric_model_distance(
+        fan,
+        r=r,
+        t=t,
+        delta_r=delta_r,
+        delta_t=delta_t,
+    )
+    laplace_kernel = -log(physical_distance) / (2.0 * pi)
+    metric_model_kernel = -log(model_distance) / (2.0 * pi)
+    return DiagonalRemainderSample(
+        delta=delta,
+        physical_distance=physical_distance,
+        model_distance=model_distance,
+        laplace_kernel=laplace_kernel,
+        metric_model_kernel=metric_model_kernel,
+        remainder=laplace_kernel - metric_model_kernel,
+    )
+
+
+def run_nearfield_template_experiment(
+    *,
+    order: int = 12,
+    scale_factor: float = 1.75,
+) -> NearfieldTemplateExperiment:
+    """Run the baseline analytic fan near-field template experiment."""
+
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    near_point_target = (0.47, 0.42)
+
+    def point_density(r: float, t: float) -> float:
+        return 1.0 + r - 0.5 * t
+
+    point_reference = point_target_laplace_potential(
+        fan,
+        near_point_target,
+        order=max(order + 8, 16),
+        source_density=point_density,
+    )
+    point_low_order = point_target_laplace_potential(
+        fan,
+        near_point_target,
+        order=max(order // 2, 2),
+        source_density=point_density,
+    )
+    scaled_target = _scale_point(near_point_target, scale_factor)
+    scaled_value = point_target_laplace_potential(
+        fan.scaled(scale_factor),
+        scaled_target,
+        order=max(order + 8, 16),
+        source_density=point_density,
+    )
+    density_mass = template_density_mass(
+        fan,
+        order=max(order + 8, 16),
+        density=point_density,
+    )
+    expected_scaled = expected_scaled_laplace_point_potential(
+        point_reference,
+        density_mass,
+        scale_factor,
+    )
+    remainders = tuple(
+        diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=delta)
+        for delta in (1.0e-1, 1.0e-2, 1.0e-3)
+    )
+    return NearfieldTemplateExperiment(
+        fan=fan,
+        order=order,
+        near_point_target=near_point_target,
+        point_target_reference=point_reference,
+        point_target_low_order=point_low_order,
+        point_target_abs_error=abs(point_low_order - point_reference),
+        scaled_near_point_target=scaled_target,
+        scaled_point_target_potential=scaled_value,
+        expected_scaled_point_target_potential=expected_scaled,
+        scaled_abs_error=abs(scaled_value - expected_scaled),
+        signed_area=fan.signed_area,
+        density_mass=density_mass,
+        scale_factor=scale_factor,
+        diagonal_remainders=remainders,
+    )
diff --git a/tests/evals/test_nearfield_templates.py b/tests/evals/test_nearfield_templates.py
new file mode 100644
index 0000000..e84eeea
--- /dev/null
+++ b/tests/evals/test_nearfield_templates.py
@@ -0,0 +1,482 @@
+from __future__ import annotations
+
+from math import cos, exp, log, pi, sin
+
+import pytest
+
+from cutkit.evals import (
+    FanTemplateMap2D,
+    diagonal_remainder_sample,
+    ewald_log_local_kernel,
+    ewald_log_local_radial_moment_2d,
+    ewald_log_local_sector_monomial_moment_2d,
+    ewald_log_smooth_kernel,
+    expected_scaled_laplace_point_potential,
+    point_target_laplace_potential,
+    point_target_smooth_ewald_potential,
+    run_curved_boundary_local_model_experiment,
+    run_dmk_split_nearfield_experiment,
+    run_graph_boundary_model_experiment,
+    run_nearfield_template_experiment,
+    run_precomputed_boundary_table_experiment,
+    run_taylor_boundary_model_experiment,
+    template_density_mass,
+)
+from cutkit.quadrature import gauss_legendre_01
+
+
+def _polar_sector_local_moment_reference(
+    *,
+    sigma: float,
+    x_power: int,
+    y_power: int,
+    theta_start: float,
+    theta_end: float,
+    radial_order: int = 220,
+    angular_order: int = 96,
+) -> float:
+    log_radius_min = log(sigma) - 32.0
+    log_radius_max = log(8.0 * sigma)
+    radial_nodes, radial_weights = gauss_legendre_01(radial_order)
+    angular_nodes, angular_weights = gauss_legendre_01(angular_order)
+    log_radius_width = log_radius_max - log_radius_min
+    theta_width = theta_end - theta_start
+    total = 0.0
+    for log_radius_node, log_radius_weight_base in zip(
+        radial_nodes, radial_weights, strict=True
+    ):
+        log_radius = log_radius_min + log_radius_width * log_radius_node
+        radius = exp(log_radius)
+        # The polar area factor contributes one radius and dr = radius d(log r).
+        radial_weight = log_radius_width * log_radius_weight_base * radius * radius
+        for theta_node, theta_weight_base in zip(
+            angular_nodes, angular_weights, strict=True
+        ):
+            theta = theta_start + theta_width * theta_node
+            theta_weight = theta_width * theta_weight_base
+            source = (radius * cos(theta), radius * sin(theta))
+            total += (
+                ewald_log_local_kernel((0.0, 0.0), source, sigma=sigma)
+                * source[0] ** x_power
+                * source[1] ** y_power
+                * radial_weight
+                * theta_weight
+            )
+    return total
+
+
+def test_fan_map_area_and_metric_are_consistent() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.25, 0.5),
+    )
+
+    assert fan.signed_area == pytest.approx(0.25)
+    metric = fan.local_metric(0.5, 0.25)
+    assert metric[0][0] > 0.0
+    assert metric[1][1] > 0.0
+    assert metric[0][1] == pytest.approx(metric[1][0])
+
+
+def test_point_target_potential_obeys_log_kernel_scale_law() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    scale_factor = 2.25
+    target = (0.8, 0.7)
+    value = point_target_laplace_potential(fan, target, order=6)
+    scaled = point_target_laplace_potential(
+        fan.scaled(scale_factor),
+        (scale_factor * target[0], scale_factor * target[1]),
+        order=6,
+    )
+    expected = expected_scaled_laplace_point_potential(
+        value,
+        template_density_mass(fan, order=6),
+        scale_factor,
+    )
+
+    assert scaled == pytest.approx(expected, abs=1.0e-13)
+
+
+def test_scaled_point_target_potential_uses_density_weighted_mass() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+
+    def source_density(r: float, t: float) -> float:
+        return 1.0 + r
+
+    scale_factor = 1.4
+    target = (0.8, 0.7)
+    value = point_target_laplace_potential(
+        fan,
+        target,
+        order=4,
+        source_density=source_density,
+    )
+    scaled = point_target_laplace_potential(
+        fan.scaled(scale_factor),
+        (scale_factor * target[0], scale_factor * target[1]),
+        order=4,
+        source_density=source_density,
+    )
+    expected = expected_scaled_laplace_point_potential(
+        value,
+        template_density_mass(fan, order=4, density=source_density),
+        scale_factor,
+    )
+
+    assert scaled == pytest.approx(expected, abs=1.0e-13)
+
+
+def test_metric_singular_remainder_shrinks_toward_diagonal() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    coarse = diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=1.0e-2)
+    fine = diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=1.0e-4)
+
+    assert abs(fine.remainder) < abs(coarse.remainder)
+    assert fine.physical_distance == pytest.approx(fine.model_distance, rel=1.0e-3)
+
+
+def test_point_target_path_converges_for_near_disjoint_target() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    target = (0.8, 0.7)
+
+    def density(r: float, t: float) -> float:
+        return 1.0 + r * t
+
+    coarse = point_target_laplace_potential(
+        fan,
+        target,
+        order=4,
+        source_density=density,
+    )
+    fine = point_target_laplace_potential(
+        fan,
+        target,
+        order=10,
+        source_density=density,
+    )
+    reference = point_target_laplace_potential(
+        fan,
+        target,
+        order=18,
+        source_density=density,
+    )
+
+    assert abs(fine - reference) < abs(coarse - reference)
+
+
+def test_point_target_path_rejects_coincident_source_node() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    nodes, _weights = gauss_legendre_01(3)
+    target = fan.point(nodes[1], nodes[1])
+
+    with pytest.raises(ValueError, match="coincides with a source quadrature node"):
+        point_target_laplace_potential(fan, target, order=3)
+
+
+def test_ewald_log_split_reconstructs_log_kernel() -> None:
+    target = (0.3, 0.4)
+    source = (0.9, -0.2)
+    sigma = 0.25
+
+    split = ewald_log_smooth_kernel(
+        target, source, sigma=sigma
+    ) + ewald_log_local_kernel(
+        target,
+        source,
+        sigma=sigma,
+    )
+
+    from cutkit.evals import laplace_log_kernel
+
+    assert split == pytest.approx(laplace_log_kernel(target, source), abs=1.0e-13)
+
+
+def test_ewald_local_full_plane_moments_match_known_values() -> None:
+    sigma = 0.17
+
+    mass = ewald_log_local_sector_monomial_moment_2d(
+        sigma=sigma,
+        x_power=0,
+        y_power=0,
+        theta_start=0.0,
+        theta_end=2.0 * pi,
+    )
+    xx_moment = ewald_log_local_sector_monomial_moment_2d(
+        sigma=sigma,
+        x_power=2,
+        y_power=0,
+        theta_start=0.0,
+        theta_end=2.0 * pi,
+    )
+    xy_moment = ewald_log_local_sector_monomial_moment_2d(
+        sigma=sigma,
+        x_power=1,
+        y_power=1,
+        theta_start=0.0,
+        theta_end=2.0 * pi,
+    )
+
+    assert mass == pytest.approx(sigma * sigma / 4.0, rel=1.0e-13)
+    assert xx_moment == pytest.approx(sigma**4 / 16.0, rel=1.0e-13)
+    assert xy_moment == pytest.approx(0.0, abs=1.0e-19)
+
+
+def test_ewald_local_half_plane_moments_match_polar_reference() -> None:
+    sigma = 0.13
+    theta_start = 0.0
+    theta_end = pi
+
+    for x_power, y_power in ((0, 0), (1, 0), (0, 1), (2, 0), (0, 2)):
+        analytic = ewald_log_local_sector_monomial_moment_2d(
+            sigma=sigma,
+            x_power=x_power,
+            y_power=y_power,
+            theta_start=theta_start,
+            theta_end=theta_end,
+        )
+        reference = _polar_sector_local_moment_reference(
+            sigma=sigma,
+            x_power=x_power,
+            y_power=y_power,
+            theta_start=theta_start,
+            theta_end=theta_end,
+        )
+        assert analytic == pytest.approx(reference, rel=5.0e-4, abs=1.0e-12)
+
+    mass = ewald_log_local_sector_monomial_moment_2d(
+        sigma=sigma,
+        x_power=0,
+        y_power=0,
+        theta_start=theta_start,
+        theta_end=theta_end,
+    )
+    assert mass == pytest.approx(sigma * sigma / 8.0, rel=1.0e-13)
+
+
+def test_ewald_local_wedge_moments_match_polar_reference() -> None:
+    sigma = 0.11
+    theta_start = -0.2
+    theta_end = 0.85 * pi
+
+    for x_power, y_power in ((0, 0), (1, 0), (0, 1), (2, 0), (1, 1)):
+        analytic = ewald_log_local_sector_monomial_moment_2d(
+            sigma=sigma,
+            x_power=x_power,
+            y_power=y_power,
+            theta_start=theta_start,
+            theta_end=theta_end,
+        )
+        reference = _polar_sector_local_moment_reference(
+            sigma=sigma,
+            x_power=x_power,
+            y_power=y_power,
+            theta_start=theta_start,
+            theta_end=theta_end,
+        )
+        assert analytic == pytest.approx(reference, rel=5.0e-4, abs=1.0e-12)
+
+    mass = ewald_log_local_sector_monomial_moment_2d(
+        sigma=sigma,
+        x_power=0,
+        y_power=0,
+        theta_start=theta_start,
+        theta_end=theta_end,
+    )
+    opening_angle = theta_end - theta_start
+    assert mass == pytest.approx(opening_angle * sigma * sigma / (8.0 * pi))
+
+
+def test_ewald_local_radial_moment_rejects_invalid_inputs() -> None:
+    with pytest.raises(ValueError, match="sigma must be positive"):
+        ewald_log_local_radial_moment_2d(sigma=0.0, monomial_degree=0)
+    with pytest.raises(ValueError, match="monomial_degree must be nonnegative"):
+        ewald_log_local_radial_moment_2d(sigma=0.1, monomial_degree=-1)
+
+
+def test_curved_boundary_local_model_beats_full_plane_on_cad_style_cell() -> None:
+    report = run_curved_boundary_local_model_experiment(
+        sigmas=(0.02, 0.04),
+        orders=(16, 32),
+        reference_order=80,
+    )
+
+    finest_by_sigma = {
+        sample.sigma: sample for sample in report.samples if sample.order == 32
+    }
+
+    assert len(finest_by_sigma) == 2
+    for sample in finest_by_sigma.values():
+        assert sample.reference_value > 0.0
+        assert sample.half_plane_abs_error < sample.full_plane_abs_error
+        assert sample.half_plane_rel_error < 0.03
+        assert sample.full_plane_rel_error > 1.0
+
+
+def test_curved_boundary_exact_model_reaches_machine_precision() -> None:
+    coarse = run_curved_boundary_local_model_experiment(
+        sigmas=(0.02, 0.04),
+        orders=(16,),
+        reference_order=24,
+    )
+    fine = run_curved_boundary_local_model_experiment(
+        sigmas=(0.02, 0.04),
+        orders=(16,),
+        reference_order=32,
+    )
+
+    coarse_by_sigma = {sample.sigma: sample for sample in coarse.samples}
+    fine_by_sigma = {sample.sigma: sample for sample in fine.samples}
+
+    for sigma, fine_sample in fine_by_sigma.items():
+        coarse_value = coarse_by_sigma[sigma].curved_boundary_value
+        assert fine_sample.curved_boundary_value == pytest.approx(
+            coarse_value,
+            rel=1.0e-10,
+            abs=1.0e-14,
+        )
+
+
+def test_high_order_taylor_boundary_model_reaches_machine_precision() -> None:
+    report = run_taylor_boundary_model_experiment(
+        sigmas=(0.02,),
+        boundary_orders=(2, 4, 6, 8),
+        angular_order=32,
+    )
+    errors = {sample.boundary_order: sample.rel_error for sample in report.samples}
+
+    assert errors[4] < errors[2]
+    assert errors[6] < errors[4]
+    assert errors[8] < 1.0e-10
+
+
+def test_high_order_taylor_boundary_model_handles_smaller_windows() -> None:
+    report = run_taylor_boundary_model_experiment(
+        sigmas=(0.005, 0.01),
+        boundary_orders=(4, 6),
+        angular_order=32,
+    )
+
+    finest_errors = {
+        sample.sigma: sample.rel_error
+        for sample in report.samples
+        if sample.boundary_order == 6
+    }
+
+    assert finest_errors[0.005] < 1.0e-10
+    assert finest_errors[0.01] < 1.0e-10
+
+
+def test_precomputed_boundary_table_interpolates_to_machine_precision() -> None:
+    report = run_precomputed_boundary_table_experiment(
+        table_sigmas=(0.004, 0.008, 0.012, 0.016, 0.020, 0.024, 0.028),
+        eval_sigmas=(0.010, 0.018, 0.026),
+        boundary_order=8,
+        angular_order=32,
+    )
+
+    assert max(sample.rel_error for sample in report.samples) < 1.0e-10
+
+
+def test_direct_graph_strip_quadrature_is_not_enough_for_machine_precision() -> None:
+    report = run_graph_boundary_model_experiment(
+        sigmas=(0.02,),
+        boundary_orders=(8,),
+        strip_orders=(64, 128),
+        angular_order=32,
+    )
+
+    errors = {sample.strip_order: sample.rel_error for sample in report.samples}
+
+    assert errors[128] < errors[64]
+    assert errors[128] > 1.0e-10
+
+
+def test_curved_boundary_folded_reference_converges_with_order() -> None:
+    report = run_curved_boundary_local_model_experiment(
+        sigmas=(0.04,),
+        orders=(10, 16, 32),
+        reference_order=80,
+    )
+
+    errors = {sample.order: sample.folded_abs_error for sample in report.samples}
+
+    assert errors[32] < errors[16]
+    assert errors[16] < errors[10]
+
+
+def test_smoothed_ewald_potential_allows_coincident_quadrature_node() -> None:
+    fan = FanTemplateMap2D(
+        vertex=(0.0, 0.0),
+        edge_start=(1.0, 0.0),
+        edge_end=(0.35, 0.9),
+    )
+    nodes, _weights = gauss_legendre_01(3)
+    target = fan.point(nodes[1], nodes[1])
+
+    value = point_target_smooth_ewald_potential(
+        fan,
+        target,
+        order=3,
+        sigma=0.2,
+    )
+
+    assert value == pytest.approx(value)
+
+
+def test_dmk_split_experiment_records_local_residual_reconstruction() -> None:
+    report = run_dmk_split_nearfield_experiment(
+        orders=(4,),
+        sigmas=(0.16,),
+        reference_order=8,
+    )
+
+    sample = report.samples[0]
+    assert sample.reconstructed_value == pytest.approx(
+        sample.smooth_value + sample.local_value
+    )
+    assert sample.analytic_reconstructed_value == pytest.approx(
+        sample.smooth_value + sample.analytic_local_value
+    )
+    assert sample.local_reference == pytest.approx(
+        sample.full_reference - sample.smooth_reference,
+        abs=1.0e-12,
+    )
+    assert sample.reconstructed_abs_error >= 0.0
+    assert sample.analytic_local_abs_error >= 0.0
+    assert sample.analytic_reconstructed_abs_error >= 0.0
+
+
+def test_baseline_experiment_records_feasibility_checks() -> None:
+    result = run_nearfield_template_experiment(order=5, scale_factor=1.5)
+
+    assert result.scaled_abs_error < 1.0e-13
+    assert result.point_target_abs_error > 0.0
+    assert result.scaled_point_target_potential == pytest.approx(
+        result.expected_scaled_point_target_potential,
+        abs=1.0e-13,
+    )
+    assert result.diagonal_remainders[-1].delta == pytest.approx(1.0e-3)
+    assert abs(result.diagonal_remainders[-1].remainder) < abs(
+        result.diagonal_remainders[0].remainder
+    )