Introduction to PNSF: Definitions and expansion variables

The PNSF package provides analytic approximations of BHPT quantities as PN expansions

<<PostNewtonianSelfForce`

Common expansion variables

1/p  (*p the semi-latus rectum, p == ((2 r_min r_max)/(M (r_min + r_max))) *)
v==Sqrt[1/p] (*Characteristic velocity*)
y==(M omegaPhi)^(2/3)
e (*Darwin eccentricity, e == ((r_max - r_min)/(r_min + r_max))*)
q==a/M (*Spin of the large black hole in the Kerr case*)

Use “PostNewtonianExpansion” with the file path to load the data structure (similar to an association) into a variable.  Some examples:

eInfEcc=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Flux","EInfinity"}];
eHorCirc=PostNewtonianExpansion[{"Schwarzschild","Circular","Flux","Horizon"}];
redshiftEcc=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Local","RedShift","!-ey"}];
spinPrecEcc=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Local","Spin"}];
tidalBCirc=PostNewtonianExpansion[{"Schwarzschild","Circular","Local","Tidal-lambdaB"}];
kerrRedshiftCirc=PostNewtonianExpansion[{"Kerr","Circular","Local","RedShift"}];

redshiftEcc=PostNewtonianData[Get["~/Library/Mathematica/Applications/PostNewtonianSelfForce/SeriesData/Schwarzschild/Eccentric/Local/RedShift.m"]];

Descriptions of the above can be easily obtained

{eInfEcc["Name"],eInfEcc["Description"]}
{eHorCirc["Name"],eHorCirc["Description"]}
{redshiftEcc["Name"],redshiftEcc["Description"]}
{spinPrecEcc["Name"],spinPrecEcc["Description"]}
{tidalBCirc["Name"],tidalBCirc["Description"]}
{kerrRedshiftCirc["Name"],kerrRedshiftCirc["Description"]}

Tutorial: Loading and manipulating eccentric orbit PNSF data for the energy flux at infinity

Ensure that the directory is in Mathematica’s file path

$Path

Load the PNSF package

<<PostNewtonianSelfForce`

Use “PostNewtonianExpansion” with the file path to load the data structure (similar to an association) into a variable

energyInf=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Flux","EInfinity"}]

The data structure comes with various pieces of information

energyInf["Name"]
energyInf["Authors"]
energyInf["Description"]
energyInf["References"]

Of course, it also comes with the full series, which is not normalized.  It is recommended to truncate this series at a desired order.

energyInf["Series"]

energyInf["Series"]+O[y]^(17/2)

Better still is to use the “PostNewtonianCoefficient[energy, order, power_of_Log]” functionality to retrieve individual PN terms.

Closed-form eccentricity functions given as such

petersMathews=PostNewtonianCoefficient[energyInf,5,0]
SeriesCoefficient[energyInf["Series"],5]

onePNTerm=PostNewtonianCoefficient[energyInf,6,0]
sixPNLogSq=PostNewtonianCoefficient[energyInf,11,2]

sixPNLogSq/.e->1/2

Infinite eccentricity series which have been calculated to finite order are shown to the highest known order.   The “ResummedSeriesData” functionality is used to store the singular behavior and series (see ResummedSeries.m)

sixP5LogTerm=PostNewtonianCoefficient[energyInf,23/2,1]

sixP5LogTerm["SeriesData"]
sixP5LogTerm["SingularFactor"]

These can be retained to lower order, with or without the singular factor included

sixP5LogTerm+O[e]^11
Series[sixP5LogTerm,{e,0,10}]

sixP5LogTerm["SeriesData"]+O[e]^11

Functions which lack a closed form but can be expanded to any arbitary order are loaded as “extra functions”

PostNewtonianCoefficient[energyInf,13/2,0]

Expanding each extra function is simply done through the function itself.  Alternatively, a more complicated manipulation can be done

oneP5PN[e,16]
ReleaseHold[PostNewtonianCoefficient[energyInf,13/2,0]/.{highestPower->16}]

oneP5PN[e,16]["SeriesData"]

Obtain numeric approximations

Normal[oneP5PN[e,30]]/.e->0.9999000000000000000000000
Normal[oneP5PN[e,30]["SeriesData"]]/.e->0.99990000000000000000000

Some subtlety arises in the 3PN non-log term: chi3PN is a special function with both closed and (convergent) infinite-series portions

