# A physical scalar effective-source calculation

Use geometrized units `G=c=1` in this benchmark. It solves the time-dependent
scalar wave equation on Schwarzschild spacetime for a unit charge on a fixed
circular geodesic at `r0/M = 6, 10`.
It uses the established public **SelfForce-1D** discontinuous Galerkin solver,
not a manufactured field. The reported quantity is the Schwarzschild
**covariant** radial self-force `Fr = q ∂r ΦR`, with `M=q=1`; force values denote
the dimensionless combination `M² Fr/q²` in these geometrized units. It is a
coordinate component, not an invariant or a gravitational self-force.
The r0=6M geodesic is marginally stable; holding its worldline fixed defines
the field problem and does not simulate an inspiral through that orbit.

The archived results have deliberately modest accuracy. The purpose is to
separate numerical convergence, startup relaxation, angular truncation and
source regularity in a reproducible physical calculation.

## Equations and normalization

Choose the equatorial orbit `θ=π/2, φ=Ωt`. The scalar equation is
`Box Φ = -4πρ`, with spherical harmonics normalized by
`∫ |Ylm|² dΩ = 1`. For `Φ = Σ ψlm(t,r) Ylm(θ,φ)/r`, the radial wave operator is

```
L = -∂t² + ∂r*² - Vℓ,
Vℓ = f [ℓ(ℓ+1)/r² + 2M/r³],  f = 1-2M/r,
dr*/dr = 1/f,  Ω = sqrt(M/r0³),  ut = (1-3M/r0)^(-1/2).
```

Within the worldtube the evolved field is `ψres = ψret - ψP` and
`L ψres = -L ψP` away from the particle. Outside it, the evolved field is
retarded. Analytic puncture jumps are imposed through the numerical fluxes at
the tube boundaries. The upstream source provider returns `L ψP`; its positive
sign on the `∂t² ψres = ... + L ψP` evolution equation is therefore correct.

The observer reconstructs the physical radial derivative as
`Σlm Re[Ylm (ψr*/(r-2M) - ψ/r²)]`. It counts `m=0` once and conjugate positive
and negative `m` pairs twice. Independent probes of the puncture verify that
its averaged worldline derivative is `Br` (independent of ℓ at fixed radius)
and its two-sided jump is `2Ar,+ L`, where `L=ℓ+1/2`.
Thus the raw residual multipoles in this package equal the averaged retarded
force multipoles minus `Br`.

## What is varied

The radial method is nodal discontinuous Galerkin with `order+1` nodes per
element and a five-stage, fourth-order low-storage Runge–Kutta integrator.
The particle lies at a shared element boundary. The diagnostic single-mode
runs vary the polynomial degree at fixed physical time step and vary the time
step at fixed degree. Full multipole runs also vary worldtube width and the
hyperboloidal transition position.

Hyperboloidal compactification includes the actual future horizon and future
null infinity. Its endpoint characteristics implement the retarded boundary
conditions; there is no artificial finite-radius outgoing boundary. Changing
`t_size` tests the coordinate transition, not the position of infinity.
`world_tube_width` counts elements on **each** side of the particle. The sharp
worldtube and its flux jumps are varied; no family of smooth spatial windows
is tested.

The upstream **second-order puncture is fixed**. Its harmonic expansion retains
radial terms through `Δr⁴` and rotated `m'=0,±1,±2`. Individual harmonic effective
sources are continuous with a derivative kink; the residual harmonics are C².
The probes approach the particle from both sides and verify this regularity.
This statement does not assign the same regularity to the unsummed 3D source.
No puncture-order refinement was performed.

Some final runs initialize the static monopole analytically through the
upstream data-reading interface. For unit charge,
`(r²f Φ00')' = -4πY00 δ(r-r0)/ut`: the inner field is constant and the outer
field is `-2πY00 log(f)/(M ut)`. The two particle grid rows keep distinct
one-sided derivatives. Higher modes start from zero and turn on smoothly.
In this code `turn_on_source_smoothly=False` exempts only modes through
`exact_initial_data_lmax=0`: see
`Src/Equations/PDE/ScalarSchwarzschild/submodule_scalar_schw_implementation.f90`,
lines 367–375, and `Src/EffSource/ScalarSchwarzschild/scalar_schw_eff_interface.cc`,
lines 81–91, at the pinned revision. The remaining modes retain the time window
and both of its derivatives.
This is relaxation toward the eternal circular-orbit solution, not a claim
about the physical transient of a uniformly switched-on charge.
The final r10 run and its tube/mapping variants all use this same initialization;
the r6 runs use zero initial data. Comparisons of a seeded full run with a
zero-start run combine relaxation and discretization effects. The separately
matched single-mode cases isolate spatial and temporal convergence.

