Notebook

$p_r$ , the radial component of the momentum vector, up to and including 3.5 PN order¶

Author: Zach Etienne¶

This notebook constructs the radial component of the momentum vector, up to and including 3.5 post-Newtonian order, and performs numerous validation tests.¶

Notebook Status: Validated

Validation Notes: All PN 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, all expressions in this notebook were validated against those in the Mathematica notebook used by Ramos-Buades, Husa, and Pratten (2018) (thanks to Toni Ramos-Buades for sharing this!) Finally, 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.ß

This notebook exists as the following Python module:¶

PN_p_r.py

This notebook & corresponding Python module depend on the following NRPy+/NRPyPN Python modules:¶

Table of Contents¶

$\label{toc}$

Part 0: Import needed Python modules
Part 1: Basic strategy for computing pr
1. Part 1.a: Construct the full Hamiltonian consistent with a binary orbiting instantaneously in the $x$ - $y$ plane
2. Part 1.b Computing $\frac{dr}{dt}$
3. Part 1.c Computing pr(drdt,...) using two approaches
  1. Part 1.c.i Approach 1: Use the full Hamiltonian
  2. Part 1.c.ii Approach 2 (default approach): use the full Hamiltonian, with higher-order terms over the course of derivation removed
Part 2: Validation Tests
1. Part 2.a: Validation of transcribed versions for $p_r$
2. Part 2.b: Comparison of two approaches for computing $p_r$
3. Part 2.c: Validation against trusted numerical values (i.e., in Table V of Ramos-Buades, Husa, and Pratten (2018))
4. Part 2.d: Validation against corresponding Python module
Part 2: LaTeX PDF output: $\LaTeX$ PDF Output

Part 0: Import needed Python modules [Back to top]¶

$\label{imports}$

In [1]:

# 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
from NRPyPN_shortcuts import Pt,Pr,nU,div # NRPyPN: shortcuts for e.g., vector operations

Part 1: Basic strategy for computing $p_r$ [Back to top]¶

$\label{strategy}$

Part 1.a: Construct the full Hamiltonian consistent with a binary orbiting instantaneously in the $x$ - $y$ plane [Back to top]¶

$\label{fullhamiltonian}$

As detailed in Ramos-Buades, Husa, and Pratten (2018), the basic strategy for computing $p_r$ first involves constructing the full Hamiltonian expression assuming that, consistent with a binary system initially on the $x$ -axis and orbiting instantaneously on the $x$ - $y$ plane, the momenta and normal vectors are given by:

$\begin{align} P_1^x &= -p_r\\ P_2^x &= +p_r\\ P_1^y &= +p_t\\ P_2^y &= -p_t\\ P_1^z = P_2^z &= 0\\ \mathbf{n} &= (1,0,0) \end{align}$

Now let's construct the full Hamiltonian, applying these assumptions as we go:

In [2]:

# Step 1: Construct full Hamiltonian
#         expression for a binary instantaneously
#         orbiting on the xy plane, store
#         result to Htot_xyplane_binary
def f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r):
    def make_replacements(expr):
        zero = sp.sympify(0)
        one = sp.sympify(1)
        return expr.subs(p1U[1], Pt).subs(p2U[1], -Pt).subs(p1U[2], zero).subs(p2U[2], zero).subs(p1U[0], -Pr).subs(
            p2U[0], Pr) \
            .subs(nU[0], one).subs(nU[1], zero).subs(nU[2], zero)

    import PN_Hamiltonian_NS as H_NS
    H_NS.f_H_Newt__H_NS_1PN__H_NS_2PN(m1, m2, p1U, n12U, r)
    H_NS.f_H_NS_3PN(m1, m2, p1U, n12U, r)

    global Htot_xyplane_binary
    Htot_xyplane_binary = make_replacements(+H_NS.H_Newt
                                            + H_NS.H_NS_1PN
                                            + H_NS.H_NS_2PN
                                            + H_NS.H_NS_3PN)

    import PN_Hamiltonian_SO as H_SO
    H_SO.f_H_SO_1p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    H_SO.f_H_SO_2p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    H_SO.f_H_SO_3p5PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    Htot_xyplane_binary += make_replacements(+H_SO.H_SO_1p5PN
                                             + H_SO.H_SO_2p5PN
                                             + H_SO.H_SO_3p5PN)

    import PN_Hamiltonian_SS as H_SS
    H_SS.f_H_SS_2PN(m1, m2, S1U, S2U, nU, r)
    H_SS.f_H_SS_S1S2_3PN(m1, m2, n12U, S1U, S2U, p1U, p2U, r)
    H_SS.f_H_SS_S1sq_S2sq_3PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    Htot_xyplane_binary += make_replacements(+H_SS.H_SS_2PN
                                             + H_SS.H_SS_S1S2_3PN
                                             + H_SS.H_SS_S1sq_S2sq_3PN)

    import PN_Hamiltonian_SSS as H_SSS
    H_SSS.f_H_SSS_3PN(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    Htot_xyplane_binary += make_replacements(+H_SSS.H_SSS_3PN)

Part 1.b: Computing $\frac{dr}{dt}$ [Back to top]¶

$\label{dr_dt}$

Hamilton's equations of motion imply that

$\frac{dr}{dt} = \frac{\partial H}{\partial p_r}.$

Next we Taylor expand $\frac{\partial H}{\partial p_r}$ in powers of $p_r$ , about $p_r=0$ :

$\begin{align} \frac{dr}{dt} = \frac{\partial H}{\partial p_r} = \left.\frac{\partial H}{\partial p_r}\right|_{p_r=0} + \frac{1}{1!} p_r \left.\frac{\partial^2 H}{\partial p_r^2}\right|_{p_r=0} + \frac{1}{2!} p_r^2 \left.\frac{\partial^3 H}{\partial p_r^3}\right|_{p_r=0} + \frac{1}{3!} p_r^3 \left.\frac{\partial^4 H}{\partial p_r^4}\right|_{p_r=0} + \mathcal{O}(p_r^4) \end{align}$

so to first order we get

$p_r \approx \left(\frac{dr}{dt} - \left.\frac{\partial H}{\partial p_r}\right|_{p_r=0} \right) \left( \left.\frac{\partial^2 H}{\partial p_r^2}\right|_{p_r=0} \right)^{-1}.$

Given the input masses and spins, we can compute $p_t$ using the formula given in this NRPyPN notebook (from Ramos-Buades, Husa, and Pratten (2018)), and the above will lead us with one equation and two unknowns: $p_r$ and $\frac{dr}{dt}$ . This is an equation we can derive directly from our Hamiltonian expression, and compare the output to the expression derived to 3.5PN order in Ramos-Buades, Husa, and Pratten (2018).

Let's next construct an expression for $\frac{dr}{dt}$ in terms of known quantities, via

$\frac{dr}{dt}=\left(\frac{dE_{\rm GW}}{dt}+\frac{dM}{dt}\right)\left[\frac{dH_{\rm circ}}{dr}\right]^{-1},$

where

$\frac{dH_{\rm circ}(r,p_t(r))}{dr} = \frac{\partial H(p_r=0)}{\partial r} + \frac{\partial H(p_r=0)}{\partial p_t} \frac{\partial p_t}{\partial r},$

and the total energy flux

$\mathcal{E}(M\Omega,...) = \left(\frac{dE_{\rm GW}}{dt}+\frac{dM}{dt}\right)$

is given in terms of input parameters (e.g., masses and spins), plus $M\Omega$ , which we also constructed in terms of the input masses, spins, and orbital separation $r$ .

In [3]:

# Function for computing dr/dt
def f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r):
    # First compute p_t
    import PN_p_t as pt
    pt.f_p_t(m1,m2, chi1U,chi2U, r)

    # Then compute dH_{circ}/dr = partial_H(p_r=0)/partial_r
    #                                  + partial_H(p_r=0)/partial_{p_t} partial_{p_t}/partial_r
    dHcirc_dr = (+sp.diff(Htot_xyplane_binary.subs(Pr,sp.sympify(0)),r)
                 +sp.diff(Htot_xyplane_binary.subs(Pr,sp.sympify(0)),Pt)*sp.diff(pt.p_t,r))

    # Then compute M\Omega
    import PN_MOmega as MOm
    MOm.f_MOmega(m1,m2, chi1U,chi2U, r)

    # Next compute dE_GW_dt_plus_dM_dt
    import PN_dE_GW_dt_and_dM_dt as dEdt
    dEdt.f_dE_GW_dt_and_dM_dt(MOm.MOmega, m1,m2, n12U, S1U,S2U)

    # Finally, compute dr/dt
    global dr_dt
    dr_dt = dEdt.dE_GW_dt_plus_dM_dt / dHcirc_dr

Part 1.c: Computing $p_r\left(\frac{dr}{dt},...\right)$ using two approaches [Back to top]¶

$\label{p_r}$

Part 1.c.i: Approach 1: Use the full Hamiltonian [Back to top]¶

$\label{approach1_fullham}$

First, we use the approximation based on the Hamiltonian computed above:

$p_r \approx \left(\frac{dr}{dt} - \left.\frac{\partial H}{\partial p_r}\right|_{p_r=0} \right) \left[ \left.\frac{\partial^2 H}{\partial p_r^2}\right|_{p_r=0} \right]^{-1}$

In [4]:

# Next we compute p_r as a function of dr_dt (unknown) and known quantities using
# p_r \approx [dr/dt - (partial_H/partial_{p_r})|_{p_r=0}] * [(partial^2_{H}/partial_{p_r^2})|_{p_r=0}]^{-1}
def f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):
    f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)

    dHdpr_przero   = sp.diff(Htot_xyplane_binary,Pr).subs(Pr,sp.sympify(0))
    d2Hdpr2_przero = sp.diff(sp.diff(Htot_xyplane_binary,Pr),Pr).subs(Pr,sp.sympify(0))
    global p_r_fullHam
    p_r_fullHam = (dr_dt - dHdpr_przero)/(d2Hdpr2_przero)

