Adiabatic binary black hole inspiral and waveform

In this example we will compute the fluxes of gravitational-wave energy for a point mass on a circular orbit in Schwarzschild spacetime, and use them to compute an adiabatic binary black hole inspiral and its associated waveform.

Installation

The first step is to install the Mathematica packages provided by the Black Hole Perturbation Toolkit. We recommend that you do so via the Paclet server using the instructions below.

Add the Black Hole Perturbation Toolkit paclet server to Mathematica’s list of paclet servers:

If[$VersionNumber≥12.1,
PacletSiteRegister["https://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"],
PacletSiteAdd["http://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"]
]

Get an updated list of packages available on the server:

If[$VersionNumber≥12.1,
PacletSiteUpdate["https://pacletserver.bhptoolkit.org"],
PacletSiteUpdate["http://pacletserver.bhptoolkit.org"]
]

Install the paclets:

PacletInstall["GeneralRelativityTensors"];
PacletInstall["KerrGeodesics"];
PacletInstall["SpinWeightedSpheroidalHarmonics"];
PacletInstall["ReggeWheeler"];
PacletInstall["Teukolsky"];
PacletInstall["PostNewtonianSelfForce"];

If a new version of a package is released and you would like to update, just run steps 2 and 3 above again and the latest available version will be installed.

Numerical Relativity waveform for comparison

We are going to produce a waveform for a binary black hole inspiral using an adiabatic self-force approximation, and compare it against a similar waveform produced in full Numerical Relativity using the Einstein Toolkit. The adiabatic approximation is valid in the limit where one of the black holes is much larger than the other, so we will choose a configuration where the black holes are very different in size. The public RIT waveform catalog has a 1:15 mass ratio configuration with quite a few gravitational-wave cycles, so let’s use that.

First, we load the gravitational wave strain for the l=2, m=2 spherical mode:

h22=;

This is a 1:15 mass ratio system with the individual masses, (Subscript[m, 1]) and (Subscript[m, 2]), of the two black holes as given in the metadata file:

m1=0.937508941265577;
m2=0.0624981715324337;

The total mass, M, is the sum of the two constituent masses:

M=m1+m2;

We now define the (small) mass ratio ϵ, the large mass ratio q, and the symmetric mass ratio, ν:

ϵ=m2/m1;
q=1/ϵ;
ν=((m1 m2)/(((m1+m2))^2));
{q,ϵ,ν}

We will set up our adiabatic calculation to match this system as closely as possible. To do so we will need to provide initial conditions for the adiabatic equations, namely the initial orbital frequency and phase. One way we obtain the from the numerical relativity data is by measuring the phase and frequency of the (2,2) mode of the waveform and dividing by 2. If you have SimulationTools installed you can do that using the following code.

Needs["SimulationTools`"];
Ωinit=-(1/2)Frequency[ToDataTable[h22]][[1]];
ϕinit=-(1/2)Phase[ToDataTable[h22]][[1]];

Doing so, we get the values below:

Ωinit=0.02810756207135854`;

ϕinit=0.004042174353145999`;

These correspond to having an initial separation given by:

r0init=(((M Ωinit))^(-2/3))

Adiabatic inspiral

In the adiabatic approximation,  we effectively treat the system as a point particle in orbit around a black hole. The orbit will gradually decay due to the loss of energy in the form of gravitational radiation. We can think of this as the point particle not quite following a geodesic of the black hole spacetime, but rather as having a small acceleration away from a geodesic. It therefore satisfies the equations for an accelerated particle,



where the force term on the right hand side is the gravitational self-force.

We are interested in a system undergoing a quasicircular inspiral, in which case the t component  of this equation gives us everything we need:

At leading order (ignoring the force term) we get a relationship between the frequency and radius of the orbit,

The adiabatic corrections from the gravitational self-force tell us that the orbit is not a constant, but slowly evolves according to

Conservation of energy tells us that any energy lost from the orbit must exactly balance the energy emitted in gravitational waves, and this energy balance law allows us to write the gravitational self-force in terms of the flux of gravitational-wave energy infinitely far from the system,

(F^t)=ℱ (dt/(dτ))

Putting this together, and using the chain rule we get an evolution equation for the frequency, Ω:

((dΩ)/dt)=((3((((Subscript[r, 0])-3M))^(3/2)))/((M^(3/2))(Subscript[r, 0])((Subscript[r, 0])-6M)))ℱ.

This is a simple ordinary differential equation that we can solve for Ω provided we have the flux, ℱ. We can then integrate this once more to get the orbital phase, (Subscript[ϕ, 0]) using

((Subscript[dϕ, 0])/dt)= Ω