## Angular sum and independent comparison

`analyze_runs.py` reads the raw residual modes and deliberately excludes the
solver's optional fitted-tail columns. It subtracts the analytic `F2`, `F4`
and `F6` terms and reports partial sums through `ℓ=8,12,16`. The terms are read
from `../scalar-mode-sum/scalar-parameters-60digits.json`; each includes its
complete angular denominator. Their infinite sums vanish. This postprocessing
improves angular convergence without changing the PDE puncture order.

An independent frequency-domain wave calculation supplies the same-cutoff
comparison in `../scalar-mode-sum/checks.json`. Both methods share the analytic
regularization parameters, so their agreement does not independently validate
those coefficients. The puncture A/B probe, independent wave-operator probe,
and analytic zero-sum identities supply separate checks. Refinement and
late-time spreads are empirical error diagnostics, not rigorous bounds.

Warburton and Wardell (2014), Table I, uses **thirty modes, ℓ=0…29**, with
higher-order subtraction. Its numbers are finite-cutoff method comparisons,
not infinite-sum reference values. At exactly that cutoff the independent
mode-sum gives `1.677283080352924e-4` at r6 and `1.3784482233828664e-5` at r10,
reproducing the published column on rounding. Our ℓ≤16 result is first compared
with the independent ℓ≤16 result; its remaining angular tail is kept separate.

## Selected results and accuracy

The final degree-14 runs give the following F2/F4/F6-subtracted sums. Extra
digits identify numerical differences; they are not a claim of that many
accurate digits in the continuum force.

| r0/M | PDE sum through ℓ=16 | Independent sum at the same cutoff | Relative difference |
|---:|---:|---:|---:|
| 6 | 1.6773351188452784e-4 | 1.677262880429488e-4 | 4.31e-5 |
| 10 | 1.3785077695498861e-5 | 1.3784462518545389e-5 | 4.46e-5 |

The final 50M time spans change these sums by 1.16e-9 and 2.60e-10. At r10,
doubling the worldtube half-width changes the fine result by 1.38e-9,
versus 2.59e-8 at the coarser time step. The hyperboloidal-transition change
is 1.58e-13. Those localization/mapping variations were performed at r10.
The independent infinite-sum results imply angular remainders above ℓ=16
of 2.05e-9 and 2.01e-11, respectively. No such inferred tail is inserted into
the time-domain output. Direct differences from the infinite-sum references
are 5.17e-9 and 5.95e-10; partial cancellation with the missing angular tail
makes these smaller than the same-cutoff differences.

The calculation reproduces the force at roughly the 1e-4 relative level.
Comparison errors, settling changes and angular remainders are distinct;
there is no rigorous total continuum-error bound. The fixed double-precision
PDE solver was not rerun at higher arithmetic precision. `results.json`
retains the separate spatial, temporal, settling, tube and mapping checks.

## Reproduce

Keep the downloaded helpers together and the two shared data files in a sibling
directory; the scripts resolve these paths relative to their own location:

```text
computations/
  scalar-effective-source/
    build_solver.py
    run_suite.py
    analyze_runs.py
    probe_source.cpp
    cases.json
  scalar-mode-sum/
    scalar-parameters-60digits.json
    checks.json
```

Requirements: Python with NumPy, Git, GNU Fortran and C++ with OpenMP, GSL,
and BLAS/LAPACK. The tested environment is in `environment.json`. The macOS
build uses the existing Accelerate framework; the Linux BLAS/LAPACK branch
of the helper has not been executed here. Nothing is installed globally.

