GW fluxes and waveforms from self-force theory

Author: Zachary Nasipak

Theoretical background

In self-force theory, the spacetime metric of a binary system with masses  is first expanded in the small mass-ratio 

	)

At leading-order the motion of the small body is simply a geodesic in the background spacetime with constants of motion 

	- : 	orbital energy
    	- : 	z-component of the angular momentum
          - : 	Carter constant

and fundamental frequencies 

	- : 	radial frequency for motion between periapsis  and apoapsis 
    	- : 	polar frequency for motion between the maximum and minimum inclination angles  and 
	- : 	azimuthal frequency for motion about the spin-axis of 

Beyond leading-order, the constants of motion (and the fundamental frequencies) are no longer constant, but are slowly-evolving mechanical parameters due to dissipation via gravitational radiation. Therefore, to solve for the motion, one can perform a multiscale expansion, in which one separates the dynamics into slowly-evolving quantities  and fast-evolving orbital / GW phases . In this multiscale self-force framework, the equations of motion for the binary inspiral then take the form

	

	

where  represents the “constants of motion” or mechanical parameters such as energy and angular momentum, while the  represent the three components of the Fourier-decomposed waveform phase. Including just the leading-order terms on the righthand-side produces the so-called adiabatic (or 0PA) approximation, while including up to the second terms produces the first post-adiabatic (or 1PA) approximation.

The inspiral waveform then takes the form

	

where  are indices labeling the Fourier harmonics of the waveform. For example, for tri-periodic orbits (i.e., eccentric, precessing orbits),  and  so that  is the discrete Fourier spectrum of the waveform.

In self-force theory, the spacetime metric of a binary system with masses  is first expanded in the small mass-ratio 

	)

At leading-order the motion of the small body is simply a geodesic in the background spacetime with constants of motion 

	- : 	orbital energy
    	- : 	z-component of the angular momentum
          - : 	Carter constant

and fundamental frequencies 

	- : 	radial frequency for motion between periapsis  and apoapsis 
    	- : 	polar frequency for motion between the maximum and minimum inclination angles  and 
	- : 	azimuthal frequency for motion about the spin-axis of 

Beyond leading-order, the constants of motion (and the fundamental frequencies) are no longer constant, but are slowly-evolving mechanical parameters due to dissipation via gravitational radiation. Therefore, to solve for the motion, one can perform a multiscale expansion, in which one separates the dynamics into slowly-evolving quantities  and fast-evolving orbital / GW phases . In this multiscale self-force framework, the equations of motion for the binary inspiral then take the form

	

	

where  represents the “constants of motion” or mechanical parameters such as energy and angular momentum, while the  represent the three components of the Fourier-decomposed waveform phase. Including just the leading-order terms on the righthand-side produces the so-called adiabatic (or 0PA) approximation, while including up to the second terms produces the first post-adiabatic (or 1PA) approximation.

The inspiral waveform then takes the form

	

where  are indices labeling the Fourier harmonics of the waveform. For example, for tri-periodic orbits (i.e., eccentric, precessing orbits),  and  so that  is the discrete Fourier spectrum of the waveform.

In this tutorial, we demonstrate how we can construct all of the quantities that go into a 0PA inspiral waveform.

Setup and quick guide

Installing the Black Hole Perturbation Toolkit

Installation instructions for Mathematica packages on the Black Hole Perturbation Toolkit: https://bhptoolkit.org/mathematica-install.html

PacletSiteRegister["https://pacletserver.bhptoolkit.org", "Black Hole Perturbation Toolkit Paclet Server"];
PacletSiteUpdate["https://pacletserver.bhptoolkit.org"];

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

<<KerrGeodesics`
<<SpinWeightedSpheroidalHarmonics`
<<Teukolsky`
<<PostNewtonianSelfForce`

A quick guide

Compute a circular orbit with radius  about a Kerr black hole spin parameter

orbit=KerrGeoOrbit[0.9,30.,0,1];

Solve the Teukolsky equation (for several modes)

l$max=6;
Do[teuk$Solutions[l,m]=TeukolskyPointParticleMode[-2,l,m,0,0,orbit],{l,2,l$max},{m,Max[l-4,1],l}];//AbsoluteTiming

Compute fluxes from the Teukolsky solutions

flux=Sum[2teuk$Solutions[l,m]["Fluxes"]["Energy"],{l,2,l$max},{m,Max[l-4,1],l}];

