In [1]:
# Always run this first
# NOTE: Do not define new basic variables in this notebook;
# define them in Variables_Q.ipynb. Use this notebook
# to define new expressions built from those variables.
from __future__ import division # This needs to be here, even though it's in Variables_Q.ipynb
import sys
sys.path.insert(0, '..') # Look for modules in directory above this one
execfile('../Utilities/ExecNotebook.ipy')
try: execnotebook(VariablesNotebook)
except: execnotebook('Variables_Q.ipynb')
The following PNCollection objects will contain all the terms in the different parts of the binding energy.
In [2]:
AngularMomentum_NoSpin = PNCollection()
AngularMomentum_Spin = PNCollection()
In this notebook, every term will be multiplied by the following coefficient.
In [3]:
AngularMomentum_NoSpin.AddDerivedVariable('L_coeff', M**2*nu/v)
Note that fractions need to be entered as, e.g., frac(3,4) so that they are not converted to finite-precision decimals.
The nonspinning orbital binding energy is known through 4pN. The 5pN term has a known relationship to the 5pN binding energy term $e_5(\nu)$, though the latter is still incomplete. We set it to zero here. These expressions come from Eq. (230) and related footnotes of Blanchet (2006). Note that his calculation is for nonspinning systems, so he writes the quantity as $J$, which is assumed to be the component along $\hat{\ell}$.
In [4]:
e_5 = 0 # Placeholder for unknown term in energy expression.
In [5]:
AngularMomentum_NoSpin.AddDerivedConstant('L_0', ellHat,
datatype=ellHat.datatype)
# L_1 is 0
AngularMomentum_NoSpin.AddDerivedConstant('L_2', (frac(3,2) + frac(1,6)*nu)*ellHat,
datatype=ellHat.datatype)
# L_3 is 0
AngularMomentum_NoSpin.AddDerivedConstant('L_4', (frac(27,8) - frac(19,8)*nu + frac(1,24)*nu**2)*ellHat,
datatype=ellHat.datatype)
# L_5 is 0
AngularMomentum_NoSpin.AddDerivedConstant('L_6', (frac(135,16) + (-frac(6889,144) + frac(41,24)*pi**2)*nu
+ frac(31,24)*nu**2 + frac(7,1296)*nu**3)*ellHat,
datatype=ellHat.datatype)
# L_7 is 0
AngularMomentum_NoSpin.AddDerivedConstant('L_8',
(frac(2835,128) - frac(5,7)*nu*((-frac(123671,5760)+frac(9037,1536)*pi**2+frac(1792,15)*ln(2)+frac(896,15)*EulerGamma)
+ (-frac(498449,3456) + frac(3157,576)*pi**2)*nu
+ frac(301,1728)*nu**2 + frac(77,31104)*nu**3) + frac(64,35)*nu)*ellHat,
datatype=ellHat.datatype)
AngularMomentum_NoSpin.AddDerivedConstant('L_lnv_8', (-frac(128,3)*nu)*ellHat,
datatype=ellHat.datatype)
# L_9 is 0
# Below are the incomplete terms
AngularMomentum_NoSpin.AddDerivedConstant('L_10',
(frac(15309,256) + nu*(-frac(2,3)*e_5 - frac(4988,945) - frac(656,135)*nu))*ellHat,
datatype=ellHat.datatype)
AngularMomentum_NoSpin.AddDerivedConstant('L_lnv_10', ((frac(9976,105) + frac(1312,15)*nu)*nu*2)*ellHat,
datatype=ellHat.datatype)
(Look for spin-squared terms.)
The spin-orbit terms in the angular momentum are complete to 3.5pN. These terms come from Eq. (4.7) of Bohé et al. (2012):
In [6]:
# Lower-order terms are 0
AngularMomentum_Spin.AddDerivedVariable('L_SO_3',
(-(35*S_ell + 15*Sigma_ell*delta)/(6*M**2))*ellHat
+ ((S_n + Sigma_n*delta)/(2*M**2))*nHat
+ (-(3*S_lambda + Sigma_lambda*delta)/M**2)*lambdaHat,
datatype=ellHat.datatype)
# L_SO_4 is 0
AngularMomentum_Spin.AddDerivedVariable('L_SO_5',
(7*(61*S_ell*nu - 99*S_ell + 30*Sigma_ell*delta*nu - 27*Sigma_ell*delta)/(72*M**2))*ellHat
+ ((-19*S_n*nu + 33*S_n - 10*Sigma_n*delta*nu + 33*Sigma_n*delta)/(24*M**2))*nHat
+ ((18*S_lambda*nu - 21*S_lambda + 8*Sigma_lambda*delta*nu - 3*Sigma_lambda*delta)/(6*M**2))*lambdaHat,
datatype=ellHat.datatype)
# L_SO_6 is 0
AngularMomentum_Spin.AddDerivedVariable('L_SO_7',
((-29*S_ell*nu**2 + 1101*S_ell*nu - 405*S_ell - 15*Sigma_ell*delta*nu**2 + 468*Sigma_ell*delta*nu
- 81*Sigma_ell*delta)/(16*M**2))*ellHat
+ ((11*S_n*nu**2 - 1331*S_n*nu + 183*S_n + 5*Sigma_n*delta*nu**2 - 734*Sigma_n*delta*nu
+ 183*Sigma_n*delta)/(48*M**2))*nHat
+ ((-32*S_lambda*nu**2 + 2*S_lambda*nu - 174*S_lambda - 16*Sigma_lambda*delta*nu**2 - 79*Sigma_lambda*delta*nu
- 12*Sigma_lambda*delta)/(24*M**2))*lambdaHat,
datatype=ellHat.datatype)
In [7]:
def AngularMomentumExpression(AngularMomentumTerms=[AngularMomentum_NoSpin, AngularMomentum_Spin], PNOrder=frac(7,2)):
# We have to play some tricks with the log terms so that `horner` works
def logterm(key,val):
if 'lnv' in val:
return logv
else:
return 1
return L_coeff*horner(sum([key*(v**n)*logterm(key,val)
for Terms in AngularMomentumTerms
for n in range(2*PNOrder+1)
for key,val in Terms.items()
if val.endswith('_{0}'.format(n))])).subs(logv, ln(v))
In [8]:
# display(AngularMomentumExpression(PNOrder=frac(8,2)))