Part 1.c.ii: Approach 2 (default approach): use the full Hamiltonian, with higher-order terms over the course of derivation removed [Back to top]¶

$\label{approach2_default}$

This approach (Eq 2.16 of Ramos-Buades, Husa, and Pratten (2018)) uses the same approximation as Approach 1

$p_r \approx \left(\frac{dr}{dt} - \left.\frac{\partial H}{\partial p_r}\right|_{p_r=0} \right) \left[ \left.\frac{\partial^2 H}{\partial p_r^2}\right|_{p_r=0} \right]^{-1},$ except it throws away higher-order terms.

To reduce the possibility of copying error, the equation for $p_r$ is taken directly from the arXiv LaTeX source code of Eq 2.16 in Ramos-Buades, Husa, and Pratten (2018), and only mildly formatted to (1) improve presentation in Jupyter notebooks and (2) to ensure some degree of consistency in notation across different terms in other NRPyPN notebooks:

$\begin{equation} \begin{split} p_r &= \left[-\frac{dr}{dt} \right.\\ &\quad\quad \left. +\frac{1}{r^{7/2}} \left( -\frac{(6 q+13) q^2 \chi _{1x} \chi _{2 y}}{4 (q+1)^4}-\frac{(6 q+1) q^2 \chi _{2 x} \chi _{2 y}}{4 (q+1)^4}+\chi _{1y} \left(-\frac{q (q+6) \chi _{1x}}{4 (q+1)^4}-\frac{q (13 q+6) \chi _{2 x}}{4 (q+1)^4}\right)\right) \right.\\ &\quad\quad \left. +\frac{1}{r^4} \left( \chi _{1z} \left(\frac{3 q (5 q+2) \chi _{1x} \chi _{2 y}}{2 (q+1)^4} -\frac{3 q^2 (2 q+5) \chi _{2 x} \chi _{2 y}}{2 (q+1)^4}\right)+\chi _{1y} \chi _{2 z} \left(\frac{3 q^2 (2 q+5) \chi _{2 x}}{2 (q+1)^4}-\frac{3 q (5 q+2) \chi _{1x}}{2 (q+1)^4}\right)\right)\right] \\ & \times \left[ -\frac{(q+1)^2}{q}-\frac{1 \left(-7 q^2-15 q-7\right)}{2 q r} \right. \\ &\quad\quad \left. -\frac{47 q^4+229 q^3+363 q^2+229 q+47}{8 q (q+1)^2 r^2} -\frac{1}{r^{5/2}}\left( \frac{\left(4 q^2+11 q+12\right) \chi _{1z}}{4 q (q+1)}+\frac{\left(12 q^2+11 q+4\right) \chi _{2 z}}{4 (q+1)} \right) \right. \\ &\quad\quad \left. \left.- \frac{1}{r^{7/2}} \left( \frac{\left(-53 q^5-357 q^4-1097 q^3-1486 q^2-842 q-144\right) \chi _{1z}}{16 q (q+1)^4}+\frac{\left(-144 q^5-842 q^4-1486 q^3-1097 q^2-357 q-53\right) \chi _{2 z}}{16 (q+1)^4} \right)\right. \right. \\ &\quad\quad \left. -\frac{1}{r^3} \left(\frac{\left(q^2+9 q+9\right) \chi _{1x}^2}{2 q (q+1)^2}+\frac{\left(3 q^2+5 q+3\right) \chi _{2 x} \chi _{1x}}{(q+1)^2}+\frac{\left(3 q^2+8 q+3\right) \chi _{1y} \chi _{2 y}}{2 (q+1)^2}-\frac{9 q^2 \chi _{2 y}^2}{4 (q+1)}+\frac{\left(3 q^2+8 q+3\right) \chi _{1z} \chi _{2 z}}{2 (q+1)^2}-\frac{9 q^2 \chi _{2 z}^2}{4 (q+1)} \right. \right. \\ &\quad\quad\quad \left. \left. +\frac{\left(9 q^3+9 q^2+q\right) \chi _{2 x}^2}{2 (q+1)^2}+\frac{-363 q^6-2608 q^5-7324 q^4-10161 q^3-7324 q^2-2608 q-363}{48 q (q+1)^4}-\frac{9 \chi _{1y}^2}{4 q (q+1)}-\frac{9 \chi _{1z}^2}{4 q (q+1)}-\frac{\pi ^2}{16} \right) \right]^{-1}. \end{split} \end{equation}$

In [5]:

# Here's the Ramos-Buades, Husa, and Pratten (2018)
#  approach for computing p_r.
# Transcribed from Eq 2.18 of
# Ramos-Buades, Husa, and Pratten (2018),
#   https://arxiv.org/abs/1810.00036
def f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):
    q = m2/m1 # It is assumed that q >= 1, so m2 >= m1.
    f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)
    chi1x = chi1U[0]
    chi1y = chi1U[1]
    chi1z = chi1U[2]
    chi2x = chi2U[0]
    chi2y = chi2U[1]
    chi2z = chi2U[2]
    p_r_num = (-dr_dt
               +(-(6*q+13)*q**2*chi1x*chi2y/(4*(q+1)**4)
                 -(6*q+ 1)*q**2*chi2x*chi2y/(4*(q+1)**4)
                 +chi1y*(-q*(   q+6)*chi1x/(4*(q+1)**4)
                         -q*(13*q+6)*chi2x/(4*(q+1)**4)))/r**div(7,2)
               +(+chi1z*(+3*q   *(5*q+2)*chi1x*chi2y/(2*(q+1)**4)
                         -3*q**2*(2*q+5)*chi2x*chi2y/(2*(q+1)**4))
                 +chi1y*chi2z*(+3*q**2*(2*q+5)*chi2x/(2*(q+1)**4)
                               -3*q   *(5*q+2)*chi1x/(2*(q+1)**4)))/r**4)
    p_r_den = (-(q+1)**2/q - (-7*q**2-15*q-7)/(2*q*r)
               -(47*q**4 + 229*q**3 + 363*q**2 + 229*q + 47)/(8*q*(q+1)**2*r**2)
               -(+( 4*q**2 + 11*q + 12)*chi1z/(4*q*(q+1))
                 +(12*q**2 + 11*q +  4)*chi2z/(4*  (q+1)))/r**div(5,2)
               -(+(- 53*q**5 - 357*q**4 - 1097*q**3 - 1486*q**2 - 842*q - 144)*chi1z/(16*q*(q+1)**4)
                 +(-144*q**5 - 842*q**4 - 1486*q**3 - 1097*q**2 - 357*q -  53)*chi2z/(16  *(q+1)**4))/r**div(7,2)
               -(+(  q**2 + 9*q + 9)*chi1x**2/(2*q*(q+1)**2)
                 +(3*q**2 + 5*q + 3)*chi2x*chi1x/((q+1)**2)
                 +(3*q**2 + 8*q + 3)*chi1y*chi2y/(2*(q+1)**2)
                 -9*q**2*chi2y**2/(4*(q+1))
                 +(3*q**2 + 8*q + 3)*chi1z*chi2z/(2*(q+1)**2)
                 -9*q**2*chi2z**2/(4*(q+1))
                 +(9*q**3 + 9*q**2 + q)*chi2x**2/(2*(q+1)**2)
                 +(-363*q**6 - 2608*q**5 - 7324*q**4 - 10161*q**3 - 7324*q**2 - 2608*q - 363)/(48*q*(q+1)**4)
                 -9*chi1y**2/(4*q*(q+1))
                 -9*chi1z**2/(4*q*(q+1)) - sp.pi**2/16)/r**3)
    global p_r
    p_r = p_r_num/p_r_den

In [6]:

# Second version, used for validation purposes only.
def f_p_r_RHP2018v2(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):
    q = m2/m1
    f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)
    chi1x = chi1U[0]
    chi1y = chi1U[1]
    chi1z = chi1U[2]
    chi2x = chi2U[0]
    chi2y = chi2U[1]
    chi2z = chi2U[2]
    p_r_num = (-dr_dt
               +(-(6*q+13)*q**2*chi1x*chi2y/(4*(q+1)**4)
                 -(6*q+1)*q**2*chi2x*chi2y/(4*(q+1)**4)
                 +chi1y*(-q*(q+6)*chi1x/(4*(q+1)**4) - q*(13*q+6)*chi2x/(4*(q+1)**4)))/r**div(7,2)
               +(+chi1z*(+div(3,2)*(q*   (5*q+2)*chi1x*chi2y)/(q+1)**4
                         -div(3,2)*(q**2*(2*q+5)*chi2x*chi2y)/(q+1)**4)
                 +chi1y*chi2z*(+div(3,2)*(q**2*(2*q+5)*chi2x)/(q+1)**4
                               -div(3,2)*(q*   (5*q+2)*chi1x)/(q+1)**4))/r**4)
    p_r_den = (-(q+1)**2/q
               -(-7*q**2-15*q-7)/(2*q*r)
               -(47*q**4 + 229*q**3 + 363*q**2 + 229*q + 47)/(8*q*(q+1)**2*r**2)
               -( (4*q**2+11*q+12)*chi1z/(4*q*(q+1)) + (12*q**2+11*q+4)*chi2z/(4*(q+1)) )/r**div(5,2)
               -(+(-53 *q**5 - 357*q**4 - 1097*q**3 - 1486*q**2 - 842*q - 144)*chi1z/(16*q*(q+1)**4)
                 +(-144*q**5 - 842*q**4 - 1486*q**3 - 1097*q**2 - 357*q - 53 )*chi2z/(16*  (q+1)**4) )/r**div(7,2)
               -(+(  q**2+9*q+9)*chi1x**2/(2*q*(q+1)**2)
                 +(3*q**2+5*q+3)*chi2x*chi1x/((q+1)**2)
                 +(3*q**2+8*q+3)*chi1y*chi2y/(2*(q+1)**2)
                 -(9*q**2*chi2y**2/(4*(q+1)))
                 +(3*q**2+8*q+3)*chi1z*chi2z/(2*(q+1)**2)
                 -(9*q**2*chi2z**2/(4*(q+1)))
                 +(9*q**3+9*q**2+q)*chi2x**2/(2*(q+1)**2)
                 +(-363*q**6 - 2608*q**5 - 7324*q**4 - 10161*q**3 - 7324*q**2 - 2608*q - 363)/(48*q*(q+1)**4)
                 -(9*chi1y**2)/(4*q*(q+1))
                 -(9*chi1z**2)/(4*q*(q+1))
                 -sp.pi**2/16) / r**3)
    global p_r_RHP2018v2
    p_r_RHP2018v2 = p_r_num/p_r_den