Compare fluxes with Post-Newtonian self-force results

flux$PN$ℐ$Series = PostNewtonianExpansion["/Kerr/Circular/Flux/Energy/Infinity"]["Series"];
flux$PN$ℋ$Series = PostNewtonianExpansion["/Kerr/Circular/Flux/Energy/Horizon"]["Series"];

y$val = 1/orbit["p"];
logy$val=Log[y$val];
a$val=orbit["a"];
flux$PN$ℐ$value =Re@ReleaseHold[Normal[flux$PN$ℐ$Series]]/.y->y$val/.Logy->logy$val/.a->a$val;
flux$PN$ℋ$value =Re@ReleaseHold[Normal[flux$PN$ℋ$Series]]/.y->y$val/.Logy->logy$val/.a->a$val;

Print["Num flux at infinity: ",flux["ℐ"], "\nPN flux at infinity:  ",flux$PN$ℐ$value];
Print["Num flux at infinity: ",flux["ℋ"], "\nPN flux at infinity:  ",flux$PN$ℋ$value];

Build waveform from the Teukolsky solutions as a function of retarded time  ( to be defined later) and viewing polar angle

u$Signal=Range[0,2000,1];
viewingAngle=Pi/3;
GW$strain =Sum[
ω$tmp=teuk$Solutions[l,m]["ω"];
h1lm[l,m]=-2teuk$Solutions[l,m]["Amplitudes"]["ℐ"]/ω$tmp^2 ;
Sslm[l,m]=teuk$Solutions[l,m]["AngularFunction"][viewingAngle,0];h1lm[l,m]Sslm[l,m]Exp[-I ω$tmp u$Signal],{l,2,l$max},{m,Max[l-4,1],l}];
GW$strain=GW$strain+Sum[
ω$tmp=teuk$Solutions[l,m]["ω"];
h1lm[l,-m]=-2(-1)^l Conjugate[teuk$Solutions[l,m]["Amplitudes"]["ℐ"]]/ω$tmp^2 ;
Sslm[l,-m]=SpinWeightedSpheroidalHarmonicS[-2,l,-m,-ω$tmp*orbit["a"]][viewingAngle,0];h1lm[l,-m]Sslm[l,-m]Exp[I ω$tmp u$Signal],{l,2,l$max},{m,Max[l-4,1],l}];

Inspect GW waveform strain (normalized by the distance to the source)

ListPlot[{Re[GW$strain],-Im[GW$strain]},
PlotTheme->"Scientific",
FrameLabel->{"u/M","(strain) x (distance)"},
ImageSize->Medium]

“Background” geodesic

We can compute the geodesic motion of a particle  in Boyer-Lindquist coordinates  using the KerrGeodesics package

