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]:
BindingEnergy_NoSpin = PNCollection()
BindingEnergy_Spin = PNCollection()
BindingEnergy_NSTides = PNCollection()
In this notebook, every term will be multiplied by the following coefficient.
In [3]:
BindingEnergy_NoSpin.AddDerivedVariable('E_coeff', -(M*nu*v**2)/2)
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 expressions through 3.5pN here come from Eq. (194) of Blanchet (2006).
The 4pN term from Eq. (5.2d) of Jaranowski and Schäfer is known exactly, now that the $\nu$-linear piece is given as Eq. (32) of Bini and Damour (2013a). The remaining terms are not known exactly, but Bini and Damour (2103b) have derived some terms, though there is incomplete information, which are noted as the constants in the following cell. Note that, though the notation is confusing, Bini and Damour claim they did not calculate the coefficient they call $a_6^{\text{ln 1}}$; but it seems to be given in their Eq. (64).
In [4]:
# Equation numbers below refer to v1 of Bini and Damour (2013b)
a_6__c1 = 0 # not yet known
a_6__ln1 = -frac(144,5) # coefficient of nu in Eq. (64)
a_65__c1 = 0 # not yet known
a_65__ln1 = 0 # not yet known
a_7__c1 = 0 # not yet known
a_7__ln1 = 0 # not yet known
In [5]:
BindingEnergy_NoSpin.AddDerivedConstant('E_0', 1)
# E_1 is 0
BindingEnergy_NoSpin.AddDerivedConstant('E_2', -frac(3,4) - frac(1,12)*nu)
# E_3 is 0
BindingEnergy_NoSpin.AddDerivedConstant('E_4', -frac(27,8) + frac(19,8)*nu - frac(1,24)*nu**2)
# E_5 is 0
BindingEnergy_NoSpin.AddDerivedConstant('E_6', -frac(675,64) + (frac(34445,576) - frac(205,96)*pi**2)*nu - frac(155,96)*nu**2
- frac(35,5184)*nu**3)
# E_7 is 0
BindingEnergy_NoSpin.AddDerivedConstant('E_8',
-frac(3969,128) + (-frac(123671,5760)+frac(9037,1536)*pi**2+frac(1792,15)*ln(2)+frac(896,15)*EulerGamma)*nu
+ (-frac(498449,3456) + frac(3157,576)*pi**2)*nu**2 + frac(301,1728)*nu**3 + frac(77,31104)*nu**4)
BindingEnergy_NoSpin.AddDerivedConstant('E_lnv_8', frac(896,15)*nu)
# E_9 is 0
# Below are the incomplete terms
BindingEnergy_NoSpin.AddDerivedConstant('E_10', -frac(45927,512)
+ (-frac(228916843,115200) - frac(9976,35)*EulerGamma + frac(729,7)*ln(3)
- frac(23672,35)*ln(2) + frac(126779,512)*pi**2)*nu
+ (-frac(21337,1024)*pi**2 + 3*a_6__c1 - frac(896,5)*ln(2)
- frac(448,5)*EulerGamma + frac(189745,576) + frac(2,3)*a_6__ln1)*nu**2
+ (-frac(1353,256)*pi**2 + frac(69423,512))*nu**3
+ frac(55,512)*nu**4
+ frac(1,512)*nu**5)
BindingEnergy_NoSpin.AddDerivedConstant('E_lnv_10', -frac(9976,35)*nu + (- frac(448,5) + 6*a_6__ln1)*nu**2)
BindingEnergy_NoSpin.AddDerivedConstant('E_11', frac(10,3)*nu * (frac(13696,525)*pi + nu*a_65__c1))
BindingEnergy_NoSpin.AddDerivedConstant('E_12',
- frac(264627,1024)+frac(2717,6718464)*nu**6+frac(5159,248832)*nu**5
+ (frac(272855,124416)*pi**2-frac(20543435,373248))*nu**4
+ (frac(1232,27)*EulerGamma+frac(6634243,110592)*pi**2-frac(11,2)*a_6__c1
-frac(75018547,51840)+frac(2464,27)*ln(2) -frac(20,9)*a_6__ln1)*nu**3
+ (frac(113594718743,14515200)+frac(18491,2304)*pi**4
+frac(246004,105)*ln(2)+frac(112772,105)*EulerGamma+frac(11,2)*a_6__c1+a_6__ln1+frac(2,3)*a_7__ln1
+ frac(11,3)*a_7__c1-frac(86017789,110592)*pi**2-frac(2673,14)*ln(3))*nu**2
+ (- frac(389727504721,43545600)+frac(74888,243)*ln(2) - frac(7128,7)*ln(3)-frac(30809603,786432)*pi**4
-frac(3934568,8505)*EulerGamma +frac(9118627045,5308416)*pi**2)*nu )
BindingEnergy_NoSpin.AddDerivedConstant('E_lnv_12',
frac(22,3)*a_7__ln1 - 2*frac(1967284,8505)*nu + 2*(frac(56386,105)+frac(11,2)*a_6__ln1)*nu**2
+ 2*(frac(616,27)-frac(11,2)*a_6__ln1)*nu**3)
(Is the following true? What about this paper?) The spin-squared terms (by which I mean both spin-spin and spin-orbit squared terms) in the energy are known only at 2pN order (from Kidder (1995) and Will and Wiseman (1996)). They are most conveniently given in Eq. (C4) of Arun et al.
In [6]:
# Lower-order terms are 0
BindingEnergy_Spin.AddDerivedVariable('E_SQ_4',
(1+delta-2*nu)*(chi1chi1+chi2chi2)/4 - 3*(chi_a_ell**2+chi_s_ell**2)/2
- delta*( chi2chi2/2 + 3*chi_a_ell*chi_s_ell ) + nu*( chi1chi2 + 6*chi_a_ell**2 ))
The spin-orbit terms in the energy are now complete to 4.0pN (the last term is zero). These terms come from Eq. (4.6) of Bohé et al. (2012):
In [7]:
# Lower-order terms are 0
BindingEnergy_Spin.AddDerivedVariable('E_SO_3', (frac(14,3)*S_ell + 2*delta*Sigma_ell)/M**2)
# E_SO_4 is 0
BindingEnergy_Spin.AddDerivedVariable('E_SO_5', ((11-61*nu/9)*S_ell + (3-10*nu/3)*delta*Sigma_ell)/M**2)
# E_SO_6 is 0
BindingEnergy_Spin.AddDerivedVariable('E_SO_7',
((frac(135,4)-frac(367,4)*nu+frac(29,12)*nu**2)*S_ell + (frac(27,4)-39*nu+frac(5,4)*nu**2)*delta*Sigma_ell)/M**2)
# E_SO_8 is 0
The tidal-coupling terms come in to the energy at relative 5pN order, and are known to 6pN order.
These terms come from Eq. (2.11) of Vines et al. (2011). Note their unusual convention for mass ratios, where $\chi_1 = m_1/m$ in their notation; in particular, $\chi$ is not a spin parameter. Also note that $\hat{\lambda} = \lambda_2 v^{10}/(m_1+m_2)^5$, and we need to add the coupling terms again with $1 \leftrightarrow 2$. Finally, note the normalization difference, where the overall factor is different by $-2$.
In [8]:
BindingEnergy_NSTides = PNCollection()
# Lower-order terms are 0
BindingEnergy_NSTides.AddDerivedConstant('E_NSTides_10', (-9*(M1/M2)*lambda2 - 9*(M2/M1)*lambda1)/M**5)
# E_NSTidal_11 is 0
BindingEnergy_NSTides.AddDerivedConstant('E_NSTides_12',
(-frac(11,2)*(M1/M2)*(3+2*M2/M+3*(M2/M)**2)*lambda2 - frac(11,2)*(M2/M1)*(3+2*M1/M+3*(M1/M)**2)*lambda1)/M**5)
In [9]:
def BindingEnergyExpression(BindingEnergyTerms=[BindingEnergy_NoSpin, BindingEnergy_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 E_coeff*horner(sum([key*(v**n)*logterm(key,val)
for Terms in BindingEnergyTerms
for n in range(2*PNOrder+1)
for key,val in Terms.items()
if val.endswith('_{0}'.format(n))])).subs(logv, ln(v))
def BindingEnergyDerivativeExpression(BindingEnergyTerms=[BindingEnergy_NoSpin, BindingEnergy_Spin], PNOrder=frac(7,2)):
Energy = BindingEnergyExpression(BindingEnergyTerms, PNOrder)
return horner(diff(Energy.subs(E_coeff, E_coeff.substitution), v)
.simplify().subs(log(v),logv)).simplify().subs(logv, ln(v))
In [10]:
# display(BindingEnergyExpression(PNOrder=frac(8,2)))
# display(BindingEnergyDerivativeExpression(PNOrder=frac(8,2)))