In [7]:

# Third version, directly from Toni Ramos-Buades' Mathematica notebook (thanks Toni!)
def f_p_r_RHP2018v3(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r):
    q = m2/m1
    f_Htot_xyplane_binary(m1, m2, n12U, n21U, S1U, S2U, p1U, p2U, r)
    f_dr_dt(Htot_xyplane_binary, m1,m2, n12U, chi1U,chi2U, S1U,S2U, r)
    chi1x = chi1U[0]
    chi1y = chi1U[1]
    chi1z = chi1U[2]
    chi2x = chi2U[0]
    chi2y = chi2U[1]
    chi2z = chi2U[2]
    Pi = sp.pi
    p_r_num = (-dr_dt
               +((chi1y*chi2z*((3*chi2x*q**2*(5 + 2*q))/(2*(1 + q)**4) - (3*chi1x*q*(2 + 5*q))/(2*(1 + q)**4)) + chi1z*((-3*chi2x*chi2y*q**2*(5 + 2*q))/(2*(1 + q)**4) + (3*chi1x*chi2y*q*(2 + 5*q))/(2*(1 + q)**4)))/
        r**4 + (-(chi2x*chi2y*q**2*(1 + 6*q))/(4*(1 + q)**4) - (chi1x*chi2y*q**2*(13 + 6*q))/(4*(1 + q)**4) + chi1y*(-(chi1x*q*(6 + q))/(4*(1 + q)**4) - (chi2x*q*(6 + 13*q))/(4*(1 + q)**4)))/r**(sp.Rational(7,2))))
    p_r_den = (-((1 + q)**2/q) - ((chi2z*(-53 - 357*q - 1097*q**2 - 1486*q**3 - 842*q**4 - 144*q**5))/(16*(1 + q)**4) + (chi1z*(-144 - 842*q - 1486*q**2 - 1097*q**3 - 357*q**4 - 53*q**5))/(16*q*(1 + q)**4))/r**(sp.Rational(7,2)) -
       (-Pi**2/16 - (9*chi1y**2)/(4*q*(1 + q)) - (9*chi1z**2)/(4*q*(1 + q)) - (9*chi2y**2*q**2)/(4*(1 + q)) - (9*chi2z**2*q**2)/(4*(1 + q)) + (chi1x**2*(9 + 9*q + q**2))/(2*q*(1 + q)**2) +
          (chi1x*chi2x*(3 + 5*q + 3*q**2))/(1 + q)**2 + (chi1y*chi2y*(3 + 8*q + 3*q**2))/(2*(1 + q)**2) + (chi1z*chi2z*(3 + 8*q + 3*q**2))/(2*(1 + q)**2) + (chi2x**2*(q + 9*q**2 + 9*q**3))/(2*(1 + q)**2) +
          (-363 - 2608*q - 7324*q**2 - 10161*q**3 - 7324*q**4 - 2608*q**5 - 363*q**6)/(48*q*(1 + q)**4))/r**3 - ((chi1z*(12 + 11*q + 4*q**2))/(4*q*(1 + q)) + (chi2z*(4 + 11*q + 12*q**2))/(4*(1 + q)))/r**(sp.Rational(5,2)) -
       (47 + 229*q + 363*q**2 + 229*q**3 + 47*q**4)/(8*q*(1 + q)**2*r**2) - (-7 - 15*q - 7*q**2)/(2*q*r))
    global p_r_RHP2018v3
    p_r_RHP2018v3 = p_r_num/p_r_den

Part 2: Validation Tests [Back to top]¶

$\label{code_validation}$

Part 2.a: Validation of transcribed versions for $p_r$ [Back to top]¶

$\label{code_validation_trans}$

As a code validation check, we verify agreement between

the SymPy expressions transcribed from the cited published work on two separate occasions, and

In [8]:

from NRPyPN_shortcuts import m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r # NRPyPN: import needed input variables

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_p_r(          m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
f_p_r_RHP2018v2(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
f_p_r_RHP2018v3(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)

if sp.simplify(p_r - p_r_RHP2018v2)           != 0: error("p_r_RHP2018v2")
if sp.simplify(p_r_RHP2018v2 - p_r_RHP2018v3) != 0: error("p_r_RHP2018v3")
print("TRANSCRIPTION TEST PASSES")

TRANSCRIPTION TEST PASSES

Part 2.b: Comparison of two approaches for computing $p_r$ [Back to top]¶

$\label{code_validation_2approaches}$

As a validation check, we inject random inputs into the two approaches for computing $p_r$ outlined above and compute the errors.

In [9]:

import random # Standard Python module: provides random number generation functionality
from NRPyPN_shortcuts import num_eval, gamma_EulerMascheroni # NRPyPN: import numerical evaluation routine & gamma constant

def eval_random(i,trusted, other, p_t):
    random.seed(i)

    qmassratio = 1.0 + 7*random.random()  # must be >= 1
    nr         = 10. + 3*random.random()  # Orbital separation
    # Choose spins so that the total spin magnitude does not exceed 1.
    nchi1x     = -0.55 + 1.1*random.random()
    nchi1y     = -0.55 + 1.1*random.random()
    nchi1z     = -0.55 + 1.1*random.random()
    nchi2x     = -0.55 + 1.1*random.random()
    nchi2y     = -0.55 + 1.1*random.random()
    nchi2z     = -0.55 + 1.1*random.random()

    nPt        = num_eval(p_t,
                          qmassratio=qmassratio, nr=nr,
                          nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,
                          nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z)
    trusted_result = num_eval(trusted,qmassratio=qmassratio, nr=nr,
                              nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,
                              nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z,
                              nPt=nPt)
    other_result= num_eval(other,qmassratio=qmassratio, nr=nr,
                           nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,
                           nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z,
                           nPt=nPt)

    relerror = abs((trusted_result-other_result)/trusted_result)
    print("%04d" % (i+1), "|", #"%.2f" % int(count+1)/int(i+1),
          "%.4e" % (relerror),"|","%.4e" % trusted_result,"%.4e" % other_result,"|",
          "%.2f" % (qmassratio),"%.2f" % nr,
          "%.2f" % (nchi1x),"%.2f" % (nchi1y),"%.2f" % (nchi1z),
          "%.2f" % (nchi2x),"%.2f" % (nchi2y),"%.2f" % (nchi2z))
    return relerror

import PN_p_t as pt
pt.f_p_t(m1,m2, chi1U,chi2U, r)

f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
f_p_r(        m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)

num_trials = 1
relerror_array = []
if num_trials > 100:
    def master_func(i):
        return eval_random(i,p_r_RHP2018,p_r,pt.p_t)

    import multiprocessing
    # ixp.zerorank1(DIM=1000) #0.0
    pool = multiprocessing.Pool()
    relerror_array = pool.map(master_func,range(1000))
    pool.terminate()
    pool.join()
else:
    for i in range(num_trials):
        relerror_array.append(eval_random(i,p_r,p_r_fullHam,pt.p_t))

avg_relerror = 0
for i in range(len(relerror_array)):
    avg_relerror += relerror_array[i]
avg_relerror /= len(relerror_array)
print("Average relative error: ","%.5e" % (avg_relerror))

# Without 3.5PN term in p_r denominator (comparing
#      Eq 2.16 of https://arxiv.org/pdf/1810.00036.pdf with
#      Eq 16 of https://arxiv.org/pdf/1702.00872.pdf):
#  average relative error in p_r over 1000 random inputs = 8.24534e-03

# With 3.5PN term in p_r:
#  average relative error in p_r over 1000 random inputs = 8.19325e-03

0001 | 3.9019e-03 | 1.0678e-04 1.0720e-04 | 6.91 12.27 -0.09 -0.27 0.01 -0.10 0.31 -0.22
Average relative error:  3.90192e-03

Part 2.c: Validation against trusted numerical values (i.e., in Table V of Ramos-Buades, Husa, and Pratten (2018)) [Back to top]¶

$\label{code_validationpub}$

We will measure success in this validation to be within the error bar of the two values of $p_r$ given for each iteration, preferably being closer to the second iteration's value than the first.

In [10]:

import PN_p_t as pt
pt.f_p_t(m1,m2, chi1U,chi2U, r)
f_p_r_fullHam(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
f_p_r(        m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
import PN_p_r as pr
pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
p_r = pr.p_r

In [11]:

# Useful function for comparing published & NRPyPN results
def compare_pub_NPN(desc,pub0,pub1,  qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z):
    nPt        = num_eval(pt.p_t,
                          qmassratio=qmassratio, nr=nr,
                          nchi1x=nchi1x,nchi1y=nchi1y,nchi1z=nchi1z,
                          nchi2x=nchi2x,nchi2y=nchi2y,nchi2z=nchi2z)

    p_r_approach1 = num_eval(p_r_fullHam, qmassratio = qmassratio, nr=nr,
                             nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,
                             nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z,
                             nPt=nPt)

    p_r_approach2 = num_eval(p_r, qmassratio = qmassratio, nr=nr,
                             nchi1x=nchi1x, nchi1y=nchi1y, nchi1z=nchi1z,
                             nchi2x=nchi2x, nchi2y=nchi2y, nchi2z=nchi2z,
                             nPt=nPt)

    print("##################################################")
    print("  "+desc)
    print("##################################################")
    print("p_r = %.9f" % pub0, " <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)")
    print("p_r = %.9f" % p_r_approach1," <- Result from NRPyPN, approach 1")
    print("p_r = %.9f" % p_r_approach2," <- Result from NRPyPN, approach 2")
    print("p_r = %.9f" % pub1, " <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)")
    relerrorpct1  = abs((pub0-p_r_approach1)/pub0)*100
    relerrorpct2  = abs((pub0-p_r_approach2)/pub0)*100
    relerror01pct = abs((pub0-pub1)/pub1)*100
    strrelerror1pct  = "%.4f" % (relerrorpct1)
    strrelerror2pct  = "%.4f" % (relerrorpct2)
    strrelerror01pct = "%.4f" % (relerror01pct)
    print("Relative error between published iteration 0 vs it. 1: "+strrelerror01pct+"%")
    resultstring1 = "Relative error between NRPyPN & published, approach 1: "+strrelerror1pct+"%"
    resultstring2 = "Relative error between NRPyPN & published, approach 2: "+strrelerror2pct+"%"
    if relerrorpct1 > relerror01pct:
        resultstring1 += " < "+strrelerror01pct+"%  <--- NOT GOOD!"
        resultstring2 += " < "+strrelerror01pct+"%  <--- NOT GOOD!"
    else:
        resultstring1 += " < "+strrelerror01pct+"%  <--- GOOD AGREEMENT!"
        resultstring2 += " < "+strrelerror01pct+"%  <--- GOOD AGREEMENT!"
    print(resultstring1)
    print(resultstring2+"\n")

In [12]:

# 1. Let's consider the case:
# * Mass ratio q=1, chi1=chi2=(0,0,0), radial separation r=12
pub_result0 = 0.53833e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
pub_result1 = 0.468113e-3 # Expected result 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
qmassratio = 1.0  # must be >= 1
nr         = 12.  # Orbital separation
# Choose spins so that the total spin magnitude does not exceed 1.
nchi1x     = +0.
nchi1y     = +0.
nchi1z     = +0.
nchi2x     = +0.
nchi2y     = +0.
nchi2z     = +0.

compare_pub_NPN("Case: q=1, nonspinning, initial separation 12",pub_result0,pub_result1,
                qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)

##################################################
  Case: q=1, nonspinning, initial separation 12
##################################################
p_r = 0.000538330  <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)
p_r = 0.000539171  <- Result from NRPyPN, approach 1
p_r = 0.000539860  <- Result from NRPyPN, approach 2
p_r = 0.000468113  <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)
Relative error between published iteration 0 vs it. 1: 15.0000%
Relative error between NRPyPN & published, approach 1: 0.1562% < 15.0000%  <--- GOOD AGREEMENT!
Relative error between NRPyPN & published, approach 2: 0.2843% < 15.0000%  <--- GOOD AGREEMENT!

In [13]:

# 2. Let's consider the case:
# * Mass ratio q=1.5, chi1= (0,0,-0.6); chi2=(0,0,0.6), radial separation r=10.8
pub_result0 = 0.699185e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
pub_result1 = 0.641051e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
qmassratio =  1.5  # must be >= 1
nr         = 10.8  # Orbital separation
nchi1x     = +0.
nchi1y     = +0.
nchi1z     = -0.6
nchi2x     = +0.
nchi2y     = +0.
nchi2z     = +0.6
compare_pub_NPN("Case: q=1.5, chi1z=-0.6, chi2z=0.6, initial separation 10.8",pub_result0,pub_result1,
                qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)

##################################################
  Case: q=1.5, chi1z=-0.6, chi2z=0.6, initial separation 10.8
##################################################
p_r = 0.000699185  <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)
p_r = 0.000675944  <- Result from NRPyPN, approach 1
p_r = 0.000677020  <- Result from NRPyPN, approach 2
p_r = 0.000641051  <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)
Relative error between published iteration 0 vs it. 1: 9.0685%
Relative error between NRPyPN & published, approach 1: 3.3240% < 9.0685%  <--- GOOD AGREEMENT!
Relative error between NRPyPN & published, approach 2: 3.1702% < 9.0685%  <--- GOOD AGREEMENT!