<<KerrGeodesics`

KerrGeoOrbit is the main function for constructing Kerr geodesics. Periodic, bound, stable orbits are parametrized in terms of the Keplerian-like parameters:

	- : 	Dimensionless Kerr spin parameter
	- : 	Semi-latus rectum
	- : 	Orbital eccentricity
	- : 	Cosine of the orbital inclination

along with the masses  and , where  is the Kerr mass parameter. For numerical calculations, we only output dimensionless quantities.
KerrGeoOrbit takes these parameters as inputs

{a$val,p$val,e$val,xI$val} = {0.9,6.2,0.4,0.8};
orbit = KerrGeoOrbit[a$val,p$val,e$val,xI$val]

We can see that the function stores useful input, such as the orbital parameters. It also tells us how the orbital solution are parametrized in time. By default, solutions are given as functions of Mino time , because this is the time parameter that decouples the equations of motion. If the motion is equatorial ( ) then one can alternatively parametrize the motion in terms of the Darwin variable  defined by .

By calling Keys we can see what information is easily accessible in the KerrGeoOrbitFunction

Keys[orbit]

For example we can call for the constants of motion—energy, angular momentum, and Carter constant—and frequencies described earlier.

orbit["ConstantsOfMotion"]
orbit["Frequencies"]

Note, however, that the frequencies are not automatically the frequencies with respect to coordinate time but Mino time. To get coordinate time frequencies, we need to divide by

<|"(Subscript[Ω, r])"->"(Subscript[Υ, r])"/"(Subscript[Υ, t])","(Subscript[Ω, θ])"->"(Subscript[Υ, θ])"/"(Subscript[Υ, t])","(Subscript[Ω, ϕ])"->"(Subscript[Υ, ϕ])"/"(Subscript[Υ, t])"|>/.orbit["Frequencies"]//N

Alternatively, we can call the KerrGeoFrequencies function

KerrGeoFrequencies[a$val,p$val,e$val,xI$val]

This function also makes other frequencies readily available using the “Time” option

KerrGeoFrequencies[a$val,p$val,e$val,xI$val,"Time"->"Mino"]
KerrGeoFrequencies[a$val,p$val,e$val,xI$val,"Time"->"Proper"]

Similarly, we do not need to call KerrGeoOrbit to get the constants of motion. We can simply call KerrGeoConstantsOfMotion or the individual functions KerrGeoEnergy, KerrGeoAngularMomentum, and KerrGeoCarterConstant

KerrGeoConstantsOfMotion[a$val,p$val,e$val,xI$val]
<|"ℰ"->KerrGeoEnergy[a$val,p$val,e$val,xI$val],"ℒ"->KerrGeoAngularMomentum[a$val,p$val,e$val,xI$val], "𝒬"->KerrGeoCarterConstant[a$val,p$val,e$val,xI$val]|>

Furthermore, the KerrGeoOrbitFunction contains the trajectories  as functions of the time parametrization (typically Mino time). We can evaluate the position of the small body at Mino time λ by simply calling the function with our time value

Tr$val=2Pi/orbit["RadialFrequency"];
λ$vals=Range[0,5Tr$val,5Tr$val/1000];
{tp$vals,rp$vals,θp$vals,ϕp$val}=orbit[λ$vals];

rp$of$t$vals=Transpose[{tp$vals,rp$vals}];
zp$of$t$vals=Transpose[{tp$vals,Cos[θp$vals]}];
ListPlot[{rp$of$t$vals,zp$of$t$vals},
PlotTheme->"Scientific",
Joined->True]

Note that by default, the orbit is begins at

orbit["InitialPhases"]

We can change the initial conditions by changing the initial phases passed to KerrGeoOrbit

orbit2 = KerrGeoOrbit[a$val,p$val,e$val,xI$val,{0,Pi,Pi/2,2.3}];

Tr$val2=2Pi/orbit2["RadialFrequency"];
λ$vals2=Range[0,5Tr$val2,5Tr$val2/1000];
{tp$vals2,rp$vals2,θp$vals2,ϕp$val2}=orbit2[λ$vals2];

rp$of$t$vals2=Transpose[{tp$vals2,rp$vals2}];
zp$of$t$vals2=Transpose[{tp$vals2,Cos[θp$vals2]}];
ListPlot[{rp$of$t$vals2,zp$of$t$vals2},
PlotTheme->"Scientific",
Joined->True]

Important to keep in mind that the functions always return the dimensionless versions of all these quantities. Therefore if one wants the dimensionful frequencies, we would need to multiply by (M^(-1)), while for the constants of motion

	 ℰ		 ℒ		 𝒬

Solving the Teukolsky equation (in the frequency domain)

Background & notation

In Boyer-Lindquist coordinates , the Teukolsky equation takes the form

    

where ,  is the (tetrad-projected) source, and the spin-weighted Teukolsky operator is given by,

	  
    
with , and

	-   for scalar perturbations;
	-   and  ; and
	-   and  for gravitational perturbations.

Transforming to the frequency domain, the fields and sources can be mode-decomposed into

	 ,

	 ,

with , leading to separated ODEs

	,

and

	

with eigenvalue  and spheroidicity .

When the source motion is bound and (multi-)periodic, the integral over frequency reduces to a sum over the discrete frequency-spectrum of the motion. For the geodesics we consider, the discrete spectrum is given by

The angular Teukolsky equation - spin weighted spheroidal harmonics

The solutions to the angular Teukolsky equation 



are known as spin-weighted spheroidal harmonics. They are generalizations of spin-weighted spherical harmonics 

, 

which are spin-weighted generalizations of (scalar) spherical harmonics . Therefore, the spheroidal eigenvalue  reduces to the spin-weighted spherical eigenvalue when the spheroidicity is zero:

.

We can access all of these quantities with the SpinWeightedSpheroidalHarmonics package

<<SpinWeightedSpheroidalHarmonics`

Let’s plot the eigenvalue  as a function of  for spin-weight -2 (associated with gravitational perturbations), and harmonic indices . We can also compare this to the spherical eigenvalue . We can access this quantity by calling SpinWeightedSpheroidalEigenvalue[s, j, m, γ]