If GSL is unavailable, its official source is
`https://ftp.gnu.org/gnu/gsl/gsl-2.8.tar.gz` (SHA-256
`6a99eeed15632c6354895b1dd542ed5a855c0f15d9ad1326c6fe2b2c9e423190`).
Build it in a scratch directory with `./configure --prefix="$GSL_PREFIX"
--disable-shared`, `make`, and `make install`, where `GSL_PREFIX` is an absolute
path in that scratch directory. The supplied helper then fetches the exact
solver revision and builds it serially, avoiding upstream first-build Fortran
module ordering races. It changes only `SConstruct` and `Src/SConscript` and
checks that all numerical source files remain unchanged.

From this package directory, with scratch paths chosen by the reader:

```sh
python3 build_solver.py --work "$WORK" --gsl-prefix "$GSL_PREFIX" --cxx g++-15 --probe
python3 run_suite.py fine_r10 --case --work "$WORK" --threads 2
python3 analyze_runs.py --runs "$WORK/runs" \
  --parameters ../scalar-mode-sum/scalar-parameters-60digits.json \
  --reference-r10 ../scalar-mode-sum/checks.json \
  --output "$WORK/results.json" fine_r10
```

`cases.json` gives all selected runs and their exact parameter overrides.
Run any other named case with `run_suite.py NAME --case --work "$WORK"`.
Each run retains its namelist, raw time series, log and completion record in
the scratch directory. Existing runs are not overwritten. Full cases are
more expensive than the single-ℓ diagnostics; run only the desired cases.
Selected raw multipoles, times, parameters and comparisons are retained in
`results.json`, rather than every temporary time-series output.

To regenerate the complete selected matrix after the build, use
`python3 run_suite.py --matrix --work "$WORK" --threads 2
--output-directory "$OUTPUT"`. Add `--resume` to reuse completed runs after
the helper checks their parameters. The default analytic data and reference
paths resolve to the sibling mode-sum package; `--parameters` and `--reference`
can override them. The output directory receives the selected results, source
checks, probe samples and environment metadata. Raw histories stay in the
scratch directory. Timings depend on concurrent machine load and are not a
performance benchmark.

With `--probe`, the build helper compiles `probe_source.cpp` against the
unchanged upstream source-provider and Wigner-matrix objects and writes
`$WORK/source-probes.csv`. Its selected output is retained as `source-probes.csv`.
Pass `--probes "$WORK/source-probes.csv" --probe-output "$WORK/source-checks.json"`
to the analysis command to recompute the A/B, operator and regularity checks.

## Sources

- SelfForce-1D, Peter Diener, Barry Wardell and Niels Warburton:
  <https://bitbucket.org/peterdiener/selfforce-1d>, revision
  `8390ca611ebcb2b644df6f59c15b1bafedf4c9f2` (ET_2026_05).
  The official [Einstein Toolkit download instructions](https://www.einsteintoolkit.org/download.html)
  identify this repository. No solver source or linked binaries are redistributed here.
- A. Heffernan, A. C. Ottewill, N. Warburton, B. Wardell and P. Diener,
  “Accelerated motion and the self-force in Schwarzschild spacetime,”
  *Classical and Quantum Gravity* **35**, 194001 (2018),
  [DOI:10.1088/1361-6382/aad420](https://doi.org/10.1088/1361-6382/aad420).
  [Open PDF, arXiv:1712.01098v3, § IV.B, p. 10](https://arxiv.org/pdf/1712.01098).
- L. M. Diaz-Rivera, E. Messaritaki, B. F. Whiting and S. Detweiler,
  “Scalar field self-force effects on orbits about a Schwarzschild black hole,”
  *Physical Review D* **70**, 124018 (2004),
  [DOI:10.1103/PhysRevD.70.124018](https://doi.org/10.1103/PhysRevD.70.124018).
  [Open PDF, Table I, p. 6](https://arxiv.org/pdf/gr-qc/0410011).
- N. Warburton and B. Wardell, “Applying the effective-source approach to
  frequency-domain self-force calculations,” *Physical Review D* **89**,
  044046 (2014), [DOI:10.1103/PhysRevD.89.044046](https://doi.org/10.1103/PhysRevD.89.044046).
  [Open PDF, Table I and adjacent text, p. 10](https://arxiv.org/pdf/1311.3104).

See `NOTICES.txt` for licensing and the sibling mode-sum package for the
provenance and regeneration of the shared analytic regularization parameters.