In [14]:

# 3. Let's consider the case:
# * Mass ratio q=4, chi1= (0,0,-0.8); chi2=(0,0,0.8), radial separation r=11
pub_result0 = 0.336564e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
pub_result1 = 0.24708e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
qmassratio =  4.0  # must be >= 1
nr         = 11.0  # Orbital separation
nchi1x     = +0.
nchi1y     = +0.
nchi1z     = -0.8
nchi2x     = +0.
nchi2y     = +0.
nchi2z     = +0.8
compare_pub_NPN("Case: q=4.0, chi1z=-0.8, chi2z=0.8, initial separation 11.0",pub_result0,pub_result1,
                qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)

##################################################
  Case: q=4.0, chi1z=-0.8, chi2z=0.8, initial separation 11.0
##################################################
p_r = 0.000336564  <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)
p_r = 0.000251240  <- Result from NRPyPN, approach 1
p_r = 0.000251143  <- Result from NRPyPN, approach 2
p_r = 0.000247080  <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)
Relative error between published iteration 0 vs it. 1: 36.2166%
Relative error between NRPyPN & published, approach 1: 25.3515% < 36.2166%  <--- GOOD AGREEMENT!
Relative error between NRPyPN & published, approach 2: 25.3804% < 36.2166%  <--- GOOD AGREEMENT!