{s,j,m}={-2,3,2};
plot1=Plot[SpinWeightedSpheroidalEigenvalue[s,j,m,γ],{γ,-10,10},
PlotTheme->"Scientific"];
plot2=Plot[j(j+1)-s(s+1),{γ,-10,10},
PlotTheme->"Scientific",
PlotStyle->{Gray,Dashed}];
Show[plot1,plot2,
ImageSize->Medium,
FrameLabel->{γ,λ}]

As expected, the spheroidal eigenvalue reduces to the spherical case at γ=0.

Furthermore, we can evaluate these functions for complex values of the spheroidicity. Let’s look at the behavior of the complex argument as we move through the complex plain

λ$ComplexSpheroidicity=Flatten[Table[{γR+γI I,SpinWeightedSpheroidalEigenvalue[s,j,m,N[γR+I γI,20]]},{γR,0,5,0.1`20},{γI,-2,2,0.1`20}],1];//AbsoluteTiming

ListDensityPlot[λ$ComplexSpheroidicity/.{x_,y_}:>{Re[x],Im[x],Arg[y]}]

Angular solutions, i.e. spin-weighted spheroidal harmonics, are accessed through SpinWeightedSpheroidalHarmonicS which takes the arguments:
	
	- 	: spin-weight
	- 	: spheroidal harmonic index
	- 	: azimuthal harmonic index
	- 	: spheroidicity

So let’s compute the harmonics for a few different harmonic modes

γ = 0.4;
Clear[Ssjm];
Do[Ssjm [s,j$,m$]= SpinWeightedSpheroidalHarmonicS[s,j$,m$,γ],{j$,Abs[s],4},{m$,-j$,j$}];

Plotting these functions we can see some of the symmetries of functions. For instance the function vanishes at the poles unless  and we have various parity symmetries as . Looking at the  case

j$vals={2,3};
m$vals={-2,2};
jm$labels=Flatten[Table[ToString[{j$,m$}],{j$,j$vals},{m$,m$vals}],1];
Plot[Evaluate[Join[Table[Ssjm[s,j$,m$][θ,0],{j$,j$vals},{m$,m$vals}]]]
,{θ,0,Pi},
PlotTheme->"Scientific",
FrameLabel->{θ,(Subscript[S, sjmγ])[θ]},PlotLegends->jm$labels]

And then for the

j$vals={2,3};
m$vals={-1,1};
jm$labels=Flatten[Table[ToString[{j$,m$}],{j$,j$vals},{m$,m$vals}],1];
Plot[Evaluate[Join[Table[Ssjm[s,j$,m$][θ,0],{j$,j$vals},{m$,m$vals}]]]
,{θ,0,Pi},
PlotTheme->"Scientific",
FrameLabel->{θ,(Subscript[S, sjmγ])[θ]},PlotLegends->jm$labels]

It is also straightforward to compute derivatives of the harmonic functions using Mathematica’s traditional derivative operators. For example, let’s plot the second polar derivative using Derivative[2,0][SpinWeightedSpheroidalHarmonic[s,j,m,γ]][θ,0]

j$vals={2,3};
m$vals={-1,1};
jm$labels=Flatten[Table[ToString[{j$,m$}],{j$,j$vals},{m$,m$vals}],1];
Plot[Evaluate[Join[Table[Derivative[2,0][Ssjm[s,j$,m$]][θ,0],{j$,j$vals},{m$,m$vals}]]]
,{θ,0,Pi},
PlotTheme->"Scientific",
FrameLabel->{θ,(Subscript[S, sjmγ])''[θ]},PlotLegends->jm$labels]

Furthermore, there are multiple methods for calculating the spin-weighted spheroidal harmonics. The default is “Leaver”, which is based on Leaver’s (1985) analytic method. Another option is “SphericalExpansion” in which the spin-weighted spheroidal harmonics are represented as series of spin-weighted spherical harmonics

	

with spherical-spheroidal coupling/mixing coefficients . This is a formally infinite sum that must be truncated in numerical calculations. Therefore, to use this method we also need to supply the “NumTerms” option

Ssjm$expansion=SpinWeightedSpheroidalHarmonicS[-2,2,2,γ,Method->{"SphericalExpansion","NumTerms"->20}]

Calling Keys to see the attributes of this function

Ssjm$expansion//Keys