Finally, we can then write the waveform in a form where it is factored into an amplitude and orbital phase.



(Note that since the orbital phase can differ slightly from the waveform phase, the amplitudes are still complex numbers.)

The key things we need, then, are the flux of gravitational waves, ℱ, and the amplitudes, (h^((1))) for a point particle on a circular orbit. The packages in the Toolkit make it very easy to compute these in the frequency domain using the Teukolsky formalism.

Compute energy flux and asymptotic waveform amplitudes

We can use the Teukolsky package to numerically compute the solutions to the Teukolsky equation. First, we load the package

<<Teukolsky`

As an example, let’s consider the (l,m)=(2,2) mode of (Subscript[ψ, 4]) for a circular orbit of radius (Subscript[r, 0])=10 about a Schwarzschild black hole.  First, we set the parameters of the case we are considering

𝓈=-2;ℓ=2;𝓂=2;

a=0; (* Schwarzschild spacetime *);
e=0;(* Zero eccentricity *);
θinc=0; (* Equatorial orbit *)
n=0;k=0; (* Circular orbits don't have frequencies related to radial or polar motion *)
p=10.0`32;(* Semi-latus rectum is just the orbital radius in this case *)

We can now use the KerrGeodesics package to compute the orbital frequency:

Ω=KerrGeoFrequencies[a,p,e,Cos[θinc]]["(Subscript[Ω, ϕ])"];

In the circular orbit case the frequency of the spherical harmonic mode is a straightforward multiple of this:

ω=𝓂 Ω;

Now, we can use the TeukolskyPointParticleMode function from the Teukolsky package to compute the solution to the Teukolsky equation in this case

orbit=KerrGeoOrbit[a,p,e,Cos[θinc]];

ψ4=TeukolskyPointParticleMode[𝓈,ℓ,𝓂,0,0,orbit]

From this, we can immediately obtain the contribution to the flux from this mode

ℱ22r10=ψ4["EnergyFlux"]

We can also obtain the metric perturbation amplitudes from the amplitudes of (Subscript[ψ, 4]) from the simple relation (Overscript[h, ..])=(Subscript[ψ, 4])

h22r10=((2ψ4["Amplitudes"]["ℐ"])/(ω^2));

Now, we could clearly run this for a range of modes and orbital radii, but in interest time let’s load some precomputed data for the (2,2) mode of the amplitude and the total (summed over modes) flux on a grid of (Subscript[r, 0]) values in the range (Subscript[r, 0])∈[6,30].

ℱ=Interpolation[,InterpolationOrder->8];

hAmp=Interpolation[,InterpolationOrder->8];

h22r10

hAmp[10]

ℱ[10]

2ℱ22r10["ℐ"]

Compute Inspiral

We are now ready to compute an inspiral by solving the equations above for Ω(t) and (Subscript[ϕ, 0])(t) and (Subscript[r, 0])(t)

{Ωsol,ϕsol,r0sol}=NDSolveValue[{
Ω0'[t]==((3 (((-3 M+r0[t]))^(3/2)) ν ℱ[r0[t]] )/((M^(3/2)) r0[t] (-6 M+r0[t]))),
ϕ0'[t]==Ω0[t],
Ω0[t]==Sqrt[(M/((r0[t])^3))],
r0[0]==r0init,
ϕ0[0]==ϕinit,
WhenEvent[r0[t]≤6.01M, "StopIntegration"]},{Ω0,ϕ0,r0},{t,0,∞},PrecisionGoal->10,AccuracyGoal->10];

{tmin,tmax}=r0sol["Domain"][[1]];

Plot solutions

Evolution of (Subscript[r, 0])

Plot[r0sol[t],{t,tmin,tmax},PlotTheme->"Detailed",ImageSize->800]

Evolution of Ω

Plot[Ωsol[t],{t,tmin,tmax},PlotTheme->"Detailed",ImageSize->800]

Evolution of (Subscript[ϕ, 0])

Plot[ϕsol[t],{t,tmin,tmax},PlotTheme->"Detailed",ImageSize->800]

Waveform

Construct waveform from amplitude and phase

We now have everything we need to construct the waveform, we just need to multiply the amplitude by the phase

h[t_]:=ν hAmp[r0sol[t]] (E^(-I (𝓂 ϕsol[t]+Arg[hAmp[r0sol[0]]])));

Let’s plot our waveform alongside the waveform from the Numerical Relativity simulation

Show[ListLinePlot[Re[h22],PlotTheme->"Detailed",FrameLabel->{"t","Re[(Subscript[h, 2, 2])]"},PlotStyle->Black],Plot[Re[h[t]],{t,tmin,tmax}],PlotRange->All]