PostNewtonianCoefficient[energyInf,8,0]

The series portion of chi3PN is saved using ResummedSeriesData.  Removing the closed portion recovers the ResummedSeries

chi3PN[e,8]

(chi3PN[e,14]+((3 (1+((85 (e^2))/6)+((5171 (e^4))/192)+((1751 (e^6))/192)+((297 (e^8))/1024)) Log[1-(e^2)])/(2 (((1-(e^2)))^(13/2)))))["SeriesData"]

Note that this series converges as the eccentricity goes to 1

Similar procedures can be followed for the angular momentum flux at infinity

Ensure the package is loaded

<<PostNewtonianSelfForce`

Load the corresponding flux using “PostNewtonianExpansion”

angMomInf=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Flux","LInfinity"}]

The same information applies

angMomInf["Name"]
angMomInf["Authors"]
angMomInf["Description"]
angMomInf["References"]

The (unnormalized) series has a form similar to its energy counterpart

angMomInf["Series"]

angMomInf["Series"]+O[y]^(7)

“PostNewtonianCoefficient[energy, order, power_of_Log]” also functions in the same manner

Closed-form eccentricity functions:

peters=PostNewtonianCoefficient[angMomInf,7/2,0]
SeriesCoefficient[angMomInf["Series"],7/2]

onePNTerm=PostNewtonianCoefficient[angMomInf,9/2,0]
sixPNLogSq=PostNewtonianCoefficient[angMomInf,19/2,2]

sixPNLogSq/.e->1/2

Infinite eccentricity series which have been calculated to finite order are shown to the highest known order.   The “ResummedSeriesData” functionality is used to store the singular behavior and series (see ResummedSeries.m)

sixP5LogTerm=PostNewtonianCoefficient[angMomInf,10,1]

sixP5LogTerm["SeriesData"]
sixP5LogTerm["SingularFactor"]

These can be retained to lower order, with or without the singular factor included

sixP5LogTerm+O[e]^11
Series[sixP5LogTerm,{e,0,10}]

sixP5LogTerm["SeriesData"]+O[e]^11

Functions which lack a closed form but can be expanded to any arbitary order are loaded as “extra functions”

PostNewtonianCoefficient[angMomInf,5,0]

As above, expanding each extra function is most easily done through the function itself.

oneP5PNTild[e,16]
ReleaseHold[PostNewtonianCoefficient[angMomInf,5,0]/.{highestPower->16}]

oneP5PNTild[e,16]["SeriesData"]

f((128 π)/5)+((550 π (e^2))/3)+((787 π (e^4))/6)+((10007 π (e^6))/1440)+((2321 π (e^8))/34560)-((237857 π (e^10))/13824000)+((182863 π (e^12))/165888000)+((4987211 π (e^14))/260112384000)-((47839147 π (e^16))/1387266048000)+((O[e])^17)

Obtain numeric approximations

Normal[oneP5PNTild[e,30]]/.e->0.9999000000000000000000000
Normal[oneP5PNTild[e,30]["SeriesData"]]/.e->0.99990000000000000000000

Like chi3PN, chiTild3PN is a special function with both closed and (convergent) infinite-series portions

PostNewtonianCoefficient[angMomInf,13/2,0]

The series portion of chi3PN is saved using ResummedSeriesData.  Removing the closed portion recovers the ResummedSeries

chiTild3PN[e,8]

(chiTild3PN[e,14]+((3 (1+((229 (e^2))/32)+((327 (e^4))/64)+((69 (e^6))/256)) Log[1-(e^2)])/(2 (((1-(e^2)))^5))))["SeriesData"]

Note that this series converges as the eccentricity goes to 1

Exercises

Basic questions:

What is the coefficient of e^8 at 1/p^4 relative order in the eccentric-orbit Schwarzschild redshift invariant?

redshiftEcc=PostNewtonianExpansion[{"Schwarzschild","Eccentric","Local","RedShift","!-ey"}];
redshiftEcc["Series"]

redshiftEcc=
PostNewtonianData[Get["~/Library/Mathematica/Applications/PostNewtonianSelfForce/SeriesData/Schwarzschild/Eccentric/Local/RedShift.m"]];
PostNewtonianCoefficient[redshiftEcc,5,0]

What is the approximate limit as e->1 of the convergent series in the 6PN chi-like function for the eccentric Schwarzschild energy flux

chi6PNLog[e,6]

(chi6PNLog[e,8]+((3 (1+((16579 (e^2))/384)+((459595 (e^4))/1536)+((847853 (e^6))/1536)+((3672745 (e^8))/12288)+((1997845 (e^10))/49152)+((41325 (e^12))/65536)) Log[1-(e^2)])/(2 (((1-(e^2)))^(19/2)))))["SeriesData"]

N[Normal[(chi6PNLog[e,#]+((3 (1+((16579 (e^2))/384)+((459595 (e^4))/1536)+((847853 (e^6))/1536)+((3672745 (e^8))/12288)+((1997845 (e^10))/49152)+((41325 (e^12))/65536)) Log[1-(e^2)])/(2 (((1-(e^2)))^(19/2)))))["SeriesData"]]/.{e->1},10]&/@Range[10,50,5]

What is the approximate value of the Schwarzschild circular-orbit energy flux for p=10 to each PN order from 10 to 20

eInfCirc=PostNewtonianData[Get["~/Library/Mathematica/Applications/PostNewtonianSelfForce/SeriesData/Schwarzschild/Circular/Flux/Infinity.m"]]

(eInfCirc["Series"]+O[y]^(#))&/@Range[6,8]

N[(Normal[(eInfCirc["Series"]+O[y]^(#))]/.{y->1/10}),15]&/@Range[21/2,51/2]
%[[#+1]]-%[[#]]&/@Range[1,Length[%]-1]

N[(Normal[(eInfCirc["Series"]+O[y]^(#))]/.{y->Sqrt[1/10]}),10]&/@Range[25/2,51/2]
%[[#+1]]-%[[#]]&/@Range[1,Length[%]-1]

Using the energy and angular momentum fluxes to describe a Schwarzschild inspiral

<<PostNewtonianSelfForce`

