NRPyPN: Validated Post-Newtonian Expressions for Input into Wolfram Mathematica, SymPy, or Highly-Optimized C Codes. Also provides low-eccentricity binary black hole initial data parameters for Bowen-York initial data (e.g., TwoPunctures)

Lead author: Zachariah B. Etienne $\leftarrow$ Please feel free to email contributions, comments, revisions, or errata!

Introduction

Post-Newtonian theory results in some of the longest and most complex mathematical expressions ever derived by humanity.

These expressions form the core of most gravitational wave data analysis codes, but generally the expressions are written in a format that is either inaccessible by others (i.e., closed-source) or written directly in e.g., C code and thus incompatible with symbolic algebra packages like Wolfram Mathematica or the Python-based SymPy.

Once in a symbolic algebra package, these expressions could be manipulated, extended, and output as more optimized C codes (e.g., using SymPy/NRPy+, thus speeding up gravitational wave data analysis).

This repository aims to provide a trusted source for validated Post-Newtonian expressions useful for gravitational wave astronomy, using the open-source SymPy computer algebra system, so that expressions can be output into Wolfram Mathematica or highly optimized C codes using the SymPy-based NRPy+.

If you are unfamiliar with using Jupyter Notebooks, first review the official Jupyter Notebook Basics Guide.

PART 0: Basic Functions for Expediting Equation Imports

PART 1: Post-Newtonian Hamiltonian for Binary Black Holes

The Post-Newtonian Hamiltonian $H$ can be written as a sum of contributions: \begin{array} \ H = H_{\rm Orb, NS} + H_{\rm SO} + H_{\rm SS} + H_{\rm SSS}, \end{array} where

  • $H_{\rm Orb, NS}$ is the non-spinning, purely orbital contribution
  • $H_{\rm SO}$ accounts for spin-orbit coupling
  • $H_{\rm SS}$ accounts for spin-spin contributions
  • $H_{\rm SSS}$ accounts for spin-spin-spin contributions

Click on any term below to open a notebook implementing that Hamiltonian component.

$H_{\rm Orb, NS}$, up to and including third post-Newtonian order (3PN)

$H_{\rm SO}$, up to and including 3.5 post-Newtonian order (3.5PN)

$H_{\rm SS}$, up to and including third post-Newtonian order (3PN)

$H_{\rm SSS}$, up to and including third post-Newtonian order (3PN)

Part 2: $\frac{dE_{\rm GW}}{dt}$, the gravitational wave flux, and $\frac{dM}{dt}$, the tidal energy injected into the black holes

Part 3: Quasicircular Orbital Parameters $p_t$ and $p_r$

Fast $p_t$ and $p_r$ for binary black hole initial data

Assumptions:

The following code cell assumes the black holes are

  • situated initially on the $x$-axis, with black hole 1 at $x>0$, and black hole 2 at $x<0$
  • orbiting instantaneously in the $x$-$y$ plane, with the center-of-mass at $x=y=z=0$.

Note that other initial coordinate configurations are possible; just permute the coordinates consistently with a right-handed coordinate system. For example,

  • perform the conversion $(x,y,z)\to(z,x,y)$ if one desires the black holes to be on the $z-axis$ instantaneously orbiting on the $z$-$x$ plane; or
  • perform the conversion $(x,y,z)\to(y,z,x)$ if one desires the black holes to be on the $y-axis$ instantaneously orbiting on the $y$-$z$ plane.

Instructions: Set the 8 parameters for the binary black hole system in the code cell below, which make the above assumptions. Then after setting the parameters press Shift+Enter. The tangential and radial components of the orbital momentum vector will be computed and output.

If using TwoPunctures: If you plan to use these parameters with TwoPunctures, you will need to set

  • TwoPunctures::give_bare_mass = no
  • TwoPunctures::target_m_minus = $q/(1+q)$
  • TwoPunctures::target_m_plus = $1/(1+q)$
  • TwoPunctures::par_P_plus[0] = -|P_r| (from code cell output below)
  • TwoPunctures::par_P_plus[1] = +|P_t| (from code cell output below)
  • TwoPunctures::par_P_minus[0] = +|P_r| (from code cell output below)
  • TwoPunctures::par_P_minus[1] = -|P_t| (from code cell output below)
  • TwoPunctures::par_S_plus[0] = nchi1x*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus (from input into code cell below)
  • TwoPunctures::par_S_plus[1] = nchi1y*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus (from input into code cell below)
  • TwoPunctures::par_S_plus[2] = nchi1z*TwoPunctures::target_m_plus*TwoPunctures::target_m_plus (from input into code cell below)
  • similarly for TwoPunctures::par_S_minus[]
In [1]:
#### INPUT PARAMETERS:
qmassratio =  1.75  # m2/m1; by convention must be >= 1
nr         =  10.8  # Orbital separation
# Dimensionless spin parameters of each black hole, assuming
#   the black holes are initially orbiting on the xy-plane
#   and situated on the x-axis, with the center-of-mass at
#   x=y=z=0. Black hole 1 is at x>0; and black hole 2 at x<0.
nchi1x     = +0.   # chi_{1x}, x-component of dimensionless spin for black hole 1 (instantaneously at x>0)
nchi1y     = +0.   # chi_{1y}, y-component of dimensionless spin for black hole 1 (instantaneously at x>0)
nchi1z     = +0.6  # chi_{1z}, z-component of dimensionless spin for black hole 1 (instantaneously at x>0)
nchi2x     = +0.   # chi_{2x}, x-component of dimensionless spin for black hole 2 (instantaneously at x<0)
nchi2y     = +0.   # chi_{2y}, y-component of dimensionless spin for black hole 2 (instantaneously at x<0)
nchi2z     = +0.6  # chi_{2z}, z-component of dimensionless spin for black hole 2 (instantaneously at x<0)

##########################################################
#### DON'T TOUCH; see output after running this cell. ####
from NRPyPN_shortcuts import eval__P_t__and__P_r

# Compute nPt, the tangential component of orbital momentum vector (output = numerical), as well as
#         nPr, the radial component of orbital momentum vector (output = numerical)
nPt, nPr = eval__P_t__and__P_r(qmassratio, nr,
                               nchi1x, nchi1y, nchi1z,
                               nchi2x, nchi2y, nchi2z)

print("If the black holes are situated on the x-axis, and orbiting instantaneously in the")
print("   xy-plane, then one could have P_y=+|P_t| and P_x=-|P_r| for the x>0 black hole,")
print("   and P_y=-|P_t| and P_x=+|P_r| for the x<0 black hole.  Other configurations are")
print("   possible provided one permutes the coordinates consistently with a right-handed")
print("   coordinate system.")
print("|P_t| = ",nPt)
print("|P_r| = ",nPr)
If the black holes are situated on the x-axis, and orbiting instantaneously in the
   xy-plane, then one could have P_y=+|P_t| and P_x=-|P_r| for the x>0 black hole,
   and P_y=-|P_t| and P_x=+|P_r| for the x<0 black hole.  Other configurations are
   possible provided one permutes the coordinates consistently with a right-handed
   coordinate system.
|P_t| =  0.0813467633310249
|P_r| =  0.000553891860334773