Recompute an EMRI frequency evolution from a published flux table

What you will reproduce

Three circular, equatorial inspirals with $q=10^{-5}$ and $M=10^6\,M_\odot$ : zero spin from $r=10$ to just outside the ISCO, spin 0.9 from $r=10$ to $r=6$ , and spin $-0.9$ from $r=20$ to $r=10$ . For each, the script reports the time and the number of wave cycles with the table flux and with the quadrupole flux, and a frequency track with nine points.

Case	Cycles, table flux	Years, table flux	Cycles, quadrupole flux	Ratio
$a=0$ , $10\to6.003$	$112\,981$	$1.417$	$112\,045$	$0.9917$
$a=0.9$ , $10\to6$	$260\,673$	$3.157$	$222\,719$	$0.8544$
$a=-0.9$ , $20\to10$	$971\,669$	$31.47$	$970\,777$	$0.9991$

Step 1: the frozen flux table

Fujita and Tagoshi’s Table VII gives the gravitational-wave luminosity, summed through $\ell=6$ , for circular equatorial orbits at seven spins from $-0.9$ to $0.9$ . The inputs file holds the rows $r=6$ to $r=20$ , transcribed by hand from arXiv version 1, with the caption quoted.

The table does not state its normalisation in a form we could quote, so the reading “energy flux to infinity per $(\mu/M)^2$ ” is our inference. Two checks support it. The mode-by-mode Table III at $r=6$ sums to the Table VII entry within $5\times10^{-5}$ . Far out, the table agrees with the Tagoshi–Sasaki post-Newtonian flux to infinity within $2\times10^{-4}$ at $r\ge50$ .

Reading: Fujita & Tagoshi (2004) · Tagoshi & Sasaki (1994)

Step 2: interpolate a smooth correction, not the flux

The flux falls by a factor of a thousand across the table, so the script interpolates the ratio of the table flux to the quadrupole flux instead. That ratio stays between 0.8 and 2 and changes smoothly with radius. The script fits a cubic through the four nearest tabulated radii and refuses any radius outside 6 to 20.

The energy slope $E'(r)$ comes from the exact Kerr circular energy by a fourth-order finite difference. Both integrals use Simpson’s rule in $\ln r$ with 2000 intervals.

t=\int_{\ln r_1}^{\ln r_0}\frac{r\,E'(r)}{\mathcal F(r)}\,d\ln r,\qquad N=\int_{\ln r_1}^{\ln r_0}\frac{\Omega_\phi(r)}{\pi}\,\frac{r\,E'(r)}{\mathcal F(r)}\,d\ln r,\qquad \mathcal F=\tfrac{32}{5}\,q\,\Omega_\phi^{10/3}\times\frac{\mathcal F_{\rm table}}{\mathcal F_{\rm quad}}

Integrating in

\ln r

spreads the points evenly where the integrand changes most.

Step 3: run and compare

Download the folder, check the checksums, run the script with the inputs file and a new output name, then compare your output with the expected one. On the machine that produced the expected file, Python 3.9.6 on macOS, a fresh run reproduces it byte for byte. Elsewhere, expect agreement to about one part in a trillion.

Time is converted to seconds with $GM_\odot/c^3=4.925490947\times10^{-6}\,\mathrm s$ and to years with the Julian year of 31,557,600 seconds. The wave frequency in the track is $f=2\Omega_\phi/(2\pi)$ , in hertz for $M=10^6\,M_\odot$ .

Step 4: break it on purpose

The output contains three deliberate changes and how much each moves the zero-spin cycle count. Linear instead of cubic interpolation of the correction changes it by −0.54%. Four Simpson intervals instead of 2000 change it by +0.088%. Replacing the Kerr energy with the Newtonian binding energy changes it by +99%.

The last number is the lesson: near the ISCO the binding energy matters as much as the flux. A fair test of a flux model has to keep the orbit exact.

What the result cannot tell you

Extrapolation is refused. Spin 0.9 has its ISCO at $2.32$ , far below the table’s first row, so that case stops at $r=6$ . The request for $r=4$ is listed in the output as refused.

Missing multipoles. The $\ell=6$ modes carry 0.29% of the zero-spin flux at $r=6$ , so the omitted $\ell\ge7$ modes are likely about 0.1% there. That is an extrapolation of the trend, not a calculation.

Missing physics. The horizon flux, eccentric and inclined orbits, resonances and every post-adiabatic term are absent. The post-adiabatic phase correction is of order one radian, so the cycle counts are right to the leading order in the mass ratio only.

Reading: Fujita & Tagoshi (2004) · Barack & Pound (2019)

Download and run

Put the six files in one folder, then run these two commands. The second writes your own output next to the expected one.

shasum -a 256 -c SHA256SUMS python3 recompute_frequency_evolution.py inputs.json my-output.json

recompute_frequency_evolution.py · The script, Python standard library only
inputs.json · Cases, units and the frozen flux table
expected-output.json · Expected output, Python 3.9.6
README.txt · Units, model and failure cases
SHA256SUMS · SHA-256 checksums of every file
LICENSE.txt · Permission and attribution

Try it in the Lab

Compare with the Lab’s quadrupole count

Choose zero spin, eccentricity zero and $p=10$ in the Lab.
Read the wave cycles left. The Lab rounds to three figures; compare it with 112,045 in the quadrupole column of the first row.

The Lab starts its count at the orbit size you set and uses the quadrupole flux only.

Exercises and explained answers

Which case shows the largest gap between the table and the quadrupole flux, and why that one?

Spin 0.9 from $r=10$ to $r=6$ , with a ratio of 0.854. On that range the quadrupole flux is 15% to 19% too large everywhere, so the error never changes sign and does not cancel. For zero spin the error changes sign near $r=8$ .

Change the zero-spin stop radius in the inputs to 5.9 and rerun. What happens?

The script stops with an error that names the radius and refuses to extrapolate, because $r=5.9$ lies below the first tabulated row. Even with a larger table it would also lie inside the zero-spin ISCO at $r=6$ , where no circular orbit is stable.

Primary sources and what we read in each

Barack & Pound (2019), Self-force and radiation reaction in general relativity, Rep. Prog. Phys., doi:10.1088/1361-6633/aae552
arXiv version 2, dated 19 November 2018. Section 5.2 and eqs. (72)–(74): the two-timescale phase expansion and what the adiabatic and post-adiabatic orders need; Section 5.3 on transient resonances; Section 6.1 and eqs. (77)–(78), the orbit-averaged balance laws, including negative horizon fluxes.
Fujita & Tagoshi (2004), New numerical methods to evaluate homogeneous solutions of the Teukolsky equation, Prog. Theor. Phys. 112, 415
arXiv version 1. Table VII, the gravitational-wave luminosity through l = 6 for circular equatorial orbits at seven spins; Table III, the mode-by-mode luminosity at r = 6M. The unit normalisation and the reading as flux to infinity are our inference, checked against the quadrupole and post-Newtonian limits.
Tagoshi & Sasaki (1994), Post-Newtonian expansion of gravitational waves from a particle in circular orbit around a Schwarzschild black hole, Prog. Theor. Phys. 92, 745
arXiv version 1. Eq. (43): the energy flux to infinity relative to the quadrupole value, through fourth post-Newtonian order; used only as an independent weak-field check of the flux table.

Keep reading

What this page links to directly

You are hereRecompute an EMRI frequency evolution from a published flux table

Experiments

Used in this articleEMRI Orbits Around a Spinning Black Hole

Articles

Extreme mass ratio inspirals · Reading path