The status of various terms in the PN expression for $\frac{dE_{\rm GW}}{dt}$ is discussed in Ajith et al (2007). Ossokine et al (2015) is the primary reference for this notebook, and their expression was corrected and extended using Blanchet (2014) and Ajith et al (2007), respectively.
Notebook Status: Validated
Validation Notes: All expressions in this notebook were transcribed twice by hand on separate occasions, and expressions were corrected as needed to ensure consistency with published work. Published work was cross-validated and typo(s) in published work were corrected. In addition, this tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. Additional validation tests may have been performed, but are as yet, undocumented.
# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexpNRPyPN as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
from NRPyPN_shortcuts import div,dot,gamma_EulerMascheroni # NRPyPN: shortcuts for e.g., vector operations
This notebook uses three sources for the terms in $\frac{dE_{\rm GW}}{dt}$:
Constants appearing within $\frac{dE_{\rm GW}}{dt}$ expression in the above sources are as follows. Note that to reduce the possibility of copying error, these equations are taken directly from Eqs A1-A13 of the LaTeX source code of Ossokine et al (2015), and only mildly formatted to (1) improve presentation in Jupyter notebooks and (2) to ensure some degree of consistency in notation across different NRPyPN notebooks:
\begin{equation} \begin{split} m &= m_{1}+m_{2},\\ \nu &= \frac{m_{1}m_{2}}{m^{2}},\\ \delta &= \frac{m_{1}-m_{2}}{m},\\ \mathbf{S} &= \mathbf{S}_{1}+\mathbf{S}_{2},\\ s_{l} &= \frac{\mathbf{S}\cdot\hat\ell}{m^{2}},\\ s_{n} &= \frac{\mathbf{S}\cdot\hat n}{m^{2}},\\ \mathbf{\Sigma} &= \frac{m}{m_2}\mathbf{S}_{2} - \frac{m}{m_1}\mathbf{S}_{1}, \end{split} \quad\quad\quad \begin{split} \sigma_{l} &=\frac{\mathbf\Sigma\cdot\hat\ell}{m^{2}}, \\ \sigma_{n} &= \frac{\mathbf\Sigma\cdot\hat n}{m^{2}}, \\ \mathbf{\chi}_{s} &= \frac{1}{2}\left(\mathbf{\chi}_{1}+\mathbf{\chi}_{2}\right),\\ \mathbf{\chi}_{a} &= \frac{1}{2}\left(\mathbf{\chi}_{1}-\mathbf{\chi}_{2}\right),\\ \mathbf{S}_0 & = \frac{m}{m_1}\mathbf{S}_{1}+\frac{m}{m_2}\mathbf{S}_{2},\\ \mathbf{s}_0 & = \frac{\mathbf{S}_0}{m^2}, \end{split} \end{equation}where $\mathbf{\ell}$ is the unit normal to the instantaneous orbital plane. Since the instantaneous orbital plane here is assumed to be the $xy$ plane, below we set $\mathbf{\ell}=\{0,0,1 \}$.
In addition the Euler-Mascheroni constant $\gamma_{E}$ is used in the expression for $\frac{dE_{\rm GW}}{dt}$.
# Constants given in Eqs A1-13 of https://arxiv.org/abs/1502.01747
def dE_GW_dt_OBKPSS2015_consts(m1,m2, n12U, S1U,S2U):
# define scalars:
m = (m1+m2)
nu = m1*m2/m**2
delta = (m1-m2)/m
# define vectors:
Stot = ixp.zerorank1()
Sigma= ixp.zerorank1()
l = ixp.zerorank1()
l[2] = sp.sympify(1)
chi1U = ixp.zerorank1()
chi2U = ixp.zerorank1()
chi_s = ixp.zerorank1()
chi_a = ixp.zerorank1()
for i in range(3):
Stot[i] = S1U[i] + S2U[i]
Sigma[i] = (m1+m2)/m2*S2U[i] - (m1+m2)/m1*S1U[i]
chi1U[i] = S1U[i]/m1**2
chi2U[i] = S2U[i]/m2**2
chi_s[i] = div(1,2) * (chi1U[i] + chi2U[i])
chi_a[i] = div(1,2) * (chi1U[i] - chi2U[i])
# define scalars that depend on vectors
s_l = dot(Stot,l) /m**2
s_n = dot(Stot,n12U)/m**2
sigma_l = dot(Sigma,l)/m**2
sigma_n = dot(Sigma,n12U)/m**2
return nu,delta, l,chi_a,chi_s, s_l,s_n,sigma_l,sigma_n
As described in Ossokine et al (2015), the gravitational wave flux is given by
$$ \frac{dE}{dt} = -\mathcal{F}, $$where
$$ \mathcal{F}(x) = \frac{32}{5}\nu^{2}x^{5}\left(1+\sum_{k=2} b_{k}x^{k/2}\right), $$and (taking the following equations directly from Eqs A22-A28 of the LaTeX source code of Ossokine et al (2015), and only mildly formatted to (1) improve presentation in Jupyter notebooks and (2) to ensure some degree of consistency in notation across different NRPyPN notebooks):
\begin{eqnarray} b_{2} &=& -\frac{1247}{336}-\frac{35}{12}\nu, \\ b_{3} &=& 4\pi - 4s_l - \frac{5}{4}\delta\sigma_{l}, \\ b_{4} &=& -\frac{44711}{9072} + \frac{9271}{504}\nu + \frac{65}{18}\nu^{2} + \left(\frac{287}{96}+\frac{\nu}{24}\right)(\mathbf{\chi}_{s}\cdot\mathbf{\ell})^{2}\nonumber\\ &&- \left(\frac{89}{96}+\frac{7\nu}{24}\right)\mathbf{\chi}_{s}^{2} + \left(\frac{287}{96}-12\nu\right)(\mathbf{\chi}_{a}\cdot\mathbf{\ell})^{2} + \left(-\frac{89}{96}+4\nu\right)\mathbf{\chi}_{a}^{2} + \frac{287}{48}\delta(\mathbf{\chi}_{s}\cdot\mathbf{\ell})(\mathbf{\chi}_{a}\cdot\mathbf{\ell}) - \frac{89}{48}\delta(\mathbf{\chi}_{s}\cdot\mathbf{\chi}_{a}), \\ b_{5} &=& -\frac{8191}{672}\pi - \frac{9}{2}s_{l} - \frac{13}{16}\delta\sigma_{l} + \nu\left[-\frac{583}{24}\pi + \frac{272}{9}s_{l} + \frac{43}{4}\delta\sigma_{l}\right],\\ b_{6} &=& \frac{6643739519}{69854400} + \frac{16}{3}\pi^{2} - \frac{1712}{105}\gamma_{E} - \frac{856}{105}\log(16x)+ \left(\frac{-134543}{7776}+\frac{41}{48}\pi^{2}\right)\nu - \frac{94403}{3024}\nu^{2}-\frac{775}{324}\nu^{3} \nonumber\\ && - 16\pi s_{l} - \frac{31\pi}{6}\delta\sigma_{l},\\ b_{7}&=& \left(\frac{476645}{6804} + \frac{6172}{189}\nu - \frac{2810}{27}\nu^{2}\right)s_{l} + \left(\frac{9535}{336}+\frac{1849}{126}\nu - \frac{1501}{36}\nu^{2}\right)\delta\sigma_{l} \nonumber\\ &&+ \left(-\frac{16285}{504} \boxed{+} \frac{214745}{1728}\nu +\frac{193385}{3024}\nu^{2}\right)\pi,\\ b_{8} &=& \left(-\frac{3485\pi}{96} + \frac{13879\pi}{72}\nu\right)s_{l} + \left(-\frac{7163\pi}{672} + \frac{130583\pi}{2016}\nu\right)\delta\sigma_l; \end{eqnarray}notice the boxed term above, which contains a missing plus sign, as confirmed by comparing nonspinning terms to those found in Blanchet (2014). Finally,
$$ x = (m\Omega)^{2/3}. $$Finally, from Eq A.14 of Ajith et al (2007), we may include the time-changing mass of the black holes $\frac{dM}{dt}$ due to injected energy from tidal heating Alvi (2001):
$$ b_{5,\rm Mdot} = \left\{ -\frac{1}{4}\, \Big[ (1 - 3\nu) \chi_{s} (1 + 3\chi_{s}^{2} + 9\chi_{a}^{2}) + (1 - \nu) \delta \chi_{a} (1 + 3\chi_{a}^{2} + 9\chi_{s}^{2}) \Big] \right\}~, $$where $\chi_{s}=\mathbf{\chi_s}\cdot\mathbf{\ell}$ and $\chi_{a}=\mathbf{\chi_a}\cdot\mathbf{\ell}$ (i.e., the components aligned with the orbital angular momentum).
# Based on Eqs A22-28 of https://arxiv.org/abs/1502.01747, with
# Eq A.14 of https://arxiv.org/abs/0709.0093 for Mdot
# and correction on b[7] term by comparison with
# https://link.springer.com/content/pdf/10.12942/lrr-2014-2.pdf
def f_dE_GW_dt_and_dM_dt(mOmega, m1,m2, n12U, S1U,S2U):
def f_compute_quantities(mOmega, m1,m2, n12U, S1U,S2U, which_quantity):
if not (which_quantity == "dM_dt" or
which_quantity == "dE_GW_dt" or
which_quantity == "dE_GW_dt_plus_dM_dt"):
print("which_quantity == "+str(which_quantity)+" not supported!")
sys.exit(1)
nu,delta, l,chi_a,chi_s, s_l,s_n,sigma_l,sigma_n = dE_GW_dt_OBKPSS2015_consts(m1,m2, n12U, S1U,S2U)
x = (mOmega)**div(2,3)
# Compute b_5_Mdot:
b_5_Mdot = (-div(1,4)*(+(1-3*nu)*dot(chi_s,l)*(1+3*dot(chi_s,l)**2+9*dot(chi_a,l)**2)
+(1- nu)*delta*dot(chi_a,l)*(1+3*dot(chi_a,l)**2+9*dot(chi_s,l)**2)))
if which_quantity == "dM_dt":
return div(32,5)*nu**2*x**5*b_5_Mdot*x**div(5,2)
b = ixp.zerorank1(DIM=10)
b[2] = -div(1247,336) - div(35,12)*nu
b[3] = +4*sp.pi - 4*s_l - div(5,4)*delta*sigma_l
b[4] =(-div(44711,9072) + div(9271,504)*nu + div(65,18)*nu**2
+(+div(287,96) + div( 1,24)*nu)*dot(chi_s,l)**2
-(+div( 89,96) + div( 7,24)*nu)*dot(chi_s,chi_s)
+(+div(287,96) - 12*nu)*dot(chi_a,l)**2
+(-div( 89,96) + 4*nu)*dot(chi_a,chi_a)
+div(287,48)*delta*dot(chi_s,l)*dot(chi_a,l) - div(89,48)*delta*dot(chi_s,chi_a))
b[5] =(-div(8191,672)*sp.pi - div(9,2)*s_l - div(13,16)*delta*sigma_l
+nu*(-div(583,24)*sp.pi + div(272,9)*s_l + div(43,4)*delta*sigma_l))
if which_quantity == "dE_GW_dt_plus_dM_dt":
b[5]+= b_5_Mdot
b[6] =(+div(6643739519,69854400) + div(16,3)*sp.pi**2 - div(1712,105)*gamma_EulerMascheroni
-div(856,105)*sp.log(16*x) + (-div(134543,7776) + div(41,48)*sp.pi**2)*nu
-div(94403,3024)*nu**2 - div(775,324)*nu**3 - 16*sp.pi*s_l - div(31,6)*sp.pi*delta*sigma_l)
b[7] =(+(+div(476645,6804) + div(6172,189)*nu - div(2810,27)*nu**2)*s_l
+(+div(9535,336) + div(1849,126)*nu - div(1501,36)*nu**2)*delta*sigma_l
+(-div(16285,504) + div(214745,1728)*nu + div(193385,3024)*nu**2)*sp.pi)
b[8] =(+(-div(3485,96)*sp.pi + div(13879,72)*sp.pi*nu)*s_l
+(-div(7163,672)*sp.pi + div(130583,2016)*sp.pi*nu)*delta*sigma_l)
b_sum = sp.sympify(1)
for k in range(9):
b_sum += b[k]*x**div(k,2)
return div(32,5)*nu**2*x**5*b_sum
global dE_GW_dt_plus_dM_dt, dE_GW_dt, dM_dt
dE_GW_dt_plus_dM_dt = \
f_compute_quantities(mOmega, m1,m2, n12U, S1U,S2U, "dE_GW_dt_plus_dM_dt")
dE_GW_dt = f_compute_quantities(mOmega, m1,m2, n12U, S1U,S2U, "dE_GW_dt")
dM_dt = f_compute_quantities(mOmega, m1,m2, n12U, S1U,S2U, "dM_dt")
# Second version, used for validation purposes only.
# Eq A15 of https://arxiv.org/pdf/1502.01747.pdf
# Eq A.14 of https://arxiv.org/pdf/0709.0093.pdf for Mdot
def f_dE_GW_dt_plus_dM_dt_v2(Omega, m1,m2, n12U, S1U,S2U):
m = (m1+m2)
nu = m1*m2/m**2
delta = (m1-m2)/m
Stot = ixp.zerorank1()
Sigma= ixp.zerorank1()
l = ixp.zerorank1()
l[2] = sp.sympify(1)
chi1U = ixp.zerorank1()
chi2U = ixp.zerorank1()
chi_s = ixp.zerorank1()
chi_a = ixp.zerorank1()
for i in range(3):
Stot[i] = S1U[i] + S2U[i]
Sigma[i] = (m1+m2)/m2*S2U[i] - (m1+m2)/m1*S1U[i]
chi1U[i] = S1U[i]/m1**2
chi2U[i] = S2U[i]/m2**2
chi_s[i] = div(1,2) * (chi1U[i] + chi2U[i])
chi_a[i] = div(1,2) * (chi1U[i] - chi2U[i])
s_l = dot(Stot,l) /m**2
s_n = dot(Stot,n12U)/m**2
b_2 = -div(1247,336)-div(35,12)*nu
sigma_l = dot(Sigma,l)/m**2
sigma_n = dot(Sigma,n12U)/m**2
# We already multiply Omega by m=m1+m2
x = (Omega)**div(2,3)
b_0 = sp.sympify(1)
b_2 = -div(1247,336) - div(35,12)*nu
b_3 = 4*sp.pi - 4*s_l - div(5,4)*delta*sigma_l
b_4 = (-div(44711,9072) + div(9271,504)*nu + div(65,18)*nu**2
+ (div(287,96) + nu/24)*dot(chi_s,l)**2
-(div(89,96)+7*nu/24)*dot(chi_s,chi_s) + (div(287,96)-12*nu)*dot(chi_a,l)**2
+(-div(89,96)+4*nu)*dot(chi_a,chi_a) + div(287,48)*delta*dot(chi_s,l)*dot(chi_a,l)
-div(89,48)*delta*dot(chi_s,chi_a))
b_5 = (-div(8191,672)*sp.pi - div(9,2)*s_l - div(13,16)*delta*sigma_l
+nu*(-div(583,24)*sp.pi + div(272,9)*s_l + div(43,4)*delta*sigma_l))
b_6 = (+div(6643739519,69854400) + div(16,3)*sp.pi**2 - div(1712,105)*gamma_EulerMascheroni
-div(856,105)*sp.log(16*x) + (-div(134543,7776)+div(41,48)*sp.pi**2)*nu
-div(94403,3024)*nu**2 - div(775,324)*nu**3 - 16*sp.pi*s_l - div(31,6)*sp.pi*delta*sigma_l)
b_7 = (+(+div(476645,6804)+div(6172,189)*nu-div(2810,27)*nu**2)*s_l
+(+div(9535,336) + div(1849,126) *nu - div(1501,36) *nu**2)*delta*sigma_l
+(-div(16285,504) + div(214745,1728)*nu + div(193385,3024)*nu**2)*sp.pi)
b_8 = (+(-div(3485,96)*sp.pi + div(13879,72)*sp.pi*nu)*s_l
+(-div(7163,672)*sp.pi + div(130583,2016)*sp.pi*nu)*delta*sigma_l)
# Eq A.14 of https://arxiv.org/pdf/0709.0093.pdf for Mdot.
# note that this expression only considers spin-aligned cases,
# so we dot the chi's with the zhat direction.
b_5_Mdot = -div(1,4)*(+(1-3*nu)*dot(chi_s,l)*(1+3*dot(chi_s,l)**2+9*dot(chi_a,l)**2)
+(1- nu)*delta*dot(chi_a,l)*(1+3*dot(chi_a,l)**2+9*dot(chi_s,l)**2))
global dE_GW_dt_plus_dM_dtv2
dE_GW_dt_plus_dM_dtv2 = div(32,5)*nu**2*x**5*(+b_0
+b_2*x**div(2,2)
+b_3*x**div(3,2)
+b_4*x**div(4,2)
+(b_5+b_5_Mdot)*x**div(5,2)
+b_6*x**div(6,2)
+b_7*x**div(7,2)
+b_8*x**div(8,2))
from NRPyPN_shortcuts import m1,m2,n12U,S1U,S2U # NRPyPN: import needed input variables.
Omega = sp.symbols('Omega',real=True)
f_dE_GW_dt_and_dM_dt(Omega, m1,m2, n12U, S1U,S2U)
def error(varname):
print("ERROR: When comparing Python module & notebook, "+varname+" was found not to match.")
sys.exit(1)
# Validation against second transcription of the expressions:
f_dE_GW_dt_plus_dM_dt_v2(Omega, m1,m2, n12U, S1U,S2U)
f_dE_GW_dt_and_dM_dt(Omega, m1,m2, n12U, S1U,S2U)
if sp.simplify(dE_GW_dt_plus_dM_dt - dE_GW_dt_plus_dM_dtv2) != 0: error("dE_GW_dt_plus_dM_dtv2")
# Validation against corresponding Python module:
import PN_dE_GW_dt_and_dM_dt as dEdt
dEdt.f_dE_GW_dt_and_dM_dt(Omega, m1,m2, n12U, S1U,S2U)
if sp.simplify(dE_GW_dt_plus_dM_dt - dEdt.dE_GW_dt_plus_dM_dt) != 0: error("dE_GW_dt_plus_dM_dt")
if sp.simplify(dE_GW_dt - dEdt.dE_GW_dt) != 0: error("dE_GW_dt")
if sp.simplify(dM_dt - dEdt.dM_dt) != 0: error("dM_dt")
# Finally, confirm that dE_GW_dt_plus_dM_dt = dE_GW_dt + dM_dt
if sp.simplify(dE_GW_dt_plus_dM_dt - (dE_GW_dt + dM_dt)) != 0: error("dE_GW_dt + dM_dt")
print("ALL TESTS PASS")
ALL TESTS PASS
The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename PN-dE_GW_dt_and_dM_dt.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import os,sys # Standard Python modules for multiplatform OS-level functions
import cmdline_helperNRPyPN as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("PN-dE_GW_dt_and_dM_dt")
Created PN-dE_GW_dt_and_dM_dt.tex, and compiled LaTeX file to PDF file PN- dE_GW_dt_and_dM_dt.pdf