# Scalar Schwarzschild mode-sum calculation For this numerical example use `G=c=1`. This calculation evaluates the self-force of a minimally coupled, massless scalar charge on an exactly circular Schwarzschild geodesic at areal radius `r/M=6` or `10`. It is a pedagogical scalar calculation, not a gravitational self-force calculation. The reported quantity is the coordinate component `M^2 F_r/q_s^2`, with `F_r=q_s partial_r Phi_R` in Schwarzschild coordinates. It is not a coordinate-invariant force magnitude. ## Download the files For fast retained-data reanalysis, save these five files in one folder: [analyze_modes.py](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/analyze_modes.py), [scalar-parameters-60digits.json](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/scalar-parameters-60digits.json), [modes-r6-wp55.json](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/modes-r6-wp55.json), [modes-r10-wp55.json](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/modes-r10-wp55.json), and [checks.json](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/checks.json). Only Python and `mpmath==1.3.0` are needed for that route. For full regeneration, also download [compute_modes.py](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/compute_modes.py), [compute_high_precision.wls](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/compute_high_precision.wls), and [export-scalar-parameters.wls](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/export-scalar-parameters.wls). Keep [NOTICES.txt](https://blackhole.org/computations/two-body-dynamics/scalar-mode-sum/NOTICES.txt) with the parameter data for its source attribution and license. The commands below specify the separately installed solver dependencies. ## Field and source normalization With `G=c=1`, set `M=q_s=1`, `Omega=r0^(-3/2)`, `u^t=(1-3/r0)^(-1/2)`, and `Delta=r(r-2)`. The conventions are ```text Box Phi = -4 pi rho Phi = sum_lm R_lm(r) Y_lm(theta,phi) exp(-i m Omega t) integral |Y_lm|^2 dOmega = 1 (Delta R')' + [omega^2 r^4/Delta - l(l+1)] R = -4 pi q_s conjugate(Y_lm(pi/2,0)) delta(r-r0)/u^t . ``` The official Black Hole Perturbation Toolkit packages supply homogeneous solutions `R_in` (ingoing at the future horizon) and `R_up` (outgoing at infinity). For the scalar Teukolsky variable, `R_in ~ exp(-i omega r*)` at the horizon and `R_up ~ exp(i omega r*)/r` at infinity. Define ```text W = Delta (R_in R_up' - R_up R_in') C_lm = -4 pi q_s^2 |Y_lm(pi/2,0)|^2/(u^t W) F_r,lm,+ = C_lm R_in(r0) R_up'(r0) F_r,lm,- = C_lm R_up(r0) R_in'(r0). ``` Terms with odd `l+m` vanish at the equator. The `m=0` term has weight one; the real part of every `m>0` term has weight two. No additional `1/r` conversion is applied to this scalar Teukolsky radial variable. The source jump is `F_r,l,+ - F_r,l,- = -(2l+1)/(u^t Delta0)`. The static monopole has `F_r,0,+ = -1/(u^t Delta0)` and `F_r,0,-=0`. ## Subtraction and the remaining tail `scalar-parameters-60digits.json` contains the official covariant scalar radial subtraction functions of orders `[-1,0,2,4,6]`. Each entry already contains its complete multipole dependence. The order `-1` entry is the outer-side `A_r L`, where `L=l+1/2`; the inner-side entry has the opposite sign. The order `0` entry is `B_r`. For these circular scalar orbits the usual `C_r` and finite remainder `D_r` vanish. Averaging the two sides cancels `A_r L`; it does not remove `B_r`. The higher-order entries are constants times ```text b_k(l) = 1/product_{j=1}^k [(2l+1)^2 - 4j^2], k=1,2,3. ``` Each `b_k` sums to zero from `l=0` to infinity, so subtracting these complete rational functions accelerates convergence without a finite add-back. They must not be replaced by their leading inverse powers without changing the add-back prescription. After order 6 subtraction the leading residual is order `L^-8`. The remaining tail is fitted to `b_4,b_5,...` over explicitly reported multipole windows. The exact sum above cutoff `N` is ```text T_k(N) = (N+1)/[(2k-1) product_{j=1}^k ((2N+2)^2-(2j-1)^2)] . ``` Indeed, `T_k(N-1)-T_k(N)=b_k(N)` and `T_k(infinity)=0`. The analysis script checks the tail against an independent direct finite sum, checks the official entries' rational basis, and varies the number of fit terms and the fit window. Small fit residuals alone are not an error bound. ## Run the calculation The tested environment is Python 3.12.14 on macOS ARM64 with `pybhpt 1.0.0`, `numpy 2.2.6`, `scipy 1.18.1`, `mpmath 1.3.0`, and Wolfram Mathematica 15.0.0. Dependency source revisions and the exact tested pybhpt wheel hash are recorded in `checks.json`. Dependencies are installed separately and are not distributed with this example. Run from this directory. Use a new virtual environment: ```sh python3.12 -m venv venv venv/bin/python -m pip install pybhpt==1.0.0 numpy==2.2.6 scipy==1.18.1 mpmath==1.3.0 venv/bin/python compute_modes.py --radius 6 --lmax 80 --flux-lmax 30 --output modes-r6-auto.json venv/bin/python compute_modes.py --radius 10 --lmax 80 --flux-lmax 30 --output modes-r10-auto.json ``` This fast calculation uses the official pybhpt `AUTO` method, checks Wronskian constancy at `r0-0.05,r0,r0+0.05`, and checks source jumps and the analytic static monopole. It also computes the covariant time component from `q_s partial_t Phi=-i omega q_s Phi`. The optional flux calculation uses the package's separate inhomogeneous point-source solver and verifies `flux_infinity+flux_horizon=F_t/u^t` at the same multipole cutoff. The scalar rest mass is constant along this helically symmetric circular orbit, and positive radiation carries away the particle's Killing energy. The time component also provides a published dissipative-force comparison. Double precision reveals the convergence pattern but eventually loses the small high-order residuals. For the higher-precision calculation, place exact checkouts of these official packages under `vendor/`: | Folder | Source revision | |---|---| | `Teukolsky` | `39ee4fb7266e780e0b6d4ae0c4a0556c91323264` | | `KerrGeodesics` | `c7c57a936236ddcce98963481440bca0c433b3c3` | | `SpinWeightedSpheroidalHarmonics` | `48cd35591f25ee5910445abad5f18d2b33c140ad` | | `RegularizationParameters` | `ad85c6f0088bc47f3a0a2ab22c8d1fd5b36b8681` | All four repositories are in the [Black Hole Perturbation Toolkit GitHub organization](https://github.com/BlackHolePerturbationToolkit). `export-scalar-parameters.wls` regenerates the distributed parameter data from the original `SchwarzschildRegularizationParameters.nb`, selecting only its first five scalar radial-force cells. It independently checks the circular `A_r` and `B_r` reductions before export. ```sh WolframKernel -noinit -noprompt -script export-scalar-parameters.wls WolframKernel -noinit -noprompt -script compute_high_precision.wls 6 40 35 modes-r6-wp35.json WolframKernel -noinit -noprompt -script compute_high_precision.wls 10 40 35 modes-r10-wp35.json WolframKernel -noinit -noprompt -script compute_high_precision.wls 6 40 55 modes-r6-wp55.json WolframKernel -noinit -noprompt -script compute_high_precision.wls 10 40 55 modes-r10-wp55.json venv/bin/python analyze_modes.py modes-r6-wp55.json --output analysis-r6.json venv/bin/python analyze_modes.py modes-r10-wp55.json --output analysis-r10.json venv/bin/python analyze_modes.py --compare-directory . --output checks-rerun.json ``` The precision argument is a base precision. Direct MST working precision is `base+l+20`. Increasing the base from 35 to 55 therefore raises every working precision by 20 digits. Every direct MST mode is checked against the independent pybhpt result. A disagreement larger than `10^-12` in an individual complex force triggers the official Mathematica numerical integration method with MST endpoint data; unresolved disagreement stops the run. This guard is necessary because direct MST evaluation at some large frequencies gave incorrect results despite increased nominal precision. The selected method and precision for every mode are retained. The comparison tolerance detects such failures; it does not establish accuracy below that threshold. The allowances below describe observed precision, cutoff and tail-fit changes. They do not bound possible stable or shared solver bias. The scripts save a checkpoint after each completed multipole. Repeating a higher-precision command with an incomplete output resumes it; incompatible inputs and already-complete outputs are rejected. Rename distributed raw files before generating an independent rerun, or choose new output names. For a quick reanalysis without either radial solver, install only `mpmath==1.3.0` and run the final `--compare-directory` command on the distributed files. It recomputes the regularization and tail fits from the retained raw 55-base-precision modes. When the original 35-base-precision and fast raw files are absent, it uses their retained one-sided multipoles and flux inputs in `checks.json`; no accepted force or fitted tail is used as input. This checks the published analysis, not a fresh solution of the radial boundary-value problem. The recorded solver environment is kept separate from the actual environment running this reanalysis. ## Reading the results `modes-r6-wp55.json` and `modes-r10-wp55.json` retain the individual complex homogeneous solutions and derivatives, Wronskians, and both one-sided force modes through `l=40`. The separate fast calculation checks cutoffs through `l=80`; `checks.json` distinguishes the two precision regimes. Mathematica decimal strings include precision marks; `analyze_modes.py` parses these without first converting to binary floats. `checks.json` records the comparison runs and empirical convergence allowances. These summarize the stated refinements; they are not total validated errors, rigorous interval enclosures or a certification of the external solver. Stable or shared method bias is not bounded. The selected tail fit uses four rational functions `b_4,...,b_7` over `l=25,...,40`. The tail-corrected point estimate and empirical convergence allowance are: | Radius | Tail-corrected `M^2 F_r/q_s^2` | Empirical convergence allowance | |---|---:|---:| | `6M` | `0.000167728340233` | `2e-15` | | `10M` | `0.000013784482575668` | `2e-18` | The allowances exceed the observed changes when the working precision is raised by 20 digits, the cutoff increases from 35 to 40, the fit window changes from 25–40 to 20–35 or 30–40, and the number of tail terms changes from four to three or five. These are independent variations: a fit with a very small residual is not accepted on that basis alone. `checks.json` retains each comparison, both one-sided sums, the low-precision multipoles, and the per-mode subtractions needed for independent analysis. The fast 80-mode calculation is a separate check; its cancellation error is too large to establish those smaller allowances as total error estimates. In `checks.json`, `radii["6"]` and `radii["10"]` each contain `accepted_force`, `estimated_abs_error`, and `multipoles`. The existing `estimated_abs_error` field stores the empirical convergence allowance, not a total validated error estimate. Every multipole records `plus`, `minus`, `average`, the complete functions `A,B,F2,F4,F6`, and the successive residuals `AB,AB2,AB24,AB246`. Values are decimal strings. The archived 55-base-precision raw files retain only `l<=40`; if selected from a longer checkpoint, their retention note identifies that fact. The benchmark concerns a fixed circular orbit. It does not evolve an inspiral, describe a finite-mass-ratio binary, or establish a gravitational gauge-invariant observable. The scalar equation has retarded horizon and infinity conditions throughout. ## Primary comparisons One comparison deliberately holds the cutoff fixed. Warburton and Wardell's 2014 Table I computes thirty multipoles, `l=0,...,29`, with higher-order regularization parameters. The present order 6-subtracted partial sums are | Radius | Present partial sum through `l=29` | Published mode-sum value | |---|---:|---:| | `6M` | `0.0001677283080352924` | `0.00016772830804` | | `10M` | `0.00001378448223382866` | `0.000013784482234` | Both reproduce the printed values at their stated precision. These are finite-cutoff sums, so their difference from a tail-corrected result is expected. Close agreement between two methods at the same cutoff is not an equally accurate estimate of the infinite sum. At `6M`, the published thirty-mode value is rounded at `1e-14` in the force; the independent full-force result cited below supplies seven significant figures. Neither validates the `2e-15` convergence allowance as a total error estimate. At `10M`, the independent Heffernan–Ottewill–Wardell full-force value agrees within `6.73e-20`, providing stronger external evidence there. - A. Heffernan, A. C. Ottewill, and B. Wardell, “High-order expansions of the Detweiler–Whiting singular field in Schwarzschild spacetime,” *Physical Review D* **86**, 104023 (2012), [doi:10.1103/PhysRevD.86.104023](https://doi.org/10.1103/PhysRevD.86.104023), § V F, p. 27, Fig. 2 and surrounding text. Their `r0=10M` result is `F_r=0.000013784482575667959(3)`; the parenthesis estimates their finite mode-tail error, not their entire numerical error. - N. Warburton and L. Barack, “Self force on a scalar charge in Kerr spacetime: Circular equatorial orbits,” *Physical Review D* **81**, 084039 (2010), [doi:10.1103/PhysRevD.81.084039](https://doi.org/10.1103/PhysRevD.81.084039). The Schwarzschild `a=0` row of Table III, p. 15 of arXiv1003.1860v3 gives `M^2 F_r/q_s^2=1.677283e-4` at `r0=6M`. Table II, p. 14 gives the time-component comparisons at both radii. Do not assign extra digits to the radial value. - 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). Table I and the following paragraph on printed p. 10 of [arXiv1311.3104v3, Open PDF](https://arxiv.org/pdf/1311.3104) specify the finite-cutoff comparison above. Section II gives the scalar source and Green-function construction. For formulas and dependency provenance, see the official [Teukolsky module](https://bhptoolkit.org/modules/teukolsky/), [pybhpt module](https://bhptoolkit.org/modules/pybhpt/), and [regularization-parameter module](https://bhptoolkit.org/modules/regularizationparameters/).