<<KerrGeodesics`

Define E, L, y

Eecc=KerrGeoEnergy[0,p,e,1]
Lecc=KerrGeoAngularMomentum[0,p,e,1]
yPN=("\!\(\*SubscriptBox[\(Ω\), \(ϕ\)]\)"/.KerrGeoFrequencies[0,p,e,1])^(2/3)

Get dp/dt and de/dt in terms of dE/dt and dL/dt

Simplify[Solve[{Edot== D[Eecc,p]pdot +D[Eecc,e]edot,Ldot== D[Lecc,p]pdot +D[Lecc,e]edot},{pdot,edot}],Assumptions->{p>6+2e,0<e<1}]

Load the PN series, with the special functions extended to high order in e

EdotPN=PostNewtonianExpansion[{"Flux","Eccentric","EInfinity"}]["Series"]+O[y]^(8);
LdotPN=PostNewtonianExpansion[{"Flux","Eccentric","LInfinity"}]["Series"]+O[y]^(13/2);

EdotPN=Normal[ReleaseHold[EdotPN/.{highestPower->80}]]/.{y->yPN};
LdotPN=Normal[ReleaseHold[LdotPN/.{highestPower->80}]]/.{y->yPN};

pdot=((2 (EdotPN (Sqrt[((4 (e^2)-(((-2+p))^2))/(3+(e^2)-p))]) (3+(e^2)-p) (p^(3/2))+LdotPN (((-4+p))^2) (Sqrt[-3-(e^2)+p])))/(4 (e^2)-(((-6+p))^2)))/.{p->p[t],e->e[t]};
edot=((EdotPN (Sqrt[((4 (e^2)-(((-2+p))^2))/(3+(e^2)-p))]) (p^(3/2)) (18+2 (e^4)-3 (e^2) (-4+p)-9 p+(p^2))+(-1+(e^2)) LdotPN (Sqrt[-3-(e^2)+p]) (12+4 (e^2)-8 p+(p^2)))/(e (4 (e^2)-(((-6+p))^2)) p))/.{p->p[t],e->e[t]};

Solve for the evolution of p and e from some initial values

p0=20;
e0=0.5;
ϵ=1/1000;

inspiral=NDSolve[{p'[t]==ϵ pdot,e'[t]==ϵ edot,p[0]==p0,e[0]==e0},{p,e},{t,0,(10^7)}][[1]]

tMax=(p/.inspiral)["Domain"][[1,2]]

Plot[p[t]/.inspiral,{t,0,tMax}]

Plot[e[t]/.inspiral,{t,0,tMax}]