In [15]:

# 4. Let's consider the case:
# * Mass ratio q=2, chi1= (0,0,0); chi2=(−0.3535, 0.3535, 0.5), radial separation r=10.8
pub_result0 = 0.531374e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
pub_result1 = 0.448148e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
qmassratio =  2.0  # must be >= 1
nr         = 10.8  # Orbital separation
nchi1x     = +0.
nchi1y     = +0.
nchi1z     = +0.
nchi2x     = -0.3535
nchi2y     = +0.3535
nchi2z     = +0.5
compare_pub_NPN("Case: q=2.0, chi2x=-0.3535, chi2y=+0.3535, chi2z=+0.5, initial separation 10.8",pub_result0,pub_result1,
                qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)

##################################################
  Case: q=2.0, chi2x=-0.3535, chi2y=+0.3535, chi2z=+0.5, initial separation 10.8
##################################################
p_r = 0.000531374  <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)
p_r = 0.000545357  <- Result from NRPyPN, approach 1
p_r = 0.000542626  <- Result from NRPyPN, approach 2
p_r = 0.000448148  <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)
Relative error between published iteration 0 vs it. 1: 18.5711%
Relative error between NRPyPN & published, approach 1: 2.6314% < 18.5711%  <--- GOOD AGREEMENT!
Relative error between NRPyPN & published, approach 2: 2.1175% < 18.5711%  <--- GOOD AGREEMENT!

In [16]:

# 5. Let's consider the case:
# * Mass ratio q=8, chi1= (0, 0, 0.5); chi2=(0, 0, 0.5), radial separation r=11
pub_result0 = 0.102969e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
pub_result1 = 0.139132e-3 # Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018) https://arxiv.org/abs/1810.00036
qmassratio =  8.0  # must be >= 1
nr         = 11.0  # Orbital separation
nchi1x     = +0.
nchi1y     = +0.
nchi1z     = +0.5
nchi2x     = +0.
nchi2y     = +0.
nchi2z     = +0.5
compare_pub_NPN("""
Case: q=8.0, chi1z=chi2z=+0.5, initial separation 11

Note: This one is weird. Clearly the value in the table
      has a typo, such that the p_r and p_t values
      should be interchanged; p_t is about 20% the
      next smallest value in the table, and the
      parameters aren't that different. We therefore
      assume that this is the case, and nonetheless
      find agreement for p_t with the published result
      to about 0.07%, and p_r to about 5%. Aiven that
      the table values seem to be clearly wrong, this
      level of agreement is an encouraging sign.
""",pub_result0,pub_result1,
                qmassratio, nr, nchi1x,nchi1y,nchi1z, nchi2x,nchi2y,nchi2z)

##################################################
  
Case: q=8.0, chi1z=chi2z=+0.5, initial separation 11

Note: This one is weird. Clearly the value in the table
      has a typo, such that the p_r and p_t values
      should be interchanged; p_t is about 20% the
      next smallest value in the table, and the
      parameters aren't that different. We therefore
      assume that this is the case, and nonetheless
      find agreement for p_t with the published result
      to about 0.07%, and p_r to about 5%. Aiven that
      the table values seem to be clearly wrong, this
      level of agreement is an encouraging sign.

##################################################
p_r = 0.000102969  <- Expected result, from Table V of Ramos-Buades, Husa, and Pratten (2018)
p_r = 0.000097694  <- Result from NRPyPN, approach 1
p_r = 0.000097574  <- Result from NRPyPN, approach 2
p_r = 0.000139132  <- Result at 2nd iteration, from Table V of Ramos-Buades, Husa, and Pratten (2018)
Relative error between published iteration 0 vs it. 1: 25.9919%
Relative error between NRPyPN & published, approach 1: 5.1225% < 25.9919%  <--- GOOD AGREEMENT!
Relative error between NRPyPN & published, approach 2: 5.2399% < 25.9919%  <--- GOOD AGREEMENT!

Part 2.d: Validation against corresponding Python module [Back to top]¶

$\label{code_validationpython}$

Here we verify agreement between

the SymPy expressions generated in this notebook, and the corresponding Python module.

In [17]:

def error(varname):
    print("ERROR: When comparing Python module & notebook, "+varname+" was found not to match.")
    sys.exit(1)

# Validation against Python module
import PN_p_r as pr
f_p_r(   m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)
pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)

if sp.simplify(p_r - pr.p_r) != 0: error("pr.p_r")
print("PYTHON MODULE TEST PASSES")

PYTHON MODULE TEST PASSES

Part 3: Output this notebook to $\LaTeX$ -formatted PDF file [Back to top]¶

$\label{latex_pdf_output}$

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-p_r.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [18]:

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-p_r")

Created PN-p_r.tex, and compiled LaTeX file to PDF file PN-p_r.pdf

prp_r, the radial component of the momentum vector, up to and including 3.5 PN order¶