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 PN expressions. 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.
As described in the nonspinning Hamiltonian notebook, the basic physical system assumes two point particles of mass $m_1$ and $m_2$ with corresponding momentum vectors $\mathbf{P}_1$ and $\mathbf{P}_2$, and displacement vectors $\mathbf{X}_1$ and $\mathbf{X}_2$ with respect to the center of mass. Here we also consider the spin vectors of each point mass $\mathbf{S}_1$ and $\mathbf{S}_2$, respectively.
Damour, Jaranowski, and Schäfer (2008) adopt the notation
\begin{align} \mathbf{r}_{12} &= (\mathbf{X}_1-\mathbf{X}_2)\\ r_{12} = r_{21} &= |\mathbf{r}_{12}|\\ \mathbf{n}_{12} &= \frac{\mathbf{r}_{12}}{r_{12}}, \end{align}and when the numbers in subscripts are flipped, the particles are interchanged.
The spin-orbit terms of the Hamiltonian up to and including 3.5 PN order are generally given by:
$$ H_{\rm SO} = \mathbf{\Omega}_1 S^i_1 + \mathbf{\Omega}_2 S^i_2, $$where we need only define $\mathbf{\Omega}_1$ within a function, as $\mathbf{\Omega}_2$ is defined simply by interchanging $1\leftrightarrow 2$ in the $\mathbf{\Omega}_1$ expression.
At 1.5PN order (as summarized in Damour, Jaranowski, and Schäfer (2008), Eq 4.11a), we have
$$ \mathbf{\Omega}_{1,SO,1.5PN} = \frac{1}{r_{12}^2}\bigg( \frac{3m_2}{2m_1}{\mathbf{n}}_{12}\times{\mathbf{p}}_1 - 2 {\mathbf{n}}_{12}\times{\mathbf{p}}_2 \bigg), $$# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os, sys # Standard Python modules for multiplatform OS-level functions
import indexedexpNRPyPN as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
from NRPyPN_shortcuts import div,dot,cross # NRPyPN: shortcuts for e.g., vector operations
# 1.5PN spin-orbit coupling term, from Eq. 4.11a of
# Damour, Jaranowski, and Schäfer (2008)
# https://arxiv.org/abs/0711.1048
def f_H_SO_1p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
def f_Omega1(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (div(3,2)*m2/m1 * cross(n12U,p1U)[i] - 2*cross(n12U,p2U)[i])/r12**2
return Omega1
global H_SO_1p5PN
Omega1 = f_Omega1(m1,m2, n12U, p1U,p2U, r12)
Omega2 = f_Omega1(m2,m1, n21U, p2U,p1U, r12)
H_SO_1p5PN = dot(Omega1,S1U) + dot(Omega2,S2U)
# Second version, used for validation purposes only.
def f_H_SO_1p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
def f_Omega_SO_1p5PN(m1,m2, n12U, p1U,p2U, r12):
Omega1U = ixp.zerorank1()
for i in range(3):
Omega1U[i] = (div(3,2)*m2/m1 * cross(n12U,p1U)[i] - 2*cross(n12U,p2U)[i])/r12**2
return Omega1U
Omega1_1p5PNU = f_Omega_SO_1p5PN(m1,m2, n12U, p1U,p2U, r12)
Omega2_1p5PNU = f_Omega_SO_1p5PN(m2,m1, n21U, p2U,p1U, r12)
global H_SO_1p5PNv2
H_SO_1p5PNv2 = dot(Omega1_1p5PNU,S1U) + dot(Omega2_1p5PNU,S2U)
To reduce the possibility of copying errors, equations are taken directly from the arXiv LaTeX source code of Eq 4.11b in Damour, Jaranowski, and Schäfer (2008), 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 Hamiltonian notebooks:
\begin{align} \mathbf{\Omega}^{\rm NLO}_{1} &= \frac{G^2}{c^4r_{12}^3} \Bigg( \bigg(-\frac{11}{2}m_2-5\frac{m_2^2}{m_1}\bigg){\mathbf{n}}_{12}\times{\mathbf{p}}_1 + \bigg(6m_1+\frac{15}{2}m_2\bigg){\mathbf{n}}_{12}\times{\mathbf{p}}_2 \Bigg)\\ &\quad + \frac{G}{c^4r_{12}^2} \Bigg( \bigg( - \frac{5m_2{\bf p}_1^2}{8m_1^3} - \frac{3({\mathbf{p}}_1\cdot{\mathbf{p}}_2)}{4m_1^2} + \frac{3{\bf p}_2^2}{4m_1m_2} - \frac{3(\mathbf{n}_{12}\cdot\mathbf{p}_1)(\mathbf{n}_{12}\cdot\mathbf{p}_2)}{4m_1^2} - \frac{3(\mathbf{n}_{12}\cdot\mathbf{p}_2)^2}{2m_1m_2} \bigg){\mathbf{n}}_{12}\times{\mathbf{p}}_1 \\ &\quad\quad\quad\quad + \bigg(\frac{({\mathbf{p}}_1\cdot{\mathbf{p}}_2)}{m_1m_2}+\frac{3(\mathbf{n}_{12}\cdot\mathbf{p}_1)(\mathbf{n}_{12}\cdot\mathbf{p}_2)}{m_1m_2}\bigg){\mathbf{n}}_{12}\times{\mathbf{p}}_2 + \bigg( \frac{3(\mathbf{n}_{12}\cdot\mathbf{p}_1)}{4m_1^2} - \frac{2(\mathbf{n}_{12}\cdot\mathbf{p}_2)}{m_1m_2} \bigg){\mathbf{p}}_1\times{\mathbf{p}}_2 \Bigg). \end{align}# 2.5PN spin-orbit coupling term, from Eq. 4.11b of
# Damour, Jaranowski, and Schäfer (2008)
# https://arxiv.org/abs/0711.1048
def f_H_SO_2p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
def f_Omega_SO_2p5PN(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(+(-div(11,2)*m2-5*m2**2/m1)*cross(n12U,p1U)[i]
+(6*m1 + div(15,2)*m2) *cross(n12U,p2U)[i])/r12**3
+(+(-div(5,8)*m2*dot(p1U,p1U)/m1**3
-div(3,4)*dot(p1U,p2U)/m1**2
+div(3,4)*dot(p2U,p2U)/(m1*m2)
-div(3,4)*dot(n12U,p1U)*dot(n12U,p2U)/m1**2
-div(3,2)*dot(n12U,p2U)**2/(m1*m2))*cross(n12U,p1U)[i]
+(dot(p1U,p2U)/(m1*m2) + 3*dot(n12U,p1U)*dot(n12U,p2U)/(m1*m2))*cross(n12U,p2U)[i]
+(div(3,4)*dot(n12U,p1U)/m1**2 - 2*dot(n12U,p2U)/(m1*m2))*cross(p1U,p2U)[i])/r12**2)
return Omega1
Omega1_2p5PNU = f_Omega_SO_2p5PN(m1,m2, n12U, p1U,p2U, r12)
Omega2_2p5PNU = f_Omega_SO_2p5PN(m2,m1, n21U, p2U,p1U, r12)
global H_SO_2p5PN
H_SO_2p5PN = dot(Omega1_2p5PNU,S1U) + dot(Omega2_2p5PNU,S2U)
# Second version, used for validation purposes only.
def f_H_SO_2p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
def f_Omega_SO_2p5PNv2(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
n12_cross_p1 = cross(n12U,p1U)
n12_cross_p2 = cross(n12U,p2U)
for i in range(3):
Omega1[i] = ( (-div(11,2)*m2 - 5*m2**2/m1)*n12_cross_p1[i] + # line 1
(6*m1 + div(15,2)*m2) *n12_cross_p2[i] ) / r12**3 # line 1
Omega1[i]+= (( -div(5,8)*m2*dot(p1U,p1U)/m1**3 # line 2
-div(3,4)*dot(p1U,p2U)/m1**2 # line 2
+div(3,4)*dot(p2U,p2U)/(m1*m2) # line 2
-div(3,4)*dot(n12U,p1U)*dot(n12U,p2U)/m1**2 # line 2
-div(3,2)*dot(n12U,p2U)**2/(m1*m2) )*n12_cross_p1[i] + # line 2
( dot(p1U,p2U)/(m1*m2) + 3*dot(n12U,p1U)*dot(n12U,p2U)/(m1*m2) )*n12_cross_p2[i] + # line 3
(+div(3,4)*dot(n12U,p1U)/m1**2 - 2*dot(n12U,p2U)/(m1*m2) )*cross(p1U,p2U)[i] )/r12**2 # line 3
return Omega1
Omega1_2p5PNU = f_Omega_SO_2p5PNv2(m1,m2, n12U, p1U,p2U, r12)
Omega2_2p5PNU = f_Omega_SO_2p5PNv2(m2,m1, n21U, p2U,p1U, r12)
global H_SO_2p5PNv2
H_SO_2p5PNv2 = dot(Omega1_2p5PNU,S1U) + dot(Omega2_2p5PNU,S2U)
To reduce the possibility of copying errors, equations are taken directly from the arXiv LaTeX source code of Eq 5 in Hartung and Steinhoff (2011), 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 Hamiltonian notebooks:
\begin{align} H^{\text{NNLO}}_{\text{SO}} & = \frac{1}{r_{12}^2} \biggl[ \biggl( \frac{7 m_2 (\mathbf{P}_1^2)^2}{16 m_1^5} + \frac{9 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_1^2}{16 m_1^4} + \frac{3 \mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^3 m_2}\nonumber\\ &\quad\quad\quad + \frac{45 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)^3}{16 m_1^2 m_2^2} + \frac{9 \mathbf{P}_1^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{16 m_1^4} - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{16 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad - \frac{3 (\mathbf{P}_1^2) (\mathbf{P}_2^2)}{16 m_1^3 m_2} - \frac{15 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2) \mathbf{P}_2^2}{16 m_1^2 m_2^2} + \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2 \mathbf{P}_2^2}{4 m_1 m_2^3}\nonumber\\ &\quad\quad\quad - \frac{3 (\mathbf{P}_1\cdot\mathbf{P}_2) \mathbf{P}_2^2}{16 m_1^2 m_2^2} - \frac{3 (\mathbf{P}_2^2)^2}{16 m_1 m_2^3} \biggr)((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1)\\ &\quad\quad\quad +\biggl( - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_1^2}{2 m_1^3 m_2}\nonumber\\ &\quad\quad\quad - \frac{15 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2(\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^2 m_2^2} + \frac{3 \mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^2 m_2^2} - \frac{\mathbf{P}_1^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{2 m_1^3 m_2} + \frac{(\mathbf{P}_1\cdot\mathbf{P}_2)^2}{2 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad + \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2 \mathbf{P}_2^2}{4 m_1^2 m_2^2} - \frac{(\mathbf{P}_1^2) (\mathbf{P}_2^2)}{4 m_1^2 m_2^2} - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_2^2}{2 m_1 m_2^3}\nonumber\\ &\quad\quad\quad - \frac{(\mathbf{P}_1\cdot\mathbf{P}_2) \mathbf{P}_2^2}{2 m_1 m_2^3} \biggr)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1)\\ &\quad\quad\quad +\biggl( - \frac{9 (\mathbf{n}_{12}\cdot\mathbf{P}_1) \mathbf{P}_1^2}{16 m_1^4} + \frac{\mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)}{m_1^3 m_2}\nonumber\\ &\quad\quad\quad + \frac{27 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{16 m_1^2 m_2^2} - \frac{(\mathbf{n}_{12}\cdot\mathbf{P}_2)(\mathbf{P}_1\cdot\mathbf{P}_2)}{8 m_1^2 m_2^2} -\frac{5 (\mathbf{n}_{12}\cdot\mathbf{P}_1) \mathbf{P}_2^2}{16 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad + \frac{(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_2^2}{m_1 m_2^3} \biggr)((\mathbf{P}_1 \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \nonumber\\ &+ \frac{1}{r_{12}^3} \biggl[ \biggl( -\frac{3 m_2 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2}{2 m_1^2} +\left( -\frac{3 m_2}{2 m_1^2} +\frac{27 m_2^2}{8 m_1^3} \right) \mathbf{P}_1^2 +\left( \frac{177}{16 m_1} +\frac{11}{m_2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2\nonumber\\ &\quad\quad\quad +\left( \frac{11}{2 m_1} +\frac{9 m_2}{2 m_1^2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_1) (\mathbf{n}_{12}\cdot\mathbf{P}_2) +\left( \frac{23}{4 m_1} +\frac{9 m_2}{2 m_1^2} \right) (\mathbf{P}_1\cdot\mathbf{P}_2)\nonumber\\ &\quad\quad\quad -\left( \frac{159}{16 m_1} +\frac{37}{8 m_2} \right) \mathbf{P}_2^2 \biggr)((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1) +\biggl( \frac{4 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2}{m_1} +\frac{13 \mathbf{P}_1^2}{2 m_1}\nonumber\\ &\quad\quad\quad +\frac{5 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{m_2} +\frac{53 \mathbf{P}_2^2}{8 m_2} - \left( \frac{211}{8 m_1} +\frac{22}{m_2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_1) (\mathbf{n}_{12}\cdot\mathbf{P}_2)\nonumber\\ &\quad\quad\quad -\left( \frac{47}{8 m_1} +\frac{5}{m_2} \right)(\mathbf{P}_1\cdot\mathbf{P}_2) \biggr)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1) +\biggl( -\left( \frac{8}{m_1} +\frac{9 m_2}{2 m_1^2} \right)(\mathbf{n}_{12}\cdot\mathbf{P}_1)\nonumber\\ &\quad\quad\quad +\left( \frac{59}{4 m_1} +\frac{27}{2 m_2} \right)(\mathbf{n}_{12}\cdot\mathbf{P}_2) \biggr)((\mathbf{P}_1 \times \mathbf{P}_2)\mathbf{S}_1) \biggr]\nonumber\\ &+\frac{1}{r_{12}^4} \biggl[ \left( \frac{181 m_1 m_2}{16} + \frac{95 m_2^2}{4} + \frac{75 m_2^3}{8 m_1} \right) ((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1)\nonumber\\ &\quad\quad\quad - \left( \frac{21 m_1^2}{2} + \frac{473 m_1 m_2}{16} + \frac{63 m_2^2}{4} \right)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1) \biggr] + (1\leftrightarrow2)\,. \end{align}Let's split the above into more bite-sized pieces. First:
\begin{align} H^a_{SO,2.5PN} &= \frac{1}{r_{12}^2} \biggl[ \biggl( \frac{7 m_2 (\mathbf{P}_1^2)^2}{16 m_1^5} + \frac{9 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_1^2}{16 m_1^4} + \frac{3 \mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^3 m_2}\nonumber\\ &\quad\quad\quad + \frac{45 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)^3}{16 m_1^2 m_2^2} + \frac{9 \mathbf{P}_1^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{16 m_1^4} - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{16 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad - \frac{3 (\mathbf{P}_1^2) (\mathbf{P}_2^2)}{16 m_1^3 m_2} - \frac{15 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2) \mathbf{P}_2^2}{16 m_1^2 m_2^2} + \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2 \mathbf{P}_2^2}{4 m_1 m_2^3}\nonumber\\ &\quad\quad\quad - \frac{3 (\mathbf{P}_1\cdot\mathbf{P}_2) \mathbf{P}_2^2}{16 m_1^2 m_2^2} - \frac{3 (\mathbf{P}_2^2)^2}{16 m_1 m_2^3} \biggr)((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1)\biggr] \end{align}# 3.5PN spin-orbit coupling term, from Eq. 5 of
# Hartung and Steinhoff (2011)
# https://arxiv.org/abs/1104.3079
# 3.5PN H_SO: Omega_1, part 1:
def HS2011_Omega_SO_3p5PN_pt1(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ((+7*m2*dot(p1U,p1U)**2/(16*m1**5)
+9*dot(n12U,p1U)*dot(n12U,p2U)*dot(p1U,p1U)/(16*m1**4)
+3*dot(p1U,p1U)*dot(n12U,p2U)**2/(4*m1**3*m2)
+45*dot(n12U,p1U)*dot(n12U,p2U)**3/(16*m1**2*m2**2)
+9*dot(p1U,p1U)*dot(p1U,p2U)/(16*m1**4)
-3*dot(n12U,p2U)**2*dot(p1U,p2U)/(16*m1**2*m2**2)
-3*dot(p1U,p1U)*dot(p2U,p2U)/(16*m1**3*m2)
-15*dot(n12U,p1U)*dot(n12U,p2U)*dot(p2U,p2U)/(16*m1**2*m2**2)
+3*dot(n12U,p2U)**2*dot(p2U,p2U)/(4*m1*m2**3)
-3*dot(p1U,p2U)*dot(p2U,p2U)/(16*m1**2*m2**2)
-3*dot(p2U,p2U)**2/(16*m1*m2**3))*cross(n12U,p1U)[i])/r12**2
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt1v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (+div(7,16)*m2*dot(p1U,p1U)**2/m1**5
+div(9,16)*dot(n12U,p1U)*dot(n12U,p2U)*dot(p1U,p1U)/m1**4
+div(3,4) *dot(p1U,p1U)*dot(n12U,p2U)**2/(m1**3*m2)
+div(45,16)*dot(n12U,p1U)*dot(n12U,p2U)**3/(m1**2*m2**2)
+div(9,16)*dot(p1U,p1U)*dot(p1U,p2U)/m1**4
-div(3,16)*dot(n12U,p2U)**2*dot(p1U,p2U)/(m1**2*m2**2)
-div(3,16)*dot(p1U,p1U)*dot(p2U,p2U)/(m1**3*m2)
-div(15,16)*dot(n12U,p1U)*dot(n12U,p2U)*dot(p2U,p2U)/(m1**2*m2**2)
+div(3,4)*dot(n12U,p2U)**2*dot(p2U,p2U)/(m1*m2**3)
-div(3,16)*dot(p1U,p2U)*dot(p2U,p2U)/(m1**2*m2**2)
-div(3,16)*dot(p2U,p2U)**2/(m1*m2**3)) * cross(n12U,p1U)[i] )/q**2
return Omega1
Next,
\begin{align} H^b_{SO,2.5PN} &= \frac{1}{r_{12}^2} \biggl[ +\biggl( - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_1^2}{2 m_1^3 m_2}\nonumber\\ &\quad\quad\quad - \frac{15 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2(\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^2 m_2^2} + \frac{3 \mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{4 m_1^2 m_2^2} - \frac{\mathbf{P}_1^2 (\mathbf{P}_1\cdot\mathbf{P}_2)}{2 m_1^3 m_2} + \frac{(\mathbf{P}_1\cdot\mathbf{P}_2)^2}{2 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad + \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2 \mathbf{P}_2^2}{4 m_1^2 m_2^2} - \frac{(\mathbf{P}_1^2) (\mathbf{P}_2^2)}{4 m_1^2 m_2^2} - \frac{3 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_2^2}{2 m_1 m_2^3}\nonumber\\ &\quad\quad\quad - \frac{(\mathbf{P}_1\cdot\mathbf{P}_2) \mathbf{P}_2^2}{2 m_1 m_2^3} \biggr)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 2:
def HS2011_Omega_SO_3p5PN_pt2(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(-3*dot(n12U,p1U)*dot(n12U,p2U)*dot(p1U,p1U)/(2*m1**3*m2)
-15*dot(n12U,p1U)**2*dot(n12U,p2U)**2/(4*m1**2*m2**2)
+3*dot(p1U,p1U)*dot(n12U,p2U)**2/(4*m1**2*m2**2)
-dot(p1U,p1U)*dot(p1U,p2U)/(2*m1**3*m2)
+dot(p1U,p2U)**2/(2*m1**2*m2**2)
+3*dot(n12U,p1U)**2*dot(p2U,p2U)/(4*m1**2*m2**2)
-dot(p1U,p1U)*dot(p2U,p2U)/(4*m1**2*m2**2)
-3*dot(n12U,p1U)*dot(n12U,p2U)*dot(p2U,p2U)/(2*m1*m2**3)
-dot(p1U,p2U)*dot(p2U,p2U)/(2*m1*m2**3))*cross(n12U,p2U)[i])/r12**2
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt2v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (-div(3,2)*dot(n12U,p1U)*dot(n12U,p2U)*dot(p1U,p1U)/(m1**3*m2)
-div(15,4)*dot(n12U,p1U)**2*dot(n12U,p2U)**2/(m1**2*m2**2)
+div(3,4)*dot(p1U,p1U)*dot(n12U,p2U)**2/(m1**2*m2**2)
-div(1,2)*dot(p1U,p1U)*dot(p1U,p2U)/(m1**3*m2)
+div(1,2)*dot(p1U,p2U)**2/(m1**2*m2**2)
+div(3,4)*dot(n12U,p1U)**2*dot(p2U,p2U)/(m1**2*m2**2)
-div(1,4)*dot(p1U,p1U)*dot(p2U,p2U)/(m1**2*m2**2)
-div(3,2)*dot(n12U,p1U)*dot(n12U,p2U)*dot(p2U,p2U)/(m1*m2**3)
-div(1,2)*dot(p1U,p2U)*dot(p2U,p2U)/(m1*m2**3))*cross(n12U,p2U)[i] )/q**2
return Omega1
Part 3:
\begin{align} H^c_{SO,2.5PN} &= \frac{1}{r_{12}^2} \biggl[ +\biggl( - \frac{9 (\mathbf{n}_{12}\cdot\mathbf{P}_1) \mathbf{P}_1^2}{16 m_1^4} + \frac{\mathbf{P}_1^2 (\mathbf{n}_{12}\cdot\mathbf{P}_2)}{m_1^3 m_2}\nonumber\\ &\quad\quad\quad + \frac{27 (\mathbf{n}_{12}\cdot\mathbf{P}_1)(\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{16 m_1^2 m_2^2} - \frac{(\mathbf{n}_{12}\cdot\mathbf{P}_2)(\mathbf{P}_1\cdot\mathbf{P}_2)}{8 m_1^2 m_2^2} -\frac{5 (\mathbf{n}_{12}\cdot\mathbf{P}_1) \mathbf{P}_2^2}{16 m_1^2 m_2^2}\nonumber\\ &\quad\quad\quad + \frac{(\mathbf{n}_{12}\cdot\mathbf{P}_2)\mathbf{P}_2^2}{m_1 m_2^3} \biggr)((\mathbf{P}_1 \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 3:
def HS2011_Omega_SO_3p5PN_pt3(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(-9*dot(n12U,p1U)*dot(p1U,p1U)/(16*m1**4)
+dot(p1U,p1U)*dot(n12U,p2U)/(m1**3*m2)
+27*dot(n12U,p1U)*dot(n12U,p2U)**2/(16*m1**2*m2**2)
-dot(n12U,p2U)*dot(p1U,p2U)/(8*m1**2*m2**2)
-5*dot(n12U,p1U)*dot(p2U,p2U)/(16*m1**2*m2**2))*cross(p1U,p2U)[i])/r12**2
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt3v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (-div(9,16)*dot(n12U,p1U)*dot(p1U,p1U)/m1**4
+ dot(p1U,p1U)*dot(n12U,p2U)/(m1**3*m2)
+div(27,16)*dot(n12U,p1U)*dot(n12U,p2U)**2/(m1**2*m2**2)
-div(1,8)*dot(n12U,p2U)*dot(p1U,p2U)/(m1**2*m2**2)
-div(5,16)*dot(n12U,p1U)*dot(p2U,p2U)/(m1**2*m2**2)
+ dot(n12U,p2U)*dot(p2U,p2U)/(m1*m2**3))*cross(p1U,p2U)[i] )/q**2
return Omega1
Part 4, the first $1/r_{12}^3$ term:
\begin{align} H^d_{SO,2.5PN} &= \frac{1}{r_{12}^3} \biggl[ \biggl( -\frac{3 m_2 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2}{2 m_1^2} +\left( -\frac{3 m_2}{2 m_1^2} +\frac{27 m_2^2}{8 m_1^3} \right) \mathbf{P}_1^2 +\left( \frac{177}{16 m_1} +\frac{11}{m_2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2\nonumber\\ &\quad\quad\quad +\left( \frac{11}{2 m_1} +\frac{9 m_2}{2 m_1^2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_1) (\mathbf{n}_{12}\cdot\mathbf{P}_2) +\left( \frac{23}{4 m_1} +\frac{9 m_2}{2 m_1^2} \right) (\mathbf{P}_1\cdot\mathbf{P}_2)\nonumber\\ &\quad\quad\quad -\left( \frac{159}{16 m_1} +\frac{37}{8 m_2} \right) \mathbf{P}_2^2 \biggr)((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 4:
def HS2011_Omega_SO_3p5PN_pt4(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(-3*m2*dot(n12U,p1U)**2/(2*m1**2)
+((-3*m2)/(2*m1**2) + 27*m2**2/(8*m1**3))*dot(p1U,p1U)
+(177/(16*m1) + 11/m2)*dot(n12U,p2U)**2
+(11/(2*m1) + 9*m2/(2*m1**2))*dot(n12U,p1U)*dot(n12U,p2U)
+(23/(4*m1) + 9*m2/(2*m1**2))*dot(p1U,p2U)
-(159/(16*m1) + 37/(8*m2))*dot(p2U,p2U))*cross(n12U,p1U)[i])/r12**3
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt4v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (-div(3,2)*m2*dot(n12U,p1U)**2/m1**2
+(-div(3,2)*m2/m1**2 + div(27,8)*m2**2/m1**3)*dot(p1U,p1U)
+(+div(177,16)/m1 + 11/m2)*dot(n12U,p2U)**2
+(+div(11,2)/m1 + div(9,2)*m2/m1**2)*dot(n12U,p1U)*dot(n12U,p2U)
+(+div(23,4)/m1 + div(9,2)*m2/m1**2)*dot(p1U,p2U)
-(+div(159,16)/m1 + div(37,8)/m2)*dot(p2U,p2U) )*cross(n12U,p1U)[i] )/q**3
return Omega1
Part 5, the second $1/r_{12}^3$ term:
\begin{align} H^e_{SO,2.5PN} &= \frac{1}{r_{12}^3} \biggl[ +\biggl( \frac{4 (\mathbf{n}_{12}\cdot\mathbf{P}_1)^2}{m_1} +\frac{13 \mathbf{P}_1^2}{2 m_1}\nonumber\\ &\quad\quad\quad +\frac{5 (\mathbf{n}_{12}\cdot\mathbf{P}_2)^2}{m_2} +\frac{53 \mathbf{P}_2^2}{8 m_2} - \left( \frac{211}{8 m_1} +\frac{22}{m_2} \right) (\mathbf{n}_{12}\cdot\mathbf{P}_1) (\mathbf{n}_{12}\cdot\mathbf{P}_2)\nonumber\\ &\quad\quad\quad -\left( \frac{47}{8 m_1} +\frac{5}{m_2} \right)(\mathbf{P}_1\cdot\mathbf{P}_2) \biggr)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 5:
def HS2011_Omega_SO_3p5PN_pt5(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(+4*dot(n12U,p1U)**2/m1
+13*dot(p1U,p1U)/(2*m1)
+5*dot(n12U,p2U)**2/m2
+53*dot(p2U,p2U)/(8*m2)
-(211/(8*m1) + 22/m2)*dot(n12U,p1U)*dot(n12U,p2U)
-(47/(8*m1) + 5/m2)*dot(p1U,p2U))*cross(n12U,p2U)[i])/r12**3
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt5v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (+4*dot(n12U,p1U)**2/m1
+div(13,2)*dot(p1U,p1U)/m1
+5*dot(n12U,p2U)**2/m2
+div(53,8)*dot(p2U,p2U)/m2
-(div(211,8)/m1+22/m2)*dot(n12U,p1U)*dot(n12U,p2U)
-(div(47,8)/m1+5/m2)*dot(p1U,p2U)) * cross(n12U,p2U)[i] )/q**3
return Omega1
Part 6, the third $1/r_{12}^3$ term:
\begin{align} H^f_{SO,2.5PN} &= \frac{1}{r_{12}^3} \biggl[ +\biggl( -\left( \frac{8}{m_1} +\frac{9 m_2}{2 m_1^2} \right)(\mathbf{n}_{12}\cdot\mathbf{P}_1)\nonumber\\ &\quad\quad\quad +\left( \frac{59}{4 m_1} +\frac{27}{2 m_2} \right)(\mathbf{n}_{12}\cdot\mathbf{P}_2) \biggr)((\mathbf{P}_1 \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 6:
def HS2011_Omega_SO_3p5PN_pt6(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(-(8/m1 + 9*m2/(2*m1**2))*dot(n12U,p1U)
+(59/(4*m1) + 27/(2*m2))*dot(n12U,p2U))*cross(p1U,p2U)[i])/r12**3
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt6v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( (-( 8/m1 + div(9,2)*m2/m1**2)*dot(n12U,p1U)
+(div(59,4)/m1 + div(27,2)/m2) *dot(n12U,p2U))*cross(p1U,p2U)[i] )/q**3
return Omega1
Finally part 7, the $1/r_{12}^4$ term:
\begin{align} H^f_{SO,2.5PN} &= \frac{1}{r_{12}^4} \biggl[ \left( \frac{181 m_1 m_2}{16} + \frac{95 m_2^2}{4} + \frac{75 m_2^3}{8 m_1} \right) ((\mathbf{n}_{12} \times \mathbf{P}_1)\mathbf{S}_1)\nonumber\\ &\quad\quad\quad - \left( \frac{21 m_1^2}{2} + \frac{473 m_1 m_2}{16} + \frac{63 m_2^2}{4} \right)((\mathbf{n}_{12} \times \mathbf{P}_2)\mathbf{S}_1) \biggr] \end{align}# 3.5PN H_SO: Omega_1, part 7:
def HS2011_Omega_SO_3p5PN_pt7(m1,m2, n12U, p1U,p2U, r12):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = (+(181*m1*m2/16 + 95*m2**2/4 + 75*m2**3/(8*m1))*cross(n12U,p1U)[i]
-(21*m1**2/2 + 473*m1*m2/16 + 63*m2**2/4)*cross(n12U,p2U)[i])/r12**4
return Omega1
# Second version, used for validation purposes only.
def HS2011_Omega_SO_3p5PN_pt7v2(m1,m2, n12U, p1U,p2U, q):
Omega1 = ixp.zerorank1()
for i in range(3):
Omega1[i] = ( +(div(181,16)*m1*m2 + div(95,4)*m2**2 + div(75,8)*m2**3/m1)*cross(n12U,p1U)[i]
-(div(21,2)*m1**2 + div(473,16)*m1*m2 + div(63,4)*m2**2 )*cross(n12U,p2U)[i] )/q**4
return Omega1
Now we put all the $\Omega$ terms together:
def f_H_SO_3p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
Omega1_3p5PNU = ixp.zerorank1()
Omega2_3p5PNU = ixp.zerorank1()
for i in range(3):
Omega1_3p5PNU[i] = HS2011_Omega_SO_3p5PN_pt1(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt3(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt4(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt5(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt6(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt7(m1,m2, n12U, p1U,p2U, r12)[i]
Omega2_3p5PNU[i] = HS2011_Omega_SO_3p5PN_pt1(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt3(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt4(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt5(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt6(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt7(m2,m1, n21U, p2U,p1U, r12)[i]
global H_SO_3p5PN
H_SO_3p5PN = dot(Omega1_3p5PNU,S1U) + dot(Omega2_3p5PNU,S2U)
# For validation purposes only:
def f_H_SO_3p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, r12):
Omega1_3p5PNU = ixp.zerorank1()
Omega2_3p5PNU = ixp.zerorank1()
for i in range(3):
Omega1_3p5PNU[i] = HS2011_Omega_SO_3p5PN_pt1v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt2v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt3v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt4v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt5v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt6v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega1_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt7v2(m1,m2, n12U, p1U,p2U, r12)[i]
Omega2_3p5PNU[i] = HS2011_Omega_SO_3p5PN_pt1v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt2v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt3v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt4v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt5v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt6v2(m2,m1, n21U, p2U,p1U, r12)[i]
Omega2_3p5PNU[i]+= HS2011_Omega_SO_3p5PN_pt7v2(m2,m1, n21U, p2U,p1U, r12)[i]
global H_SO_3p5PNv2
H_SO_3p5PNv2 = dot(Omega1_3p5PNU,S1U) + dot(Omega2_3p5PNU,S2U)
As a code validation check, we verify agreement between
from NRPyPN_shortcuts import m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q # Import needed input variables
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
f_H_SO_1p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
f_H_SO_2p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
f_H_SO_3p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
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_H_SO_1p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
f_H_SO_2p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
f_H_SO_3p5PNv2(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
if sp.simplify(H_SO_1p5PN - H_SO_1p5PNv2) != 0: error("H_SO_1p5PNv2")
if sp.simplify(H_SO_2p5PN - H_SO_2p5PNv2) != 0: error("H_SO_2p5PNv2")
if sp.simplify(H_SO_3p5PN - H_SO_3p5PNv2) != 0: error("H_SO_3p5PNv2")
# Validation against corresponding Python module:
import PN_Hamiltonian_SO as HSO
HSO.f_H_SO_1p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
HSO.f_H_SO_2p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
HSO.f_H_SO_3p5PN(m1,m2, n12U,n21U, S1U, S2U, p1U,p2U, q)
if sp.simplify(H_SO_1p5PN - HSO.H_SO_1p5PN) != 0: error("H_SO_1p5PN")
if sp.simplify(H_SO_2p5PN - HSO.H_SO_2p5PN) != 0: error("H_SO_2p5PN")
if sp.simplify(H_SO_3p5PN - HSO.H_SO_3p5PN) != 0: error("H_SO_3p5PN")
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-Hamiltonian-Spin-Orbit.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helperNRPyPN as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("PN-Hamiltonian-Spin-Orbit")
Created PN-Hamiltonian-Spin-Orbit.tex, and compiled LaTeX file to PDF file PN-Hamiltonian-Spin-Orbit.pdf