The output function now contains an attribute called “ExpansionCoefficients” which includes all of the coefficients   that were calculated. We can see how quickly these coefficients decay, and we see the series converges quite rapidly

ListLogPlot[Abs@Ssjm$expansion["ExpansionCoefficients"],
PlotTheme->"Scientific",
FrameLabel->{"\!\(\*TemplateBox[<|"boxes" -> FormBox[
StyleBox["l", "TI"], TraditionalForm], "errors" -> {}, "input" -> "l", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)","\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubsuperscriptBox[
StyleBox["b", "TI"], 
RowBox[{
StyleBox["s", "TI"], 
StyleBox["j", "TI"], 
StyleBox["m", "TI"], "γ"}], 
StyleBox["l", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "b^l_{sjm\\gamma}", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"}]

The radial Teukolsky equation - homogeneous solutions

The radial Teukolsky equation for parameters:

	- 		: spin-weight 
	- 	: harmonic numbers
	-  		: frequency
	-  		: Kerr spin

is given by



where ,  ,    is the spheroidal eigenvalue (separation constant), and  is the radial decomposition of the source. We construct these solutions with the Teukolsky package.

<<Teukolsky`

Let’s focus first on the homogeneous equation .

By default, this package constructs the “In” and “Up” homogeneous solutions  and , which correspond to the asymptotic boundary conditions

    ,

    ,

    ,

    ,

where , , and the tortoise coordinate is given by the differential relation , which can be immediately integrated, leading to

,

where  and .

So let’s build those solutions! First let’s define our radial Teukolsky parameters

{s$val,j$val,m$val,a$val,ω$val}={-2,0.9`40,4,2,0.5`40};

To construct solutions, we use the TeukolskyRadial[s,j,m,a,ω] function

R$teuk=TeukolskyRadial[s$val,j$val,m$val,a$val,ω$val]

We see it returns an Association with the Keys corresponding to the “In” and “Up” homogeneous solutions

RIn=R$teuk["In"];
RUp=R$teuk["Up"];

We can now inspect the attributes of these solutions

RUp//Keys

Let’s inspect some of the attributes, namely “BoundaryConditions”, “Method”, and “Domain”

{RIn["BoundaryConditions"],RUp["BoundaryConditions"]}
{RIn["Method"],RUp["Method"]}
{RIn["Domain"],RUp["Domain"]}

“BoundaryConditions” reminds us of what boundary conditions each homogeneous solution satisfies, “Method” tells us that the solver used the MST method to construct the solutions, and “Domain” reports that both solutions are valid over the full physical radial domain . (The MST method refers to a semi-analytic representation of the homogeneous solutions found by Mano, Suzuki, and Takasugi (1996).)

We can also look at the scattering amplitudes of the homogeneous solutions, as defined above:

	,

RUp["Amplitudes"]

We can also plot the solutions. Just as an example, let’s plot just the “Up” solution. First evaluate the MST solutions

r$vals=Range[10., 150, 0.5];
RUp$vals=RUp[r$vals];//AbsoluteTiming

Then we also construct the asymptotic solution by defining

RStarOfR[a_,r_]:=With[{κ=Sqrt[1-a^2]},r+(1+κ)/κ Log[(r-(1+κ))/2]-(1-κ)/κ Log[(r-(1-κ))/2]];
RUp$asymp$vals=r$vals^-(2 s$val+1)Exp[I ω$val RStarOfR[a$val,r$vals]];

Now let’s plot the two functions

plotData={
Transpose[{r$vals,Re@RUp$vals}],
Transpose[{r$vals,Re@RUp$asymp$vals}]
};
ListPlot[plotData,
Joined->True,
PlotRange->All,
PlotStyle->{Automatic,Dashed},
PlotTheme->"Scientific",
FrameLabel->{"r/M","Re[(R^up)]"}]

Note that while the MST method has the advantage of being semi-analytic, it can be quite slow when evaluating the solutions at many points. Therefore, one can instead construct solutions via the “NumericalIntegration” method. This uses MST solutions to implement boundary conditions and uses Mathematica’s ODE solvers to compute points between the boundary points.

R$teuk$int=TeukolskyRadial[s$val,j$val,m$val,a$val,ω$val,"Method"->{"NumericalIntegration","Domain"->{"In"->{2,20},"Up"->{4,200}}}];//AbsoluteTiming
RUp$vals$int=R$teuk$int["Up"][r$vals];//AbsoluteTiming

Now we can look at the relative difference of the two methods

plotData=Transpose[{
r$vals,
Abs[1-RUp$vals$int/RUp$vals]
}];
ListPlot[plotData,
Joined->False,
PlotRange->All,
PlotTheme->"Scientific",
FrameLabel->{"r/M","abs error"}]

The radial Teukolsky equation - inhomogeneous amplitudes

Next we turn to inhomogeneous solutions with a point particle source .

Rather than constructing the inhomogeneous solutions over the full domain, we typically construct inhomogeneous solutions just in vacuum regions, i.e.,  and . 

In these regions, the solution has the asymptotic form





where  are known as the Teukolsky amplitudes, which define the complex amplitude of the inhomogeneous solution at the horizon (ℋ) and infinity (ℐ). These amplitudes are important for computing the first-order amplitudes of the waveform and the leading-order (0PA) fluxes.

To get these amplitudes, we first define the source using the KerrGeodesics package. So let’s consider a small body on an orbit parametrized by .

orbit=KerrGeoOrbit[0.9,6.2,0.4,0.8];

Recall that for bound motion, the frequency-domain representation of the field reduces to a discrete sum over frequency modes

	
	
Therefore, we can construct Teukolsky amplitudes on a (frequency) mode-by-mode basis via TeukolskyPointParticleMode[s,j,n,k,orbit]. Modes are labeled in the order . Below we compute the  for .

AbsoluteTiming[Ψ$teuk=TeukolskyPointParticleMode[-2,2,2,0,-1,orbit]]

We can look at the attributes of this output

Ψ$teuk//Keys

Immediately we see it contains information regarding the frequency of this mode, eigenvalue, radial and angular homogeneous solutions described in the preceding sections, the Teukolsky amplitudes. Let’s look at the amplitudes ,  the frequency, along with the homogeneous radial and angular mode functions

Ψ$teuk["Amplitudes"]
Ψ$teuk["ω"]
Ψ$teuk["RadialFunctions"]
Ψ$teuk["AngularFunction"]

Additionally, we see that there are attributes related to the fluxes. This is because the 0PA energy and angular momentum fluxes can be written in terms of the Teukolsky amplitudes:

	
	
The first term in the sum gives the flux through the horizon (“ℋ”), while the second is flux at infinity (“ℐ”)

Ψ$teuk["Fluxes"]

THE DISCLAIMER

Precision, precision, precision. The Teukolsky amplitudes are constructed via the typical Greens function/variation of parameters approach, where the homogeneous solutions are integrated against the source. As the mode numbers increase, these source integrals become highly oscillatory and are prone to precision loss.

To see this in action, let’s first confirm the accuracy of our earlier calculation by recalculating it with twice as much numerical precision.

orbit$highPrec30=KerrGeoOrbit[0.9`30,6.2`30,0.4`30,0.8`30];
AbsoluteTiming[Ψ$teuk$highPrec30=TeukolskyPointParticleMode[-2,2,2,0,-1,orbit$highPrec30];]

We see the higher precision increases the computational time

Ψ$teuk$highPrec30["Amplitudes"]
Ψ$teuk$highPrec30["Fluxes"]
Print["Relative differences ",Abs[1- %%/Ψ$teuk["Amplitudes"]]]
Print["Relative differences ",Abs[1- %%/Ψ$teuk["Fluxes"]]]

It also highlights that our horizon solutions were only accurate to about 5-6 digits.

Now let’s repeat this experiment but with the mode number . Let’s run it at several different precision levels. We will notice quite a few warnings from Mathematica

Ψ$teuk$Unstable=TeukolskyPointParticleMode[-2,18,10,-2,8,orbit];
Ψ$teuk$Unstable$highPrec30=TeukolskyPointParticleMode[-2,18,10,-2,8,orbit$highPrec30];
orbit$highPrec40=KerrGeoOrbit[0.9`40,6.2`40,0.4`40,0.8`40];
Ψ$teuk$Unstable$highPrec40=TeukolskyPointParticleMode[-2,18,10,-2,8,orbit$highPrec40];
orbit$highPrec50=KerrGeoOrbit[0.9`50,6.2`50,0.4`50,0.8`50];
Ψ$teuk$Unstable$highPrec50=TeukolskyPointParticleMode[-2,18,10,-2,8,orbit$highPrec50];
orbit$highPrec100=KerrGeoOrbit[0.9`100,6.2`100,0.4`100,0.8`100];
Ψ$teuk$Unstable$highPrec100=TeukolskyPointParticleMode[-2,18,10,-2,8,orbit$highPrec100];

Now let’s compare the Teukolsky amplitudes across all of these calculations

Ψ$teuk$Unstable["Amplitudes"]
Ψ$teuk$Unstable$highPrec30["Amplitudes"]
Ψ$teuk$Unstable$highPrec40["Amplitudes"]
Ψ$teuk$Unstable$highPrec50["Amplitudes"]
Ψ$teuk$Unstable$highPrec100["Amplitudes"]

We see the machine precision result does not agree with any of the extended precision results. 30 digits of precision totally fails internally. 40 digits leads to a completely new answer. Finally as we go beyond 50 digits, we seem to asymptote to some consistent result, though it is misleading because the integration routine reports a relative error of 1. So for higher mode numbers, we need to be careful with the amount of precision we use or whether we should trust the result.

Finally, we can then “extend” these solutions into the source region, i.e. , using the so-called method of extended homogeneous solutions for a point-particle source on a bound periodic geodesic,

    ,

    ,

where we have the mode frequencies , the discrete spheroidicity , and the extended homogeneous radial solutions

“Snapshot” waveforms and fluxes

Now we have all of the ingredients we need to construct our waveforms. Recall at  0PA and 1SF, the waveform takes the form

	
	

where the the evolution of the phases and amplitudes dictated by the 0PA equations of motion

	,

	

From the previous sections, we saw that we can get the leading-order frequencies  from the KerrGeodesics package and we can compute the 0PA fluxes  using the Teukolsky amplitudes in the Teukolsky package.  (Note at this order we can equate  to the orbital Keplerian elements  via their geodesic relations, so  and likewise for fluxes.) We can also compute the first-order amplitudes using the asymptotic relation

	
	
Therefore, the waveform amplitudes are related to the Teukolsky amplitudes by

	  	(spherical harmonic basis)
	
				(spheroidal harmonic basis)


where we need to account for the fact that we usually want to produce waveforms on a spherical harmonic basis, but Kerr perturbations naturally separate on a spheroidal harmonic basis. The coupling coefficients help (re-)project between the two bases.

We now have everything we need to calculate waveforms. For simplicity, we will look at “snapshot” waveforms and fluxes. This is what we get from setting  in the equations of motion; in other words, we ignore the effects of radiation reaction. In fact, this is what we initially calculated in our opening quick guide. So let’s return to that example now that we  understand the internal mechanics of the calculation.

A quick guide - revisited

Compute a circular orbit with radius  about a Kerr black hole spin parameter

orbit=KerrGeoOrbit[0.99,5.,0,1];

Solve the Teukolsky equation (for several modes)

l$max=8;
Do[
teuk$Solutions[l,m]=TeukolskyPointParticleMode[-2,l,m,0,0,orbit],
{l,2,l$max},
{m,Max[l-4,1],l}];//AbsoluteTiming

Compute fluxes from the Teukolsky solutions

flux=Sum[
2teuk$Solutions[l,m]["Fluxes"]["Energy"],
{l,2,l$max},
{m,Max[l-4,1],l}];

Compare fluxes with Post-Newtonian self-force results

flux$PN$ℐ$Series = PostNewtonianExpansion["/Kerr/Circular/Flux/Energy/Infinity"]["Series"];
flux$PN$ℋ$Series = PostNewtonianExpansion["/Kerr/Circular/Flux/Energy/Horizon"]["Series"];

y$val = 1/orbit["p"];
logy$val=Log[y$val];
a$val=orbit["a"];
flux$PN$ℐ$value =Re@ReleaseHold[Normal[flux$PN$ℐ$Series]]/.y->y$val/.Logy->logy$val/.a->a$val;
flux$PN$ℋ$value =Re@ReleaseHold[Normal[flux$PN$ℋ$Series]]/.y->y$val/.Logy->logy$val/.a->a$val;

Print["Num flux at infinity: ",flux["ℐ"], "\nPN flux at infinity:  ",flux$PN$ℐ$value];
Print["Num flux at horizon: ",flux["ℋ"], "\nPN flux at horizon:  ",flux$PN$ℋ$value];

Build waveform from the Teukolsky solutions as a function of retarded time  and viewing polar angle

u$Signal=Range[0,200,1];
viewingAngle=Pi/3;
GW$strain =Sum[
ω$tmp=teuk$Solutions[l,m]["ω"];
h1lm[l,m]=-2teuk$Solutions[l,m]["Amplitudes"]["ℐ"]/ω$tmp^2 ;
Sslm[l,m]=teuk$Solutions[l,m]["AngularFunction"][viewingAngle,0];h1lm[l,m]Sslm[l,m]Exp[-I ω$tmp u$Signal],
{l,2,l$max},
{m,Max[l-4,1],l}];
GW$strain=GW$strain+Sum[
ω$tmp=teuk$Solutions[l,m]["ω"];
h1lm[l,-m]=-2(-1)^l Conjugate[teuk$Solutions[l,m]["Amplitudes"]["ℐ"]]/ω$tmp^2 ;
Sslm[l,-m]=SpinWeightedSpheroidalHarmonicS[-2,l,-m,-ω$tmp*orbit["a"]][viewingAngle,0];h1lm[l,-m]Sslm[l,-m]Exp[I ω$tmp u$Signal],
{l,2,l$max},
{m,Max[l-4,1],l}];

Inspect GW waveform strain (normalized by the distance to the source)

ListPlot[{Re[GW$strain],-Im[GW$strain]},
PlotTheme->"Scientific",
FrameLabel->{"u/M","(strain) x (distance) / μ"},
ImageSize->Medium,
Joined->True]

Calculating adiabatic (0PA) inspirals

a0$val=0.7;
rISCO$val=KerrGeoSeparatrix[a0$val,0,1];
x0$vals=Range[(15^(3/2)+a0$val)^(-2/3),(rISCO$val^(3/2)+a0$val)^(-2/3),0.03];
r0$vals=(x0$vals^(-3/2)-a0$val)^(2/3);

Clear[circ$fluxes,circ$amp,circ$freq];
Monitor[Do[
orbit=KerrGeoOrbit[a0$val,r0,0,1];
Do[
teuk$Solutions[l,m]=TeukolskyPointParticleMode[-2,l,m,0,0,orbit];
Zinf= teuk$Solutions[l,m]["Amplitudes"]["ℐ"];
circ$freq[r0]=teuk$Solutions[l,m]["ω"];
circ$amp[r0][l,m]=-2 Zinf/circ$freq[r0]^2,
{l,2,6},{m,Max[l-2,1],l}];
circ$fluxes[r0]=Sum[2teuk$Solutions[l,m]["Fluxes"]["Energy"],{l,2,6},{m,Max[l-2,1],l}];,{r0,r0$vals}];,r0]

log$fluxes$data=Table[{r0,Log@Total[circ$fluxes[r0]]},{r0,r0$vals}];
Do[amps$data[l,m]=Table[{r0,circ$amp[r0][l,m]},{r0,r0$vals}],
{l,2,6},{m,Max[l-2,1],l}];

Do[amp$interp[l,m]=Interpolation[amps$data[l,m]],{l,2,6},{m,Max[l-2,1],l}];
LogEdot$interp=Interpolation[fluxes$data];

FullSimplify[D[KerrGeoEnergy[a0$val,r0,0,1],r0],{Im[r0]==0,r0>=6}]

r0dot[r0_]:=If[r0>=Min[r0$vals],-Exp[LogEdot$interp[r0]]/(((-0.735`+2.7999999999999994` (Sqrt[r0])-3.` r0+3.3306690738754696`*^-16 (r0^(3/2))+0.5` (r0^2))/((1.4` +(-3.`+r0) (Sqrt[r0])) r0 (Sqrt[1.4` (r0^(3/2))+(-3+r0) (r0^2)])))),0];

traj=Quiet@NDSolve[{r0'[t]==r0dot[r0[t]],r0[0]==10,Φ'[t]==If[r0dot[r0[t]]!=0,1/(r0[t]^(3/2)+a0$val),0],Φ[0]==0},{r0,Φ},{t,0,220}][[1]];
r0$traj=r0/.traj;
Φ$traj=Φ/.traj;
Plot[{r0$traj[t],Φ$traj[t]},{t,0,220},PlotTheme->"Scientific",FrameLabel->{"ϵ t/M"},PlotLegends->{"\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["r", "TI"], "0"], TraditionalForm], "errors" -> {}, "input" -> "r_0", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)","\!\(\*TemplateBox[<|"boxes" -> FormBox["ϵΦ", TraditionalForm], "errors" -> {}, "input" -> "\\epsilon \\Phi", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"}]